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