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