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