src/gromacs/gmxana/gmx_disre.c

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