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