src/gromacs/gmxana/gmx_nmr.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstdlib>
  42 #include <cstring>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/commandline/pargs.h"
  47 #include "gromacs/commandline/viewit.h"
  48 #include "gromacs/correlationfunctions/autocorr.h"
  49 #include "gromacs/fileio/enxio.h"
  50 #include "gromacs/fileio/gmxfio.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/math/functions.h"
  57 #include "gromacs/math/units.h"
  58 #include "gromacs/math/vec.h"
  59 #include "gromacs/mdlib/energyoutput.h"
  60 #include "gromacs/mdtypes/inputrec.h"
  61 #include "gromacs/mdtypes/md_enums.h"
  62 #include "gromacs/topology/ifunc.h"
  63 #include "gromacs/topology/mtop_lookup.h"
  64 #include "gromacs/topology/mtop_util.h"
  65 #include "gromacs/topology/topology.h"
  66 #include "gromacs/trajectory/energyframe.h"
  67 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  68 #include "gromacs/utility/arraysize.h"
  69 #include "gromacs/utility/cstringutil.h"
  70 #include "gromacs/utility/fatalerror.h"
  71 #include "gromacs/utility/gmxassert.h"
  72 #include "gromacs/utility/pleasecite.h"
  73 #include "gromacs/utility/smalloc.h"
  74 #include "gromacs/utility/strconvert.h"
  75
  76 static real minthird = -1.0 / 3.0, minsixth = -1.0 / 6.0;
  77
  78 static double mypow(double x, double y)
  79 {
  80     if (x > 0)
  81     {
  82         return std::pow(x, y);
  83     }
  84     else
  85     {
  86         return 0.0;
  87     }
  88 }
  89
  90 static real blk_value(t_enxblock* blk, int sub, int index)
  91 {
  92     range_check(index, 0, blk->sub[sub].nr);
  93     if (blk->sub[sub].type == xdr_datatype_float)
  94     {
  95         return blk->sub[sub].fval[index];
  96     }
  97     else if (blk->sub[sub].type == xdr_datatype_double)
  98     {
  99         return blk->sub[sub].dval[index];
 100     }
 101     else
 102     {
 103         gmx_incons("Unknown datatype in t_enxblock");
 104     }
 105 }
 106
 107 static int* select_it(int nre, char* nm[], int* nset)
 108 {
 109     gmx_bool* bE;
 110     int       n, k, j, i;
 111     int*      set;
 112     gmx_bool  bVerbose = TRUE;
 113
 114     if ((getenv("GMX_ENER_VERBOSE")) != nullptr)
 115     {
 116         bVerbose = FALSE;
 117     }
 118
 119     fprintf(stderr, "Select the terms you want from the following list\n");
 120     fprintf(stderr, "End your selection with 0\n");
 121
 122     if (bVerbose)
 123     {
 124         for (k = 0; (k < nre);)
 125         {
 126             for (j = 0; (j < 4) && (k < nre); j++, k++)
 127             {
 128                 fprintf(stderr, " %3d=%14s", k + 1, nm[k]);
 129             }
 130             fprintf(stderr, "\n");
 131         }
 132     }
 133
 134     snew(bE, nre);
 135     do
 136     {
 137         if (1 != scanf("%d", &n))
 138         {
 139             gmx_fatal(FARGS, "Error reading user input");
 140         }
 141         if ((n > 0) && (n <= nre))
 142         {
 143             bE[n - 1] = TRUE;
 144         }
 145     } while (n != 0);
 146
 147     snew(set, nre);
 148     for (i = (*nset) = 0; (i < nre); i++)
 149     {
 150         if (bE[i])
 151         {
 152             set[(*nset)++] = i;
 153         }
 154     }
 155
 156     sfree(bE);
 157
 158     return set;
 159 }
 160
 161 static void get_orires_parms(const char* topnm, t_inputrec* ir, int* nor, int* nex, int** label, real** obs)
 162 {
 163     gmx_mtop_t mtop;
 164     t_topology top;
 165     t_iparams* ip;
 166     int        natoms, i;
 167     t_iatom*   iatom;
 168     int        nb;
 169     matrix     box;
 170
 171     read_tpx(topnm, ir, box, &natoms, nullptr, nullptr, &mtop);
 172     top = gmx_mtop_t_to_t_topology(&mtop, FALSE);
 173
 174     ip    = top.idef.iparams;
 175     iatom = top.idef.il[F_ORIRES].iatoms;
 176
 177     /* Count how many distance restraint there are... */
 178     nb = top.idef.il[F_ORIRES].nr;
 179     if (nb == 0)
 180     {
 181         gmx_fatal(FARGS, "No orientation restraints in topology!\n");
 182     }
 183
 184     *nor = nb / 3;
 185     *nex = 0;
 186     snew(*label, *nor);
 187     snew(*obs, *nor);
 188     for (i = 0; i < nb; i += 3)
 189     {
 190         (*label)[i / 3] = ip[iatom[i]].orires.label;
 191         (*obs)[i / 3]   = ip[iatom[i]].orires.obs;
 192         if (ip[iatom[i]].orires.ex >= *nex)
 193         {
 194             *nex = ip[iatom[i]].orires.ex + 1;
 195         }
 196     }
 197     fprintf(stderr, "Found %d orientation restraints with %d experiments", *nor, *nex);
 198     done_top_mtop(&top, &mtop);
 199 }
 200
 201 static int get_bounds(real** bounds, int** index, int** dr_pair, int* npairs, const InteractionDefinitions& idef)
 202 {
 203     int   i, j, k, type, ftype, natom;
 204     real* b;
 205     int * ind, *pair;
 206     int   nb, label1;
 207
 208     gmx::ArrayRef<const t_functype> functype = idef.functype;
 209     gmx::ArrayRef<const t_iparams>  iparams  = idef.iparams;
 210
 211     /* Count how many distance restraint there are... */
 212     nb = idef.il[F_DISRES].size();
 213     if (nb == 0)
 214     {
 215         gmx_fatal(FARGS, "No distance restraints in topology!\n");
 216     }
 217
 218     /* Allocate memory */
 219     snew(b, nb);
 220     snew(ind, nb);
 221     snew(pair, nb + 1);
 222
 223     /* Fill the bound array */
 224     nb = 0;
 225     for (gmx::index i = 0; i < functype.ssize(); i++)
 226     {
 227         ftype = functype[i];
 228         if (ftype == F_DISRES)
 229         {
 230
 231             label1  = iparams[i].disres.label;
 232             b[nb]   = iparams[i].disres.up1;
 233             ind[nb] = label1;
 234             nb++;
 235         }
 236     }
 237     *bounds = b;
 238
 239     /* Fill the index array */
 240     label1                          = -1;
 241     const InteractionList&   disres = idef.il[F_DISRES];
 242     gmx::ArrayRef<const int> iatom  = disres.iatoms;
 243     for (i = j = k = 0; (i < disres.size());)
 244     {
 245         type  = iatom[i];
 246         ftype = functype[type];
 247         natom = interaction_function[ftype].nratoms + 1;
 248         if (label1 != iparams[type].disres.label)
 249         {
 250             pair[j] = k;
 251             label1  = iparams[type].disres.label;
 252             j++;
 253         }
 254         k++;
 255         i += natom;
 256     }
 257     pair[j] = k;
 258     *npairs = k;
 259     if (j != nb)
 260     {
 261         gmx_incons("get_bounds for distance restraints");
 262     }
 263
 264     *index   = ind;
 265     *dr_pair = pair;
 266
 267     return nb;
 268 }
 269
 270 static void
 271 calc_violations(real rt[], real rav3[], int nb, const int index[], real bounds[], real* viol, double* st, double* sa)
 272 {
 273     const real sixth = 1.0 / 6.0;
 274     int        i, j;
 275     double     rsum, rav, sumaver, sumt;
 276
 277     sumaver = 0;
 278     sumt    = 0;
 279     for (i = 0; (i < nb); i++)
 280     {
 281         rsum = 0.0;
 282         rav  = 0.0;
 283         for (j = index[i]; (j < index[i + 1]); j++)
 284         {
 285             if (viol)
 286             {
 287                 viol[j] += mypow(rt[j], -3.0);
 288             }
 289             rav += gmx::square(rav3[j]);
 290             rsum += mypow(rt[j], -6);
 291         }
 292         rsum = std::max(0.0, mypow(rsum, -sixth) - bounds[i]);
 293         rav  = std::max(0.0, mypow(rav, -sixth) - bounds[i]);
 294
 295         sumt += rsum;
 296         sumaver += rav;
 297     }
 298     *st = sumt;
 299     *sa = sumaver;
 300 }
 301
 302 static void analyse_disre(const char*             voutfn,
 303                           int                     nframes,
 304                           real                    violaver[],
 305                           real                    bounds[],
 306                           int                     index[],
 307                           int                     pair[],
 308                           int                     nbounds,
 309                           const gmx_output_env_t* oenv)
 310 {
 311     FILE*  vout;
 312     double sum, sumt, sumaver;
 313     int    i, j;
 314
 315     /* Subtract bounds from distances, to calculate violations */
 316     calc_violations(violaver, violaver, nbounds, pair, bounds, nullptr, &sumt, &sumaver);
 317
 318 #ifdef DEBUG
 319     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n", sumaver);
 320     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sumt);
 321 #endif
 322     vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm", oenv);
 323     sum  = 0.0;
 324     sumt = 0.0;
 325     for (i = 0; (i < nbounds); i++)
 326     {
 327         /* Do ensemble averaging */
 328         sumaver = 0;
 329         for (j = pair[i]; (j < pair[i + 1]); j++)
 330         {
 331             sumaver += gmx::square(violaver[j] / real(nframes));
 332         }
 333         sumaver = std::max(0.0, mypow(sumaver, minsixth) - bounds[i]);
 334
 335         sumt += sumaver;
 336         sum = std::max(sum, sumaver);
 337         fprintf(vout, "%10d  %10.5e\n", index[i], sumaver);
 338     }
 339 #ifdef DEBUG
 340     for (j = 0; (j < dr.ndr); j++)
 341     {
 342         fprintf(vout, "%10d  %10.5e\n", j, mypow(violaver[j] / real(nframes), minthird));
 343     }
 344 #endif
 345     xvgrclose(vout);
 346
 347     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n", sumt);
 348     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
 349
 350     do_view(oenv, voutfn, "-graphtype bar");
 351 }
 352
 353 static void print_time(FILE* fp, double t)
 354 {
 355     fprintf(fp, "%12.6f", t);
 356 }
 357
 358 int gmx_nmr(int argc, char* argv[])
 359 {
 360     const char* desc[] = {
 361         "[THISMODULE] extracts distance or orientation restraint",
 362         "data from an energy file. The user is prompted to interactively",
 363         "select the desired terms.[PAR]",
 364
 365         "When the [TT]-viol[tt] option is set, the time averaged",
 366         "violations are plotted and the running time-averaged and",
 367         "instantaneous sum of violations are recalculated. Additionally",
 368         "running time-averaged and instantaneous distances between",
 369         "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
 370
 371         "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
 372         "[TT]-odt[tt] are used for analyzing orientation restraint data.",
 373         "The first two options plot the orientation, the last three the",
 374         "deviations of the orientations from the experimental values.",
 375         "The options that end on an 'a' plot the average over time",
 376         "as a function of restraint. The options that end on a 't'",
 377         "prompt the user for restraint label numbers and plot the data",
 378         "as a function of time. Option [TT]-odr[tt] plots the RMS",
 379         "deviation as a function of restraint.",
 380         "When the run used time or ensemble averaged orientation restraints,",
 381         "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
 382         "not ensemble-averaged orientations and deviations instead of",
 383         "the time and ensemble averages.[PAR]",
 384
 385         "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
 386         "tensor for each orientation restraint experiment. With option",
 387         "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
 388
 389
 390     };
 391     static gmx_bool bPrAll = FALSE;
 392     static gmx_bool bDp = FALSE, bOrinst = FALSE, bOvec = FALSE;
 393     static int      skip = 0;
 394     t_pargs         pa[] = {
 395         { "-dp", FALSE, etBOOL, { &bDp }, "Print energies in high precision" },
 396         { "-skip", FALSE, etINT, { &skip }, "Skip number of frames between data points" },
 397         { "-aver",
 398           FALSE,
 399           etBOOL,
 400           { &bPrAll },
 401           "Also print the exact average and rmsd stored in the energy frames (only when 1 term is "
 402           "requested)" },
 403         { "-orinst", FALSE, etBOOL, { &bOrinst }, "Analyse instantaneous orientation data" },
 404         { "-ovec", FALSE, etBOOL, { &bOvec }, "Also plot the eigenvectors with [TT]-oten[tt]" }
 405     };
 406     const char* drleg[] = { "Running average", "Instantaneous" };
 407
 408     FILE /* *out     = NULL,*/ *out_disre = nullptr, *fp_pairs = nullptr, *fort = nullptr,
 409                                *fodt = nullptr, *foten = nullptr;
 410     ener_file_t  fp;
 411     int          timecheck = 0;
 412     gmx_enxnm_t* enm       = nullptr;
 413     t_enxframe   fr;
 414     int          nre, teller, teller_disre;
 415     int          nor = 0, nex = 0, norfr = 0, enx_i = 0;
 416     real *bounds = nullptr, *violaver = nullptr, *oobs = nullptr, *orient = nullptr, *odrms = nullptr;
 417     int *       index = nullptr, *pair = nullptr, norsel = 0, *orsel = nullptr, *or_label = nullptr;
 418     int         nbounds = 0, npairs;
 419     gmx_bool    bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN;
 420     gmx_bool    bCont;
 421     double      sumaver, sumt;
 422     int *       set = nullptr, i, j, k, nset, sss;
 423     char **     pairleg, **odtleg, **otenleg;
 424     char**      leg = nullptr;
 425     const char *anm_j, *anm_k, *resnm_j, *resnm_k;
 426     int         resnr_j, resnr_k;
 427     const char* orinst_sub = "@ subtitle \"instantaneous\"\n";
 428     char        buf[256];
 429     gmx_output_env_t* oenv;
 430     t_enxblock*       blk_disre = nullptr;
 431     int               ndisre    = 0;
 432
 433     t_filenm fnm[] = { { efEDR, "-f", nullptr, ffREAD },
 434                        { efEDR, "-f2", nullptr, ffOPTRD },
 435                        { efTPR, "-s", nullptr, ffOPTRD },
 436                        // { efXVG, "-o",    "energy",  ffWRITE },
 437                        { efXVG, "-viol", "violaver", ffOPTWR },
 438                        { efXVG, "-pairs", "pairs", ffOPTWR },
 439                        { efXVG, "-ora", "orienta", ffOPTWR },
 440                        { efXVG, "-ort", "orientt", ffOPTWR },
 441                        { efXVG, "-oda", "orideva", ffOPTWR },
 442                        { efXVG, "-odr", "oridevr", ffOPTWR },
 443                        { efXVG, "-odt", "oridevt", ffOPTWR },
 444                        { efXVG, "-oten", "oriten", ffOPTWR } };
 445 #define NFILE asize(fnm)
 446     int npargs;
 447
 448     npargs = asize(pa);
 449     if (!parse_common_args(
 450                 &argc, argv, PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END, NFILE, fnm, npargs, pa, asize(desc), desc, 0, nullptr, &oenv))
 451     {
 452         return 0;
 453     }
 454
 455     bDRAll = opt2bSet("-pairs", NFILE, fnm);
 456     bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
 457     bORA   = opt2bSet("-ora", NFILE, fnm);
 458     bORT   = opt2bSet("-ort", NFILE, fnm);
 459     bODA   = opt2bSet("-oda", NFILE, fnm);
 460     bODR   = opt2bSet("-odr", NFILE, fnm);
 461     bODT   = opt2bSet("-odt", NFILE, fnm);
 462     bORIRE = bORA || bORT || bODA || bODR || bODT;
 463     bOTEN  = opt2bSet("-oten", NFILE, fnm);
 464     if (!(bDRAll || bDisRe || bORA || bORT || bODA || bODR || bODT || bORIRE || bOTEN))
 465     {
 466         printf("No output selected. Run with -h to see options. Terminating program.\n");
 467         return 0;
 468     }
 469     nset = 0;
 470
 471     fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
 472     do_enxnms(fp, &nre, &enm);
 473     free_enxnms(nre, enm);
 474
 475     t_inputrec  irInstance;
 476     t_inputrec* ir = &irInstance;
 477     init_enxframe(&fr);
 478     gmx::TopologyInformation        topInfo;
 479     std::unique_ptr<gmx_localtop_t> top;
 480     if (!bDisRe)
 481     {
 482         if (bORIRE || bOTEN)
 483         {
 484             get_orires_parms(ftp2fn(efTPR, NFILE, fnm), ir, &nor, &nex, &or_label, &oobs);
 485         }
 486
 487         if (bORIRE)
 488         {
 489             if (bOrinst)
 490             {
 491                 enx_i = enxORI;
 492             }
 493             else
 494             {
 495                 enx_i = enxOR;
 496             }
 497
 498             if (bORA || bODA)
 499             {
 500                 snew(orient, nor);
 501             }
 502             if (bODR)
 503             {
 504                 snew(odrms, nor);
 505             }
 506             if (bORT || bODT)
 507             {
 508                 fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n");
 509                 fprintf(stderr, "End your selection with 0\n");
 510                 j     = -1;
 511                 orsel = nullptr;
 512                 do
 513                 {
 514                     j++;
 515                     srenew(orsel, j + 1);
 516                     if (1 != scanf("%d", &(orsel[j])))
 517                     {
 518                         gmx_fatal(FARGS, "Error reading user input");
 519                     }
 520                 } while (orsel[j] > 0);
 521                 if (orsel[0] == -1)
 522                 {
 523                     fprintf(stderr, "Selecting all %d orientation restraints\n", nor);
 524                     norsel = nor;
 525                     srenew(orsel, nor);
 526                     for (i = 0; i < nor; i++)
 527                     {
 528                         orsel[i] = i;
 529                     }
 530                 }
 531                 else
 532                 {
 533                     /* Build the selection */
 534                     norsel = 0;
 535                     for (i = 0; i < j; i++)
 536                     {
 537                         for (k = 0; k < nor; k++)
 538                         {
 539                             if (or_label[k] == orsel[i])
 540                             {
 541                                 orsel[norsel] = k;
 542                                 norsel++;
 543                                 break;
 544                             }
 545                         }
 546                         if (k == nor)
 547                         {
 548                             fprintf(stderr, "Orientation restraint label %d not found\n", orsel[i]);
 549                         }
 550                     }
 551                 }
 552                 snew(odtleg, norsel);
 553                 for (i = 0; i < norsel; i++)
 554                 {
 555                     snew(odtleg[i], 256);
 556                     sprintf(odtleg[i], "%d", or_label[orsel[i]]);
 557                 }
 558                 if (bORT)
 559                 {
 560                     fort = xvgropen(
 561                             opt2fn("-ort", NFILE, fnm), "Calculated orientations", "Time (ps)", "", oenv);
 562                     if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 563                     {
 564                         fprintf(fort, "%s", orinst_sub);
 565                     }
 566                     xvgr_legend(fort, norsel, odtleg, oenv);
 567                 }
 568                 if (bODT)
 569                 {
 570                     fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
 571                                     "Orientation restraint deviation",
 572                                     "Time (ps)",
 573                                     "",
 574                                     oenv);
 575                     if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 576                     {
 577                         fprintf(fodt, "%s", orinst_sub);
 578                     }
 579                     xvgr_legend(fodt, norsel, odtleg, oenv);
 580                 }
 581                 for (i = 0; i < norsel; i++)
 582                 {
 583                     sfree(odtleg[i]);
 584                 }
 585                 sfree(odtleg);
 586             }
 587         }
 588         if (bOTEN)
 589         {
 590             foten = xvgropen(opt2fn("-oten", NFILE, fnm), "Order tensor", "Time (ps)", "", oenv);
 591             snew(otenleg, bOvec ? nex * 12 : nex * 3);
 592             for (i = 0; i < nex; i++)
 593             {
 594                 for (j = 0; j < 3; j++)
 595                 {
 596                     sprintf(buf, "eig%d", j + 1);
 597                     otenleg[(bOvec ? 12 : 3) * i + j] = gmx_strdup(buf);
 598                 }
 599                 if (bOvec)
 600                 {
 601                     for (j = 0; j < 9; j++)
 602                     {
 603                         sprintf(buf, "vec%d%s", j / 3 + 1, j % 3 == 0 ? "x" : (j % 3 == 1 ? "y" : "z"));
 604                         otenleg[12 * i + 3 + j] = gmx_strdup(buf);
 605                     }
 606                 }
 607             }
 608             xvgr_legend(foten, bOvec ? nex * 12 : nex * 3, otenleg, oenv);
 609             for (j = 0; j < 3; j++)
 610             {
 611                 sfree(otenleg[j]);
 612             }
 613             sfree(otenleg);
 614         }
 615     }
 616     else
 617     {
 618         {
 619             topInfo.fillFromInputFile(ftp2fn(efTPR, NFILE, fnm));
 620             top = std::make_unique<gmx_localtop_t>(topInfo.mtop()->ffparams);
 621             gmx_mtop_generate_local_top(*topInfo.mtop(), top.get(), ir->efep != efepNO);
 622         }
 623         nbounds = get_bounds(&bounds, &index, &pair, &npairs, top->idef);
 624         snew(violaver, npairs);
 625         out_disre = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations", "Time (ps)", "nm", oenv);
 626         xvgr_legend(out_disre, 2, drleg, oenv);
 627         if (bDRAll)
 628         {
 629             fp_pairs = xvgropen(
 630                     opt2fn("-pairs", NFILE, fnm), "Pair Distances", "Time (ps)", "Distance (nm)", oenv);
 631             if (output_env_get_print_xvgr_codes(oenv))
 632             {
 633                 fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n", ir->dr_tau);
 634             }
 635         }
 636     }
 637
 638     /* Initiate counters */
 639     teller       = 0;
 640     teller_disre = 0;
 641     do
 642     {
 643         /* This loop searches for the first frame (when -b option is given),
 644          * or when this has been found it reads just one energy frame
 645          */
 646         do
 647         {
 648             bCont = do_enx(fp, &fr);
 649             if (bCont)
 650             {
 651                 timecheck = check_times(fr.t);
 652             }
 653         } while (bCont && (timecheck < 0));
 654
 655         if ((timecheck == 0) && bCont)
 656         {
 657             /*
 658              * Define distance restraint legends. Can only be done after
 659              * the first frame has been read... (Then we know how many there are)
 660              */
 661             blk_disre = find_block_id_enxframe(&fr, enxDISRE, nullptr);
 662             if (bDisRe && bDRAll && !leg && blk_disre)
 663             {
 664                 const InteractionList&   ilist = top->idef.il[F_DISRES];
 665                 gmx::ArrayRef<const int> fa    = ilist.iatoms;
 666                 const t_iparams*         ip    = top->idef.iparams.data();
 667                 if (blk_disre->nsub != 2 || (blk_disre->sub[0].nr != blk_disre->sub[1].nr))
 668                 {
 669                     gmx_incons("Number of disre sub-blocks not equal to 2");
 670                 }
 671
 672                 ndisre = blk_disre->sub[0].nr;
 673                 if (ndisre != ilist.size() / 3)
 674                 {
 675                     gmx_fatal(FARGS,
 676                               "Number of disre pairs in the energy file (%d) does not match the "
 677                               "number in the run input file (%d)\n",
 678                               ndisre,
 679                               ilist.size() / 3);
 680                 }
 681                 snew(pairleg, ndisre);
 682                 int molb = 0;
 683                 for (i = 0; i < ndisre; i++)
 684                 {
 685                     snew(pairleg[i], 30);
 686                     j = fa[3 * i + 1];
 687                     k = fa[3 * i + 2];
 688                     mtopGetAtomAndResidueName(topInfo.mtop(), j, &molb, &anm_j, &resnr_j, &resnm_j, nullptr);
 689                     mtopGetAtomAndResidueName(topInfo.mtop(), k, &molb, &anm_k, &resnr_k, &resnm_k, nullptr);
 690                     sprintf(pairleg[i],
 691                             "%d %s %d %s (%d)",
 692                             resnr_j,
 693                             anm_j,
 694                             resnr_k,
 695                             anm_k,
 696                             ip[fa[3 * i]].disres.label);
 697                 }
 698                 set = select_it(ndisre, pairleg, &nset);
 699                 snew(leg, 2 * nset);
 700                 for (i = 0; (i < nset); i++)
 701                 {
 702                     snew(leg[2 * i], 32);
 703                     sprintf(leg[2 * i], "a %s", pairleg[set[i]]);
 704                     snew(leg[2 * i + 1], 32);
 705                     sprintf(leg[2 * i + 1], "i %s", pairleg[set[i]]);
 706                 }
 707                 xvgr_legend(fp_pairs, 2 * nset, leg, oenv);
 708             }
 709
 710             /*
 711              * Printing time, only when we do not want to skip frames
 712              */
 713             if (!skip || teller % skip == 0)
 714             {
 715                 if (bDisRe)
 716                 {
 717                     /*******************************************
 718                      * D I S T A N C E   R E S T R A I N T S
 719                      *******************************************/
 720                     if (ndisre > 0)
 721                     {
 722                         GMX_RELEASE_ASSERT(blk_disre != nullptr,
 723                                            "Trying to dereference NULL blk_disre pointer");
 724 #if !GMX_DOUBLE
 725                         float* disre_rt     = blk_disre->sub[0].fval;
 726                         float* disre_rm3tav = blk_disre->sub[1].fval;
 727 #else
 728                         double* disre_rt     = blk_disre->sub[0].dval;
 729                         double* disre_rm3tav = blk_disre->sub[1].dval;
 730 #endif
 731
 732                         print_time(out_disre, fr.t);
 733                         if (violaver == nullptr)
 734                         {
 735                             snew(violaver, ndisre);
 736                         }
 737
 738                         /* Subtract bounds from distances, to calculate violations */
 739                         calc_violations(
 740                                 disre_rt, disre_rm3tav, nbounds, pair, bounds, violaver, &sumt, &sumaver);
 741
 742                         fprintf(out_disre, "  %8.4f  %8.4f\n", sumaver, sumt);
 743                         if (bDRAll)
 744                         {
 745                             print_time(fp_pairs, fr.t);
 746                             for (i = 0; (i < nset); i++)
 747                             {
 748                                 sss = set[i];
 749                                 fprintf(fp_pairs, "  %8.4f", mypow(disre_rm3tav[sss], minthird));
 750                                 fprintf(fp_pairs, "  %8.4f", disre_rt[sss]);
 751                             }
 752                             fprintf(fp_pairs, "\n");
 753                         }
 754                         teller_disre++;
 755                     }
 756                 }
 757
 758                 /*******************************************
 759                  * O R I R E S
 760                  *******************************************/
 761                 else
 762                 {
 763                     t_enxblock* blk = find_block_id_enxframe(&fr, enx_i, nullptr);
 764                     if (bORIRE && blk)
 765                     {
 766                         if (blk->nsub != 1)
 767                         {
 768                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly.");
 769                         }
 770
 771                         if (blk->sub[0].nr != nor)
 772                         {
 773                             gmx_fatal(FARGS,
 774                                       "Number of orientation restraints in energy file (%d) does "
 775                                       "not match with the topology (%d)",
 776                                       blk->sub[0].nr,
 777                                       nor);
 778                         }
 779                         if (bORA || bODA)
 780                         {
 781                             for (i = 0; i < nor; i++)
 782                             {
 783                                 orient[i] += blk_value(blk, 0, i);
 784                             }
 785                         }
 786                         if (bODR)
 787                         {
 788                             for (i = 0; i < nor; i++)
 789                             {
 790                                 real v = blk_value(blk, 0, i);
 791                                 odrms[i] += gmx::square(v - oobs[i]);
 792                             }
 793                         }
 794                         if (bORT)
 795                         {
 796                             fprintf(fort, "  %10f", fr.t);
 797                             for (i = 0; i < norsel; i++)
 798                             {
 799                                 fprintf(fort, " %g", blk_value(blk, 0, orsel[i]));
 800                             }
 801                             fprintf(fort, "\n");
 802                         }
 803                         if (bODT)
 804                         {
 805                             fprintf(fodt, "  %10f", fr.t);
 806                             for (i = 0; i < norsel; i++)
 807                             {
 808                                 fprintf(fodt, " %g", blk_value(blk, 0, orsel[i]) - oobs[orsel[i]]);
 809                             }
 810                             fprintf(fodt, "\n");
 811                         }
 812                         norfr++;
 813                     }
 814                     blk = find_block_id_enxframe(&fr, enxORT, nullptr);
 815                     if (bOTEN && blk)
 816                     {
 817                         if (blk->nsub != 1)
 818                         {
 819                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
 820                         }
 821
 822                         if (blk->sub[0].nr != nex * 12)
 823                         {
 824                             gmx_fatal(FARGS,
 825                                       "Number of orientation experiments in energy file (%d) does "
 826                                       "not match with the topology (%d)",
 827                                       blk->sub[0].nr / 12,
 828                                       nex);
 829                         }
 830                         fprintf(foten, "  %10f", fr.t);
 831                         for (i = 0; i < nex; i++)
 832                         {
 833                             for (j = 0; j < (bOvec ? 12 : 3); j++)
 834                             {
 835                                 fprintf(foten, " %g", blk_value(blk, 0, i * 12 + j));
 836                             }
 837                         }
 838                         fprintf(foten, "\n");
 839                     }
 840                 }
 841             }
 842             teller++;
 843         }
 844     } while (bCont && (timecheck == 0));
 845     free_enxframe(&fr);
 846
 847     fprintf(stderr, "\n");
 848     done_ener_file(fp);
 849     if (out_disre)
 850     {
 851         xvgrclose(out_disre);
 852     }
 853
 854     if (bDRAll)
 855     {
 856         xvgrclose(fp_pairs);
 857     }
 858
 859     if (bORT)
 860     {
 861         xvgrclose(fort);
 862     }
 863     if (bODT)
 864     {
 865         xvgrclose(fodt);
 866     }
 867     if (bORA)
 868     {
 869         FILE* out = xvgropen(
 870                 opt2fn("-ora", NFILE, fnm), "Average calculated orientations", "Restraint label", "", oenv);
 871         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 872         {
 873             fprintf(out, "%s", orinst_sub);
 874         }
 875         for (i = 0; i < nor; i++)
 876         {
 877             fprintf(out, "%5d  %g\n", or_label[i], orient[i] / real(norfr));
 878         }
 879         xvgrclose(out);
 880     }
 881     if (bODA)
 882     {
 883         FILE* out = xvgropen(
 884                 opt2fn("-oda", NFILE, fnm), "Average restraint deviation", "Restraint label", "", oenv);
 885         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 886         {
 887             fprintf(out, "%s", orinst_sub);
 888         }
 889         for (i = 0; i < nor; i++)
 890         {
 891             fprintf(out, "%5d  %g\n", or_label[i], orient[i] / real(norfr) - oobs[i]);
 892         }
 893         xvgrclose(out);
 894     }
 895     if (bODR)
 896     {
 897         FILE* out = xvgropen(opt2fn("-odr", NFILE, fnm),
 898                              "RMS orientation restraint deviations",
 899                              "Restraint label",
 900                              "",
 901                              oenv);
 902         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 903         {
 904             fprintf(out, "%s", orinst_sub);
 905         }
 906         for (i = 0; i < nor; i++)
 907         {
 908             fprintf(out, "%5d  %g\n", or_label[i], std::sqrt(odrms[i] / real(norfr)));
 909         }
 910         xvgrclose(out);
 911     }
 912     // Clean up orires variables.
 913     sfree(or_label);
 914     sfree(oobs);
 915     sfree(orient);
 916     sfree(odrms);
 917     sfree(orsel);
 918     if (bOTEN)
 919     {
 920         xvgrclose(foten);
 921     }
 922
 923     if (bDisRe)
 924     {
 925         analyse_disre(opt2fn("-viol", NFILE, fnm), teller_disre, violaver, bounds, index, pair, nbounds, oenv);
 926     }
 927     {
 928         const char* nxy = "-nxy";
 929
 930         do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
 931         do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
 932         do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
 933         do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
 934         do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
 935         do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
 936     }
 937     output_env_done(oenv);
 938
 939     return 0;
 940 }