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