src/gromacs/gmxana/gmx_disre.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <math.h>
  40 #include <stdlib.h>
  41 #include <string.h>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/matio.h"
  46 #include "gromacs/fileio/pdbio.h"
  47 #include "gromacs/fileio/tpxio.h"
  48 #include "gromacs/fileio/trxio.h"
  49 #include "gromacs/fileio/xvgr.h"
  50 #include "gromacs/gmxana/gmx_ana.h"
  51 #include "gromacs/gmxana/gstat.h"
  52 #include "gromacs/legacyheaders/disre.h"
  53 #include "gromacs/legacyheaders/force.h"
  54 #include "gromacs/legacyheaders/macros.h"
  55 #include "gromacs/legacyheaders/main.h"
  56 #include "gromacs/legacyheaders/mdatoms.h"
  57 #include "gromacs/legacyheaders/mdrun.h"
  58 #include "gromacs/legacyheaders/names.h"
  59 #include "gromacs/legacyheaders/nrnb.h"
  60 #include "gromacs/legacyheaders/typedefs.h"
  61 #include "gromacs/legacyheaders/viewit.h"
  62 #include "gromacs/math/do_fit.h"
  63 #include "gromacs/math/vec.h"
  64 #include "gromacs/pbcutil/ishift.h"
  65 #include "gromacs/pbcutil/mshift.h"
  66 #include "gromacs/pbcutil/pbc.h"
  67 #include "gromacs/pbcutil/rmpbc.h"
  68 #include "gromacs/topology/index.h"
  69 #include "gromacs/topology/mtop_util.h"
  70 #include "gromacs/utility/fatalerror.h"
  71 #include "gromacs/utility/smalloc.h"
  72
  73 typedef struct {
  74     int  n;
  75     real v;
  76 } t_toppop;
  77
  78 t_toppop *top  = NULL;
  79 int       ntop = 0;
  80
  81 typedef struct {
  82     int   nv, nframes;
  83     real  sumv, averv, maxv;
  84     real *aver1, *aver2, *aver_3, *aver_6;
  85 } t_dr_result;
  86
  87 static void init5(int n)
  88 {
  89     ntop = n;
  90     snew(top, ntop);
  91 }
  92
  93 static void reset5(void)
  94 {
  95     int i;
  96
  97     for (i = 0; (i < ntop); i++)
  98     {
  99         top[i].n = -1;
 100         top[i].v = 0;
 101     }
 102 }
 103
 104 int tpcomp(const void *a, const void *b)
 105 {
 106     t_toppop *tpa;
 107     t_toppop *tpb;
 108
 109     tpa = (t_toppop *)a;
 110     tpb = (t_toppop *)b;
 111
 112     return  (1e7*(tpb->v-tpa->v));
 113 }
 114
 115 static void add5(int ndr, real viol)
 116 {
 117     int i, mini;
 118
 119     mini = 0;
 120     for (i = 1; (i < ntop); i++)
 121     {
 122         if (top[i].v < top[mini].v)
 123         {
 124             mini = i;
 125         }
 126     }
 127     if (viol > top[mini].v)
 128     {
 129         top[mini].v = viol;
 130         top[mini].n = ndr;
 131     }
 132 }
 133
 134 static void print5(FILE *fp)
 135 {
 136     int i;
 137
 138     qsort(top, ntop, sizeof(top[0]), tpcomp);
 139     fprintf(fp, "Index:");
 140     for (i = 0; (i < ntop); i++)
 141     {
 142         fprintf(fp, " %6d", top[i].n);
 143     }
 144     fprintf(fp, "\nViol: ");
 145     for (i = 0; (i < ntop); i++)
 146     {
 147         fprintf(fp, " %6.3f", top[i].v);
 148     }
 149     fprintf(fp, "\n");
 150 }
 151
 152 static void check_viol(FILE *log,
 153                        t_ilist *disres, t_iparams forceparams[],
 154                        rvec x[], rvec f[],
 155                        t_pbc *pbc, t_graph *g, t_dr_result dr[],
 156                        int clust_id, int isize, atom_id index[], real vvindex[],
 157                        t_fcdata *fcd)
 158 {
 159     t_iatom         *forceatoms;
 160     int              i, j, nat, n, type, nviol, ndr, label;
 161     real             ener, rt, mviol, tviol, viol, lam, dvdl, drt;
 162     rvec            *fshift;
 163     static  gmx_bool bFirst = TRUE;
 164
 165     lam   = 0;
 166     dvdl  = 0;
 167     tviol = 0;
 168     nviol = 0;
 169     mviol = 0;
 170     ndr   = 0;
 171     if (ntop)
 172     {
 173         reset5();
 174     }
 175     forceatoms = disres->iatoms;
 176     for (j = 0; (j < isize); j++)
 177     {
 178         vvindex[j] = 0;
 179     }
 180     nat = interaction_function[F_DISRES].nratoms+1;
 181     for (i = 0; (i < disres->nr); )
 182     {
 183         type  = forceatoms[i];
 184         n     = 0;
 185         label = forceparams[type].disres.label;
 186         if (debug)
 187         {
 188             fprintf(debug, "DISRE: ndr = %d, label = %d  i=%d, n =%d\n",
 189                     ndr, label, i, n);
 190         }
 191         if (ndr != label)
 192         {
 193             gmx_fatal(FARGS, "tpr inconsistency. ndr = %d, label = %d\n", ndr, label);
 194         }
 195         do
 196         {
 197             n += nat;
 198         }
 199         while (((i+n) < disres->nr) &&
 200                (forceparams[forceatoms[i+n]].disres.label == label));
 201
 202         calc_disres_R_6(n, &forceatoms[i], forceparams,
 203                         (const rvec*)x, pbc, fcd, NULL);
 204
 205         if (fcd->disres.Rt_6[0] <= 0)
 206         {
 207             gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[0]);
 208         }
 209
 210         rt = pow(fcd->disres.Rt_6[0], -1.0/6.0);
 211         dr[clust_id].aver1[ndr]  += rt;
 212         dr[clust_id].aver2[ndr]  += sqr(rt);
 213         drt = pow(rt, -3.0);
 214         dr[clust_id].aver_3[ndr] += drt;
 215         dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[0];
 216
 217         snew(fshift, SHIFTS);
 218         ener = interaction_function[F_DISRES].ifunc(n, &forceatoms[i], forceparams,
 219                                                     (const rvec*)x, f, fshift,
 220                                                     pbc, g, lam, &dvdl, NULL, fcd, NULL);
 221         sfree(fshift);
 222         viol = fcd->disres.sumviol;
 223
 224         if (viol > 0)
 225         {
 226             nviol++;
 227             if (ntop)
 228             {
 229                 add5(forceparams[type].disres.label, viol);
 230             }
 231             if (viol > mviol)
 232             {
 233                 mviol = viol;
 234             }
 235             tviol += viol;
 236             for (j = 0; (j < isize); j++)
 237             {
 238                 if (index[j] == forceparams[type].disres.label)
 239                 {
 240                     vvindex[j] = pow(fcd->disres.Rt_6[0], -1.0/6.0);
 241                 }
 242             }
 243         }
 244         ndr++;
 245         i   += n;
 246     }
 247     dr[clust_id].nv    = nviol;
 248     dr[clust_id].maxv  = mviol;
 249     dr[clust_id].sumv  = tviol;
 250     dr[clust_id].averv = tviol/ndr;
 251     dr[clust_id].nframes++;
 252
 253     if (bFirst)
 254     {
 255         fprintf(stderr, "\nThere are %d restraints and %d pairs\n",
 256                 ndr, disres->nr/nat);
 257         bFirst = FALSE;
 258     }
 259     if (ntop)
 260     {
 261         print5(log);
 262     }
 263 }
 264
 265 typedef struct {
 266     int      index;
 267     gmx_bool bCore;
 268     real     up1, r, rT3, rT6, viol, violT3, violT6;
 269 } t_dr_stats;
 270
 271 static int drs_comp(const void *a, const void *b)
 272 {
 273     t_dr_stats *da, *db;
 274
 275     da = (t_dr_stats *)a;
 276     db = (t_dr_stats *)b;
 277
 278     if (da->viol > db->viol)
 279     {
 280         return -1;
 281     }
 282     else if (da->viol < db->viol)
 283     {
 284         return 1;
 285     }
 286     else
 287     {
 288         return 0;
 289     }
 290 }
 291
 292 static void dump_dump(FILE *log, int ndr, t_dr_stats drs[])
 293 {
 294     static const char *core[] = { "All restraints", "Core restraints" };
 295     static const char *tp[]   = { "linear", "third power", "sixth power" };
 296     real               viol_tot, viol_max, viol = 0;
 297     gmx_bool           bCore;
 298     int                nviol, nrestr;
 299     int                i, kkk;
 300
 301     for (bCore = FALSE; (bCore <= TRUE); bCore++)
 302     {
 303         for (kkk = 0; (kkk < 3); kkk++)
 304         {
 305             viol_tot  = 0;
 306             viol_max  = 0;
 307             nviol     = 0;
 308             nrestr    = 0;
 309             for (i = 0; (i < ndr); i++)
 310             {
 311                 if (!bCore || (bCore && drs[i].bCore))
 312                 {
 313                     switch (kkk)
 314                     {
 315                         case 0:
 316                             viol = drs[i].viol;
 317                             break;
 318                         case 1:
 319                             viol = drs[i].violT3;
 320                             break;
 321                         case 2:
 322                             viol = drs[i].violT6;
 323                             break;
 324                         default:
 325                             gmx_incons("Dumping violations");
 326                     }
 327                     viol_max     = max(viol_max, viol);
 328                     if (viol > 0)
 329                     {
 330                         nviol++;
 331                     }
 332                     viol_tot  += viol;
 333                     nrestr++;
 334                 }
 335             }
 336             if ((nrestr > 0) || (bCore && (nrestr < ndr)))
 337             {
 338                 fprintf(log, "\n");
 339                 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
 340                 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
 341                 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
 342                 if (nrestr > 0)
 343                 {
 344                     fprintf(log, "Average violation: %8.3f nm\n", viol_tot/nrestr);
 345                 }
 346                 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
 347                 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
 348             }
 349         }
 350     }
 351 }
 352
 353 static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
 354 {
 355     int i;
 356
 357     fprintf(log, "Restr. Core     Up1     <r>   <rT3>   <rT6>  <viol><violT3><violT6>\n");
 358     for (i = 0; (i < ndr); i++)
 359     {
 360         if (bLinear  && (drs[i].viol == 0))
 361         {
 362             break;
 363         }
 364         fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
 365                 drs[i].index, yesno_names[drs[i].bCore],
 366                 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
 367                 drs[i].viol, drs[i].violT3, drs[i].violT6);
 368     }
 369 }
 370
 371 static gmx_bool is_core(int i, int isize, atom_id index[])
 372 {
 373     int      kk;
 374     gmx_bool bIC = FALSE;
 375
 376     for (kk = 0; !bIC && (kk < isize); kk++)
 377     {
 378         bIC = (index[kk] == i);
 379     }
 380
 381     return bIC;
 382 }
 383
 384 static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
 385                        t_iparams ip[], t_dr_result *dr,
 386                        int isize, atom_id index[], t_atoms *atoms)
 387 {
 388     int         i, j, nra;
 389     t_dr_stats *drs;
 390
 391     fprintf(log, "\n");
 392     fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
 393     snew(drs, ndr);
 394     j         = 0;
 395     nra       = interaction_function[F_DISRES].nratoms+1;
 396     for (i = 0; (i < ndr); i++)
 397     {
 398         /* Search for the appropriate restraint */
 399         for (; (j < disres->nr); j += nra)
 400         {
 401             if (ip[disres->iatoms[j]].disres.label == i)
 402             {
 403                 break;
 404             }
 405         }
 406         drs[i].index  = i;
 407         drs[i].bCore  = is_core(i, isize, index);
 408         drs[i].up1    = ip[disres->iatoms[j]].disres.up1;
 409         drs[i].r      = dr->aver1[i]/nsteps;
 410         drs[i].rT3    = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
 411         drs[i].rT6    = pow(dr->aver_6[i]/nsteps, -1.0/6.0);
 412         drs[i].viol   = max(0, drs[i].r-drs[i].up1);
 413         drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
 414         drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
 415         if (atoms)
 416         {
 417             int j1 = disres->iatoms[j+1];
 418             int j2 = disres->iatoms[j+2];
 419             atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
 420             atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
 421         }
 422     }
 423     dump_viol(log, ndr, drs, FALSE);
 424
 425     fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
 426     qsort(drs, ndr, sizeof(drs[0]), drs_comp);
 427     dump_viol(log, ndr, drs, TRUE);
 428
 429     dump_dump(log, ndr, drs);
 430
 431     sfree(drs);
 432 }
 433
 434 static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
 435                              t_iparams ip[], t_blocka *clust, t_dr_result dr[],
 436                              char *clust_name[], int isize, atom_id index[])
 437 {
 438     int         i, j, k, nra, mmm = 0;
 439     double      sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
 440     t_dr_stats *drs;
 441
 442     fprintf(fp, "\n");
 443     fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
 444     fprintf(fp, "Cluster  NFrames    SumV      MaxV     SumVT     MaxVT     SumVS     MaxVS\n");
 445
 446     snew(drs, ndr);
 447
 448     for (k = 0; (k < clust->nr); k++)
 449     {
 450         if (dr[k].nframes == 0)
 451         {
 452             continue;
 453         }
 454         if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
 455         {
 456             gmx_fatal(FARGS, "Inconsistency in cluster %s.\n"
 457                       "Found %d frames in trajectory rather than the expected %d\n",
 458                       clust_name[k], dr[k].nframes,
 459                       clust->index[k+1]-clust->index[k]);
 460         }
 461         if (!clust_name[k])
 462         {
 463             gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
 464         }
 465         j         = 0;
 466         nra       = interaction_function[F_DISRES].nratoms+1;
 467         sumV      = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
 468         for (i = 0; (i < ndr); i++)
 469         {
 470             /* Search for the appropriate restraint */
 471             for (; (j < disres->nr); j += nra)
 472             {
 473                 if (ip[disres->iatoms[j]].disres.label == i)
 474                 {
 475                     break;
 476                 }
 477             }
 478             drs[i].index  = i;
 479             drs[i].bCore  = is_core(i, isize, index);
 480             drs[i].up1    = ip[disres->iatoms[j]].disres.up1;
 481             drs[i].r      = dr[k].aver1[i]/dr[k].nframes;
 482             if ((dr[k].aver_3[i] <= 0) || (dr[k].aver_3[i] != dr[k].aver_3[i]))
 483             {
 484                 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
 485             }
 486             drs[i].rT3    = pow(dr[k].aver_3[i]/dr[k].nframes, -1.0/3.0);
 487             drs[i].rT6    = pow(dr[k].aver_6[i]/dr[k].nframes, -1.0/6.0);
 488             drs[i].viol   = max(0, drs[i].r-drs[i].up1);
 489             drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
 490             drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
 491             sumV         += drs[i].viol;
 492             sumVT3       += drs[i].violT3;
 493             sumVT6       += drs[i].violT6;
 494             maxV          = max(maxV, drs[i].viol);
 495             maxVT3        = max(maxVT3, drs[i].violT3);
 496             maxVT6        = max(maxVT6, drs[i].violT6);
 497         }
 498         if (strcmp(clust_name[k], "1000") == 0)
 499         {
 500             mmm++;
 501         }
 502         fprintf(fp, "%-10s%6d%8.3f  %8.3f  %8.3f  %8.3f  %8.3f  %8.3f\n",
 503                 clust_name[k],
 504                 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
 505
 506     }
 507     fflush(fp);
 508     sfree(drs);
 509 }
 510
 511 static void init_dr_res(t_dr_result *dr, int ndr)
 512 {
 513     snew(dr->aver1, ndr+1);
 514     snew(dr->aver2, ndr+1);
 515     snew(dr->aver_3, ndr+1);
 516     snew(dr->aver_6, ndr+1);
 517     dr->nv      = 0;
 518     dr->nframes = 0;
 519     dr->sumv    = 0;
 520     dr->maxv    = 0;
 521     dr->averv   = 0;
 522 }
 523
 524 static void dump_disre_matrix(const char *fn, t_dr_result *dr, int ndr,
 525                               int nsteps, t_idef *idef, gmx_mtop_t *mtop,
 526                               real max_dr, int nlevels, gmx_bool bThird)
 527 {
 528     FILE      *fp;
 529     int       *resnr;
 530     int        n_res, a_offset, mb, mol, a;
 531     t_atoms   *atoms;
 532     int        iii, i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
 533     atom_id    ai, aj, *ptr;
 534     real     **matrix, *t_res, hi, *w_dr, rav, rviol;
 535     t_rgb      rlo = { 1, 1, 1 };
 536     t_rgb      rhi = { 0, 0, 0 };
 537     if (fn == NULL)
 538     {
 539         return;
 540     }
 541
 542     snew(resnr, mtop->natoms);
 543     n_res    = 0;
 544     a_offset = 0;
 545     for (mb = 0; mb < mtop->nmolblock; mb++)
 546     {
 547         atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
 548         for (mol = 0; mol < mtop->molblock[mb].nmol; mol++)
 549         {
 550             for (a = 0; a < atoms->nr; a++)
 551             {
 552                 resnr[a_offset+a] = n_res + atoms->atom[a].resind;
 553             }
 554             n_res    += atoms->nres;
 555             a_offset += atoms->nr;
 556         }
 557     }
 558
 559     snew(t_res, n_res);
 560     for (i = 0; (i < n_res); i++)
 561     {
 562         t_res[i] = i+1;
 563     }
 564     snew(matrix, n_res);
 565     for (i = 0; (i < n_res); i++)
 566     {
 567         snew(matrix[i], n_res);
 568     }
 569     nratoms = interaction_function[F_DISRES].nratoms;
 570     nra     = (idef->il[F_DISRES].nr/(nratoms+1));
 571     snew(ptr, nra+1);
 572     index   = 0;
 573     nlabel  = 0;
 574     ptr[0]  = 0;
 575     snew(w_dr, ndr);
 576     for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms+1)
 577     {
 578         tp       = idef->il[F_DISRES].iatoms[i];
 579         label    = idef->iparams[tp].disres.label;
 580
 581         if (label != index)
 582         {
 583             /* Set index pointer */
 584             ptr[index+1] = i;
 585             if (nlabel <= 0)
 586             {
 587                 gmx_fatal(FARGS, "nlabel is %d, label = %d", nlabel, label);
 588             }
 589             if (index >= ndr)
 590             {
 591                 gmx_fatal(FARGS, "ndr = %d, index = %d", ndr, index);
 592             }
 593             /* Update the weight */
 594             w_dr[index] = 1.0/nlabel;
 595             index       = label;
 596             nlabel      = 1;
 597         }
 598         else
 599         {
 600             nlabel++;
 601         }
 602     }
 603     printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
 604     hi = 0;
 605     for (i = 0; (i < ndr); i++)
 606     {
 607         for (j = ptr[i]; (j < ptr[i+1]); j += nratoms+1)
 608         {
 609             tp  = idef->il[F_DISRES].iatoms[j];
 610             ai  = idef->il[F_DISRES].iatoms[j+1];
 611             aj  = idef->il[F_DISRES].iatoms[j+2];
 612
 613             ri = resnr[ai];
 614             rj = resnr[aj];
 615             if (bThird)
 616             {
 617                 rav = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
 618             }
 619             else
 620             {
 621                 rav = dr->aver1[i]/nsteps;
 622             }
 623             if (debug)
 624             {
 625                 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
 626             }
 627             rviol           = max(0, rav-idef->iparams[tp].disres.up1);
 628             matrix[ri][rj] += w_dr[i]*rviol;
 629             matrix[rj][ri] += w_dr[i]*rviol;
 630             hi              = max(hi, matrix[ri][rj]);
 631             hi              = max(hi, matrix[rj][ri]);
 632         }
 633     }
 634
 635     sfree(resnr);
 636
 637     if (max_dr > 0)
 638     {
 639         if (hi > max_dr)
 640         {
 641             printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n", max_dr, hi);
 642         }
 643         hi = max_dr;
 644     }
 645     printf("Highest level in the matrix will be %g\n", hi);
 646     fp = gmx_ffopen(fn, "w");
 647     write_xpm(fp, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue",
 648               n_res, n_res, t_res, t_res, matrix, 0, hi, rlo, rhi, &nlevels);
 649     gmx_ffclose(fp);
 650 }
 651
 652 int gmx_disre(int argc, char *argv[])
 653 {
 654     const char     *desc[] = {
 655         "[THISMODULE] computes violations of distance restraints.",
 656         "If necessary, all protons can be added to a protein molecule ",
 657         "using the [gmx-protonate] program.[PAR]",
 658         "The program always",
 659         "computes the instantaneous violations rather than time-averaged,",
 660         "because this analysis is done from a trajectory file afterwards",
 661         "it does not make sense to use time averaging. However,",
 662         "the time averaged values per restraint are given in the log file.[PAR]",
 663         "An index file may be used to select specific restraints for",
 664         "printing.[PAR]",
 665         "When the optional [TT]-q[tt] flag is given a [TT].pdb[tt] file coloured by the",
 666         "amount of average violations.[PAR]",
 667         "When the [TT]-c[tt] option is given, an index file will be read",
 668         "containing the frames in your trajectory corresponding to the clusters",
 669         "(defined in another manner) that you want to analyze. For these clusters",
 670         "the program will compute average violations using the third power",
 671         "averaging algorithm and print them in the log file."
 672     };
 673     static int      ntop      = 0;
 674     static int      nlevels   = 20;
 675     static real     max_dr    = 0;
 676     static gmx_bool bThird    = TRUE;
 677     t_pargs         pa[]      = {
 678         { "-ntop", FALSE, etINT,  {&ntop},
 679           "Number of large violations that are stored in the log file every step" },
 680         { "-maxdr", FALSE, etREAL, {&max_dr},
 681           "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
 682         { "-nlevels", FALSE, etINT, {&nlevels},
 683           "Number of levels in the matrix output" },
 684         { "-third", FALSE, etBOOL, {&bThird},
 685           "Use inverse third power averaging or linear for matrix output" }
 686     };
 687
 688     FILE           *out = NULL, *aver = NULL, *numv = NULL, *maxxv = NULL, *xvg = NULL;
 689     t_tpxheader     header;
 690     t_inputrec      ir;
 691     gmx_mtop_t      mtop;
 692     rvec           *xtop;
 693     gmx_localtop_t *top;
 694     t_atoms        *atoms = NULL;
 695     t_fcdata        fcd;
 696     t_nrnb          nrnb;
 697     t_graph        *g;
 698     int             ntopatoms, natoms, i, j, kkk;
 699     t_trxstatus    *status;
 700     real            t;
 701     rvec           *x, *f, *xav = NULL;
 702     matrix          box;
 703     gmx_bool        bPDB;
 704     int             isize;
 705     atom_id        *index = NULL, *ind_fit = NULL;
 706     char           *grpname;
 707     t_cluster_ndx  *clust = NULL;
 708     t_dr_result     dr, *dr_clust = NULL;
 709     char          **leg;
 710     real           *vvindex = NULL, *w_rls = NULL;
 711     t_mdatoms      *mdatoms;
 712     t_pbc           pbc, *pbc_null;
 713     int             my_clust;
 714     FILE           *fplog;
 715     output_env_t    oenv;
 716     gmx_rmpbc_t     gpbc = NULL;
 717
 718     t_filenm        fnm[] = {
 719         { efTPR, NULL, NULL, ffREAD },
 720         { efTRX, "-f", NULL, ffREAD },
 721         { efXVG, "-ds", "drsum",  ffWRITE },
 722         { efXVG, "-da", "draver", ffWRITE },
 723         { efXVG, "-dn", "drnum",  ffWRITE },
 724         { efXVG, "-dm", "drmax",  ffWRITE },
 725         { efXVG, "-dr", "restr",  ffWRITE },
 726         { efLOG, "-l",  "disres", ffWRITE },
 727         { efNDX, NULL,  "viol",   ffOPTRD },
 728         { efPDB, "-q",  "viol",   ffOPTWR },
 729         { efNDX, "-c",  "clust",  ffOPTRD },
 730         { efXPM, "-x",  "matrix", ffOPTWR }
 731     };
 732 #define NFILE asize(fnm)
 733
 734     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 735                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 736     {
 737         return 0;
 738     }
 739
 740     fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
 741
 742     if (ntop)
 743     {
 744         init5(ntop);
 745     }
 746
 747     read_tpxheader(ftp2fn(efTPR, NFILE, fnm), &header, FALSE, NULL, NULL);
 748     snew(xtop, header.natoms);
 749     read_tpx(ftp2fn(efTPR, NFILE, fnm), &ir, box, &ntopatoms, xtop, NULL, NULL, &mtop);
 750     bPDB = opt2bSet("-q", NFILE, fnm);
 751     if (bPDB)
 752     {
 753         snew(xav, ntopatoms);
 754         snew(ind_fit, ntopatoms);
 755         snew(w_rls, ntopatoms);
 756         for (kkk = 0; (kkk < ntopatoms); kkk++)
 757         {
 758             w_rls[kkk]   = 1;
 759             ind_fit[kkk] = kkk;
 760         }
 761
 762         snew(atoms, 1);
 763         *atoms = gmx_mtop_global_atoms(&mtop);
 764
 765         if (atoms->pdbinfo == NULL)
 766         {
 767             snew(atoms->pdbinfo, atoms->nr);
 768         }
 769     }
 770
 771     top = gmx_mtop_generate_local_top(&mtop, &ir);
 772
 773     g        = NULL;
 774     pbc_null = NULL;
 775     if (ir.ePBC != epbcNONE)
 776     {
 777         if (ir.bPeriodicMols)
 778         {
 779             pbc_null = &pbc;
 780         }
 781         else
 782         {
 783             g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
 784         }
 785     }
 786
 787     if (ftp2bSet(efNDX, NFILE, fnm))
 788     {
 789         rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 790         xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
 791                        "nm", oenv);
 792         snew(vvindex, isize);
 793         snew(leg, isize);
 794         for (i = 0; (i < isize); i++)
 795         {
 796             index[i]++;
 797             snew(leg[i], 12);
 798             sprintf(leg[i], "index %d", index[i]);
 799         }
 800         xvgr_legend(xvg, isize, (const char**)leg, oenv);
 801     }
 802     else
 803     {
 804         isize = 0;
 805     }
 806
 807     ir.dr_tau = 0.0;
 808     init_disres(fplog, &mtop, &ir, NULL, &fcd, NULL, FALSE);
 809
 810     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 811     snew(f, 5*natoms);
 812
 813     init_dr_res(&dr, fcd.disres.nres);
 814     if (opt2bSet("-c", NFILE, fnm))
 815     {
 816         clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
 817         snew(dr_clust, clust->clust->nr+1);
 818         for (i = 0; (i <= clust->clust->nr); i++)
 819         {
 820             init_dr_res(&dr_clust[i], fcd.disres.nres);
 821         }
 822     }
 823     else
 824     {
 825         out = xvgropen(opt2fn("-ds", NFILE, fnm),
 826                        "Sum of Violations", "Time (ps)", "nm", oenv);
 827         aver = xvgropen(opt2fn("-da", NFILE, fnm),
 828                         "Average Violation", "Time (ps)", "nm", oenv);
 829         numv = xvgropen(opt2fn("-dn", NFILE, fnm),
 830                         "# Violations", "Time (ps)", "#", oenv);
 831         maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
 832                          "Largest Violation", "Time (ps)", "nm", oenv);
 833     }
 834
 835     mdatoms = init_mdatoms(fplog, &mtop, ir.efep != efepNO);
 836     atoms2md(&mtop, &ir, 0, NULL, mtop.natoms, mdatoms);
 837     update_mdatoms(mdatoms, ir.fepvals->init_lambda);
 838     init_nrnb(&nrnb);
 839     if (ir.ePBC != epbcNONE)
 840     {
 841         gpbc = gmx_rmpbc_init(&top->idef, ir.ePBC, natoms);
 842     }
 843
 844     j = 0;
 845     do
 846     {
 847         if (ir.ePBC != epbcNONE)
 848         {
 849             if (ir.bPeriodicMols)
 850             {
 851                 set_pbc(&pbc, ir.ePBC, box);
 852             }
 853             else
 854             {
 855                 gmx_rmpbc(gpbc, natoms, box, x);
 856             }
 857         }
 858
 859         if (clust)
 860         {
 861             if (j > clust->maxframe)
 862             {
 863                 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
 864             }
 865             my_clust = clust->inv_clust[j];
 866             range_check(my_clust, 0, clust->clust->nr);
 867             check_viol(fplog, &(top->idef.il[F_DISRES]),
 868                        top->idef.iparams,
 869                        x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
 870         }
 871         else
 872         {
 873             check_viol(fplog, &(top->idef.il[F_DISRES]),
 874                        top->idef.iparams,
 875                        x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
 876         }
 877         if (bPDB)
 878         {
 879             reset_x(atoms->nr, ind_fit, atoms->nr, NULL, x, w_rls);
 880             do_fit(atoms->nr, w_rls, x, x);
 881             if (j == 0)
 882             {
 883                 /* Store the first frame of the trajectory as 'characteristic'
 884                  * for colouring with violations.
 885                  */
 886                 for (kkk = 0; (kkk < atoms->nr); kkk++)
 887                 {
 888                     copy_rvec(x[kkk], xav[kkk]);
 889                 }
 890             }
 891         }
 892         if (!clust)
 893         {
 894             if (isize > 0)
 895             {
 896                 fprintf(xvg, "%10g", t);
 897                 for (i = 0; (i < isize); i++)
 898                 {
 899                     fprintf(xvg, "  %10g", vvindex[i]);
 900                 }
 901                 fprintf(xvg, "\n");
 902             }
 903             fprintf(out,  "%10g  %10g\n", t, dr.sumv);
 904             fprintf(aver, "%10g  %10g\n", t, dr.averv);
 905             fprintf(maxxv, "%10g  %10g\n", t, dr.maxv);
 906             fprintf(numv, "%10g  %10d\n", t, dr.nv);
 907         }
 908         j++;
 909     }
 910     while (read_next_x(oenv, status, &t, x, box));
 911     close_trj(status);
 912     if (ir.ePBC != epbcNONE)
 913     {
 914         gmx_rmpbc_done(gpbc);
 915     }
 916
 917     if (clust)
 918     {
 919         dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
 920                          top->idef.iparams, clust->clust, dr_clust,
 921                          clust->grpname, isize, index);
 922     }
 923     else
 924     {
 925         dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
 926                    top->idef.iparams, &dr, isize, index,
 927                    bPDB ? atoms : NULL);
 928         if (bPDB)
 929         {
 930             write_sto_conf(opt2fn("-q", NFILE, fnm),
 931                            "Coloured by average violation in Angstrom",
 932                            atoms, xav, NULL, ir.ePBC, box);
 933         }
 934         dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
 935                           j, &top->idef, &mtop, max_dr, nlevels, bThird);
 936         gmx_ffclose(out);
 937         gmx_ffclose(aver);
 938         gmx_ffclose(numv);
 939         gmx_ffclose(maxxv);
 940         if (isize > 0)
 941         {
 942             gmx_ffclose(xvg);
 943             do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");
 944         }
 945         do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
 946         do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
 947         do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
 948         do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
 949     }
 950
 951     gmx_log_close(fplog);
 952
 953     return 0;
 954 }