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