src/gromacs/gmxana/gmx_disre.cpp

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