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