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