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