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