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, 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/mdlib/mdrun.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/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(void)
  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         interaction_function[F_DISRES].ifunc(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     t_tpxheader       header;
 662     gmx_mtop_t        mtop;
 663     rvec             *xtop;
 664     gmx_localtop_t   *top;
 665     t_atoms          *atoms = nullptr;
 666     t_fcdata          fcd;
 667     t_nrnb            nrnb;
 668     t_graph          *g;
 669     int               ntopatoms, natoms, i, j, kkk;
 670     t_trxstatus      *status;
 671     real              t;
 672     rvec             *x, *xav = nullptr;
 673     rvec4            *f;
 674     matrix            box;
 675     gmx_bool          bPDB;
 676     int               isize;
 677     int              *index = nullptr, *ind_fit = nullptr;
 678     char             *grpname;
 679     t_cluster_ndx    *clust = nullptr;
 680     t_dr_result       dr, *dr_clust = nullptr;
 681     char            **leg;
 682     real             *vvindex = nullptr, *w_rls = nullptr;
 683     t_pbc             pbc, *pbc_null;
 684     int               my_clust;
 685     FILE             *fplog;
 686     gmx_output_env_t *oenv;
 687     gmx_rmpbc_t       gpbc = nullptr;
 688
 689     t_filenm          fnm[] = {
 690         { efTPR, nullptr, nullptr, ffREAD },
 691         { efTRX, "-f", nullptr, ffREAD },
 692         { efXVG, "-ds", "drsum",  ffWRITE },
 693         { efXVG, "-da", "draver", ffWRITE },
 694         { efXVG, "-dn", "drnum",  ffWRITE },
 695         { efXVG, "-dm", "drmax",  ffWRITE },
 696         { efXVG, "-dr", "restr",  ffWRITE },
 697         { efLOG, "-l",  "disres", ffWRITE },
 698         { efNDX, nullptr,  "viol",   ffOPTRD },
 699         { efPDB, "-q",  "viol",   ffOPTWR },
 700         { efNDX, "-c",  "clust",  ffOPTRD },
 701         { efXPM, "-x",  "matrix", ffOPTWR }
 702     };
 703 #define NFILE asize(fnm)
 704
 705     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 706                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 707     {
 708         return 0;
 709     }
 710
 711     fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
 712
 713     if (ntop)
 714     {
 715         init5(ntop);
 716     }
 717
 718     t_inputrec      irInstance;
 719     t_inputrec     *ir = &irInstance;
 720
 721     read_tpxheader(ftp2fn(efTPR, NFILE, fnm), &header, FALSE);
 722     snew(xtop, header.natoms);
 723     read_tpx(ftp2fn(efTPR, NFILE, fnm), ir, box, &ntopatoms, xtop, nullptr, &mtop);
 724     bPDB = opt2bSet("-q", NFILE, fnm);
 725     if (bPDB)
 726     {
 727         snew(xav, ntopatoms);
 728         snew(ind_fit, ntopatoms);
 729         snew(w_rls, ntopatoms);
 730         for (kkk = 0; (kkk < ntopatoms); kkk++)
 731         {
 732             w_rls[kkk]   = 1;
 733             ind_fit[kkk] = kkk;
 734         }
 735
 736         snew(atoms, 1);
 737         *atoms = gmx_mtop_global_atoms(&mtop);
 738
 739         if (atoms->pdbinfo == nullptr)
 740         {
 741             snew(atoms->pdbinfo, atoms->nr);
 742         }
 743         atoms->havePdbInfo = TRUE;
 744     }
 745
 746     top = gmx_mtop_generate_local_top(&mtop, ir->efep != efepNO);
 747
 748     g        = nullptr;
 749     pbc_null = nullptr;
 750     if (ir->ePBC != epbcNONE)
 751     {
 752         if (ir->bPeriodicMols)
 753         {
 754             pbc_null = &pbc;
 755         }
 756         else
 757         {
 758             g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
 759         }
 760     }
 761
 762     if (ftp2bSet(efNDX, NFILE, fnm))
 763     {
 764         /* TODO: Nothing is written to this file if -c is provided, but it is
 765          * still opened... */
 766         rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 767         xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
 768                        "nm", oenv);
 769         snew(vvindex, isize);
 770         snew(leg, isize);
 771         for (i = 0; (i < isize); i++)
 772         {
 773             index[i]++;
 774             snew(leg[i], 12);
 775             sprintf(leg[i], "index %d", index[i]);
 776         }
 777         xvgr_legend(xvg, isize, leg, oenv);
 778     }
 779     else
 780     {
 781         isize = 0;
 782     }
 783
 784     ir->dr_tau = 0.0;
 785     init_disres(fplog, &mtop, ir, nullptr, nullptr, &fcd, nullptr, FALSE);
 786
 787     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 788     snew(f, 5*natoms);
 789
 790     init_dr_res(&dr, fcd.disres.nres);
 791     if (opt2bSet("-c", NFILE, fnm))
 792     {
 793         clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
 794         snew(dr_clust, clust->clust->nr+1);
 795         for (i = 0; (i <= clust->clust->nr); i++)
 796         {
 797             init_dr_res(&dr_clust[i], fcd.disres.nres);
 798         }
 799     }
 800     else
 801     {
 802         out = xvgropen(opt2fn("-ds", NFILE, fnm),
 803                        "Sum of Violations", "Time (ps)", "nm", oenv);
 804         aver = xvgropen(opt2fn("-da", NFILE, fnm),
 805                         "Average Violation", "Time (ps)", "nm", oenv);
 806         numv = xvgropen(opt2fn("-dn", NFILE, fnm),
 807                         "# Violations", "Time (ps)", "#", oenv);
 808         maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
 809                          "Largest Violation", "Time (ps)", "nm", oenv);
 810     }
 811
 812     auto mdAtoms = gmx::makeMDAtoms(fplog, mtop, *ir, false);
 813     atoms2md(&mtop, ir, -1, nullptr, mtop.natoms, mdAtoms.get());
 814     update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda);
 815     init_nrnb(&nrnb);
 816     if (ir->ePBC != epbcNONE)
 817     {
 818         gpbc = gmx_rmpbc_init(&top->idef, ir->ePBC, natoms);
 819     }
 820
 821     j = 0;
 822     do
 823     {
 824         if (ir->ePBC != epbcNONE)
 825         {
 826             if (ir->bPeriodicMols)
 827             {
 828                 set_pbc(&pbc, ir->ePBC, box);
 829             }
 830             else
 831             {
 832                 gmx_rmpbc(gpbc, natoms, box, x);
 833             }
 834         }
 835
 836         if (clust)
 837         {
 838             if (j > clust->maxframe)
 839             {
 840                 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
 841             }
 842             my_clust = clust->inv_clust[j];
 843             range_check(my_clust, 0, clust->clust->nr);
 844             check_viol(fplog, &(top->idef.il[F_DISRES]),
 845                        top->idef.iparams,
 846                        x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
 847         }
 848         else
 849         {
 850             check_viol(fplog, &(top->idef.il[F_DISRES]),
 851                        top->idef.iparams,
 852                        x, f, pbc_null, g, &dr, 0, 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     }
 887     while (read_next_x(oenv, status, &t, x, box));
 888     close_trx(status);
 889     if (ir->ePBC != epbcNONE)
 890     {
 891         gmx_rmpbc_done(gpbc);
 892     }
 893
 894     if (clust)
 895     {
 896         dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
 897                          top->idef.iparams, clust->clust, dr_clust,
 898                          clust->grpname, isize, index);
 899     }
 900     else
 901     {
 902         dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
 903                    top->idef.iparams, &dr, isize, index,
 904                    bPDB ? atoms : nullptr);
 905         if (bPDB)
 906         {
 907             write_sto_conf(opt2fn("-q", NFILE, fnm),
 908                            "Coloured by average violation in Angstrom",
 909                            atoms, xav, nullptr, ir->ePBC, box);
 910         }
 911         dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
 912                           j, &top->idef, &mtop, max_dr, nlevels, bThird);
 913         xvgrclose(out);
 914         xvgrclose(aver);
 915         xvgrclose(numv);
 916         xvgrclose(maxxv);
 917         do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
 918         do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
 919         do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
 920         do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
 921     }
 922     if (isize > 0)
 923     {
 924         xvgrclose(xvg);
 925         if (!clust)
 926         {
 927             do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");
 928         }
 929     }
 930
 931     gmx_ffclose(fplog);
 932
 933     return 0;
 934 }