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,2018,2019,2020, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40 #include <cstdlib>
  41 #include <cstring>
  42
  43 #include <algorithm>
  44 #include <unordered_map>
  45
  46 #include "gromacs/commandline/pargs.h"
  47 #include "gromacs/commandline/viewit.h"
  48 #include "gromacs/fileio/confio.h"
  49 #include "gromacs/fileio/matio.h"
  50 #include "gromacs/fileio/pdbio.h"
  51 #include "gromacs/fileio/tpxio.h"
  52 #include "gromacs/fileio/trxio.h"
  53 #include "gromacs/fileio/xvgr.h"
  54 #include "gromacs/gmxana/gmx_ana.h"
  55 #include "gromacs/gmxana/gstat.h"
  56 #include "gromacs/listed_forces/disre.h"
  57 #include "gromacs/math/do_fit.h"
  58 #include "gromacs/math/functions.h"
  59 #include "gromacs/math/vec.h"
  60 #include "gromacs/mdlib/force.h"
  61 #include "gromacs/mdlib/mdatoms.h"
  62 #include "gromacs/mdtypes/fcdata.h"
  63 #include "gromacs/mdtypes/inputrec.h"
  64 #include "gromacs/mdtypes/md_enums.h"
  65 #include "gromacs/pbcutil/ishift.h"
  66 #include "gromacs/pbcutil/mshift.h"
  67 #include "gromacs/pbcutil/pbc.h"
  68 #include "gromacs/pbcutil/rmpbc.h"
  69 #include "gromacs/topology/index.h"
  70 #include "gromacs/topology/mtop_util.h"
  71 #include "gromacs/topology/topology.h"
  72 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  73 #include "gromacs/utility/arraysize.h"
  74 #include "gromacs/utility/fatalerror.h"
  75 #include "gromacs/utility/futil.h"
  76 #include "gromacs/utility/smalloc.h"
  77
  78 typedef struct
  79 {
  80     int  n;
  81     real v;
  82 } t_toppop;
  83
  84 static t_toppop* top  = nullptr;
  85 static int       ntop = 0;
  86
  87 typedef struct
  88 {
  89     int   nv, nframes;
  90     real  sumv, averv, maxv;
  91     real *aver1, *aver2, *aver_3, *aver_6;
  92 } t_dr_result;
  93
  94 static void init5(int n)
  95 {
  96     ntop = n;
  97     snew(top, ntop);
  98 }
  99
 100 static void reset5()
 101 {
 102     int i;
 103
 104     for (i = 0; (i < ntop); i++)
 105     {
 106         top[i].n = -1;
 107         top[i].v = 0;
 108     }
 109 }
 110
 111 static void add5(int ndr, real viol)
 112 {
 113     int i, mini;
 114
 115     mini = 0;
 116     for (i = 1; (i < ntop); i++)
 117     {
 118         if (top[i].v < top[mini].v)
 119         {
 120             mini = i;
 121         }
 122     }
 123     if (viol > top[mini].v)
 124     {
 125         top[mini].v = viol;
 126         top[mini].n = ndr;
 127     }
 128 }
 129
 130 static void print5(FILE* fp)
 131 {
 132     int i;
 133
 134     std::sort(top, top + ntop, [](const t_toppop& a, const t_toppop& b) { return a.v > b.v; }); // reverse sort
 135     fprintf(fp, "Index:");
 136     for (i = 0; (i < ntop); i++)
 137     {
 138         fprintf(fp, " %6d", top[i].n);
 139     }
 140     fprintf(fp, "\nViol: ");
 141     for (i = 0; (i < ntop); i++)
 142     {
 143         fprintf(fp, " %6.3f", top[i].v);
 144     }
 145     fprintf(fp, "\n");
 146 }
 147
 148 static void check_viol(FILE*       log,
 149                        t_ilist*    disres,
 150                        t_iparams   forceparams[],
 151                        rvec        x[],
 152                        rvec4       f[],
 153                        t_pbc*      pbc,
 154                        t_graph*    g,
 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     t_iatom*        forceatoms;
 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     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->nr); 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->nr);)
 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->nr)
 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, x, f, fshift, pbc, g, lam, &dvdl, nullptr, 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->nr / 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 t_ilist*      disres,
 376                        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->nr; 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 t_ilist*      disres,
 426                              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                               t_idef*           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].nr / (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].nr); 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     gmx_localtop_t    top;
 702     t_fcdata          fcd;
 703     t_graph*          g;
 704     int               i, j, kkk;
 705     t_trxstatus*      status;
 706     real              t;
 707     rvec *            x, *xav = nullptr;
 708     rvec4*            f;
 709     matrix            box;
 710     gmx_bool          bPDB;
 711     int               isize;
 712     int *             index = nullptr, *ind_fit = nullptr;
 713     char*             grpname;
 714     t_cluster_ndx*    clust = nullptr;
 715     t_dr_result       dr, *dr_clust = nullptr;
 716     char**            leg;
 717     real *            vvindex = nullptr, *w_rls = nullptr;
 718     t_pbc             pbc, *pbc_null;
 719     int               my_clust;
 720     FILE*             fplog;
 721     gmx_output_env_t* oenv;
 722     gmx_rmpbc_t       gpbc = nullptr;
 723
 724     t_filenm fnm[] = { { efTPR, nullptr, nullptr, ffREAD }, { efTRX, "-f", nullptr, ffREAD },
 725                        { efXVG, "-ds", "drsum", ffWRITE },  { efXVG, "-da", "draver", ffWRITE },
 726                        { efXVG, "-dn", "drnum", ffWRITE },  { efXVG, "-dm", "drmax", ffWRITE },
 727                        { efXVG, "-dr", "restr", ffWRITE },  { efLOG, "-l", "disres", ffWRITE },
 728                        { efNDX, nullptr, "viol", ffOPTRD }, { efPDB, "-q", "viol", ffOPTWR },
 729                        { efNDX, "-c", "clust", ffOPTRD },   { efXPM, "-x", "matrix", ffOPTWR } };
 730 #define NFILE asize(fnm)
 731
 732     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa,
 733                            asize(desc), desc, 0, nullptr, &oenv))
 734     {
 735         return 0;
 736     }
 737
 738     fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
 739
 740     if (ntop)
 741     {
 742         init5(ntop);
 743     }
 744
 745     t_inputrec  irInstance;
 746     t_inputrec* ir = &irInstance;
 747
 748     gmx::TopologyInformation topInfo;
 749     topInfo.fillFromInputFile(ftp2fn(efTPR, NFILE, fnm));
 750     int          ntopatoms = topInfo.mtop()->natoms;
 751     AtomsDataPtr atoms;
 752     bPDB = opt2bSet("-q", NFILE, fnm);
 753     if (bPDB)
 754     {
 755         snew(xav, ntopatoms);
 756         snew(ind_fit, ntopatoms);
 757         snew(w_rls, ntopatoms);
 758         for (kkk = 0; (kkk < ntopatoms); kkk++)
 759         {
 760             w_rls[kkk]   = 1;
 761             ind_fit[kkk] = kkk;
 762         }
 763
 764         atoms = topInfo.copyAtoms();
 765
 766         if (atoms->pdbinfo == nullptr)
 767         {
 768             snew(atoms->pdbinfo, atoms->nr);
 769         }
 770         atoms->havePdbInfo = TRUE;
 771     }
 772
 773     gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO);
 774
 775     g        = nullptr;
 776     pbc_null = nullptr;
 777     if (ir->ePBC != epbcNONE)
 778     {
 779         if (ir->bPeriodicMols)
 780         {
 781             pbc_null = &pbc;
 782         }
 783         else
 784         {
 785             g = mk_graph(fplog, &top.idef, 0, ntopatoms, FALSE, FALSE);
 786         }
 787     }
 788
 789     if (ftp2bSet(efNDX, NFILE, fnm))
 790     {
 791         /* TODO: Nothing is written to this file if -c is provided, but it is
 792          * still opened... */
 793         rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 794         xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)", "nm", oenv);
 795         snew(vvindex, isize);
 796         snew(leg, isize);
 797         for (i = 0; (i < isize); i++)
 798         {
 799             index[i]++;
 800             snew(leg[i], 12);
 801             sprintf(leg[i], "index %d", index[i]);
 802         }
 803         xvgr_legend(xvg, isize, leg, oenv);
 804     }
 805     else
 806     {
 807         isize = 0;
 808     }
 809
 810     ir->dr_tau = 0.0;
 811     init_disres(fplog, topInfo.mtop(), ir, nullptr, nullptr, &fcd, nullptr, FALSE);
 812
 813     int natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 814     snew(f, 5 * natoms);
 815
 816     init_dr_res(&dr, fcd.disres.nres);
 817     if (opt2bSet("-c", NFILE, fnm))
 818     {
 819         clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
 820         snew(dr_clust, clust->clust->nr + 1);
 821         for (i = 0; (i <= clust->clust->nr); i++)
 822         {
 823             init_dr_res(&dr_clust[i], fcd.disres.nres);
 824         }
 825     }
 826     else
 827     {
 828         out   = xvgropen(opt2fn("-ds", NFILE, fnm), "Sum of Violations", "Time (ps)", "nm", oenv);
 829         aver  = xvgropen(opt2fn("-da", NFILE, fnm), "Average Violation", "Time (ps)", "nm", oenv);
 830         numv  = xvgropen(opt2fn("-dn", NFILE, fnm), "# Violations", "Time (ps)", "#", oenv);
 831         maxxv = xvgropen(opt2fn("-dm", NFILE, fnm), "Largest Violation", "Time (ps)", "nm", oenv);
 832     }
 833
 834     auto mdAtoms = gmx::makeMDAtoms(fplog, *topInfo.mtop(), *ir, false);
 835     atoms2md(topInfo.mtop(), ir, -1, nullptr, ntopatoms, mdAtoms.get());
 836     update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda);
 837     if (ir->ePBC != epbcNONE)
 838     {
 839         gpbc = gmx_rmpbc_init(&top.idef, ir->ePBC, natoms);
 840     }
 841
 842     j = 0;
 843     do
 844     {
 845         if (ir->ePBC != epbcNONE)
 846         {
 847             if (ir->bPeriodicMols)
 848             {
 849                 set_pbc(&pbc, ir->ePBC, box);
 850             }
 851             else
 852             {
 853                 gmx_rmpbc(gpbc, natoms, box, x);
 854             }
 855         }
 856
 857         if (clust)
 858         {
 859             if (j > clust->maxframe)
 860             {
 861                 gmx_fatal(FARGS,
 862                           "There are more frames in the trajectory than in the cluster index file. "
 863                           "t = %8f\n",
 864                           t);
 865             }
 866             my_clust = clust->inv_clust[j];
 867             range_check(my_clust, 0, clust->clust->nr);
 868             check_viol(fplog, &(top.idef.il[F_DISRES]), top.idef.iparams, x, f, pbc_null, g,
 869                        dr_clust, my_clust, isize, index, vvindex, &fcd);
 870         }
 871         else
 872         {
 873             check_viol(fplog, &(top.idef.il[F_DISRES]), top.idef.iparams, x, f, pbc_null, g, &dr, 0,
 874                        isize, index, vvindex, &fcd);
 875         }
 876         if (bPDB)
 877         {
 878             reset_x(atoms->nr, ind_fit, atoms->nr, nullptr, x, w_rls);
 879             do_fit(atoms->nr, w_rls, x, x);
 880             if (j == 0)
 881             {
 882                 /* Store the first frame of the trajectory as 'characteristic'
 883                  * for colouring with violations.
 884                  */
 885                 for (kkk = 0; (kkk < atoms->nr); kkk++)
 886                 {
 887                     copy_rvec(x[kkk], xav[kkk]);
 888                 }
 889             }
 890         }
 891         if (!clust)
 892         {
 893             if (isize > 0)
 894             {
 895                 fprintf(xvg, "%10g", t);
 896                 for (i = 0; (i < isize); i++)
 897                 {
 898                     fprintf(xvg, "  %10g", vvindex[i]);
 899                 }
 900                 fprintf(xvg, "\n");
 901             }
 902             fprintf(out, "%10g  %10g\n", t, dr.sumv);
 903             fprintf(aver, "%10g  %10g\n", t, dr.averv);
 904             fprintf(maxxv, "%10g  %10g\n", t, dr.maxv);
 905             fprintf(numv, "%10g  %10d\n", t, dr.nv);
 906         }
 907         j++;
 908     } while (read_next_x(oenv, status, &t, x, box));
 909     close_trx(status);
 910     if (ir->ePBC != epbcNONE)
 911     {
 912         gmx_rmpbc_done(gpbc);
 913     }
 914
 915     if (clust)
 916     {
 917         dump_clust_stats(fplog, fcd.disres, &(top.idef.il[F_DISRES]), top.idef.iparams,
 918                          clust->clust, dr_clust, clust->grpname, isize, index);
 919     }
 920     else
 921     {
 922         dump_stats(fplog, j, fcd.disres, &(top.idef.il[F_DISRES]), top.idef.iparams, &dr, isize,
 923                    index, bPDB ? atoms.get() : nullptr);
 924         if (bPDB)
 925         {
 926             write_sto_conf(opt2fn("-q", NFILE, fnm), "Coloured by average violation in Angstrom",
 927                            atoms.get(), xav, nullptr, ir->ePBC, box);
 928         }
 929         dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres, j, &top.idef,
 930                           topInfo.mtop(), max_dr, nlevels, bThird);
 931         xvgrclose(out);
 932         xvgrclose(aver);
 933         xvgrclose(numv);
 934         xvgrclose(maxxv);
 935         do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
 936         do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
 937         do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
 938         do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
 939     }
 940     if (isize > 0)
 941     {
 942         xvgrclose(xvg);
 943         if (!clust)
 944         {
 945             do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");
 946         }
 947     }
 948
 949     gmx_ffclose(fplog);
 950
 951     return 0;
 952 }