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