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