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