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