src/gromacs/gmxana/gmx_disre.c

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