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