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,2021, 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(
 622                     *topInfo.mtop(), top.get(), ir->efep != FreeEnergyPerturbationType::No);
 623         }
 624         nbounds = get_bounds(&bounds, &index, &pair, &npairs, top->idef);
 625         snew(violaver, npairs);
 626         out_disre = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations", "Time (ps)", "nm", oenv);
 627         xvgr_legend(out_disre, 2, drleg, oenv);
 628         if (bDRAll)
 629         {
 630             fp_pairs = xvgropen(
 631                     opt2fn("-pairs", NFILE, fnm), "Pair Distances", "Time (ps)", "Distance (nm)", oenv);
 632             if (output_env_get_print_xvgr_codes(oenv))
 633             {
 634                 fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n", ir->dr_tau);
 635             }
 636         }
 637     }
 638
 639     /* Initiate counters */
 640     teller       = 0;
 641     teller_disre = 0;
 642     do
 643     {
 644         /* This loop searches for the first frame (when -b option is given),
 645          * or when this has been found it reads just one energy frame
 646          */
 647         do
 648         {
 649             bCont = do_enx(fp, &fr);
 650             if (bCont)
 651             {
 652                 timecheck = check_times(fr.t);
 653             }
 654         } while (bCont && (timecheck < 0));
 655
 656         if ((timecheck == 0) && bCont)
 657         {
 658             /*
 659              * Define distance restraint legends. Can only be done after
 660              * the first frame has been read... (Then we know how many there are)
 661              */
 662             blk_disre = find_block_id_enxframe(&fr, enxDISRE, nullptr);
 663             if (bDisRe && bDRAll && !leg && blk_disre)
 664             {
 665                 const InteractionList&   ilist = top->idef.il[F_DISRES];
 666                 gmx::ArrayRef<const int> fa    = ilist.iatoms;
 667                 const t_iparams*         ip    = top->idef.iparams.data();
 668                 if (blk_disre->nsub != 2 || (blk_disre->sub[0].nr != blk_disre->sub[1].nr))
 669                 {
 670                     gmx_incons("Number of disre sub-blocks not equal to 2");
 671                 }
 672
 673                 ndisre = blk_disre->sub[0].nr;
 674                 if (ndisre != ilist.size() / 3)
 675                 {
 676                     gmx_fatal(FARGS,
 677                               "Number of disre pairs in the energy file (%d) does not match the "
 678                               "number in the run input file (%d)\n",
 679                               ndisre,
 680                               ilist.size() / 3);
 681                 }
 682                 snew(pairleg, ndisre);
 683                 int molb = 0;
 684                 for (i = 0; i < ndisre; i++)
 685                 {
 686                     snew(pairleg[i], 30);
 687                     j = fa[3 * i + 1];
 688                     k = fa[3 * i + 2];
 689                     mtopGetAtomAndResidueName(topInfo.mtop(), j, &molb, &anm_j, &resnr_j, &resnm_j, nullptr);
 690                     mtopGetAtomAndResidueName(topInfo.mtop(), k, &molb, &anm_k, &resnr_k, &resnm_k, nullptr);
 691                     sprintf(pairleg[i],
 692                             "%d %s %d %s (%d)",
 693                             resnr_j,
 694                             anm_j,
 695                             resnr_k,
 696                             anm_k,
 697                             ip[fa[3 * i]].disres.label);
 698                 }
 699                 set = select_it(ndisre, pairleg, &nset);
 700                 snew(leg, 2 * nset);
 701                 for (i = 0; (i < nset); i++)
 702                 {
 703                     snew(leg[2 * i], 32);
 704                     sprintf(leg[2 * i], "a %s", pairleg[set[i]]);
 705                     snew(leg[2 * i + 1], 32);
 706                     sprintf(leg[2 * i + 1], "i %s", pairleg[set[i]]);
 707                 }
 708                 xvgr_legend(fp_pairs, 2 * nset, leg, oenv);
 709             }
 710
 711             /*
 712              * Printing time, only when we do not want to skip frames
 713              */
 714             if (!skip || teller % skip == 0)
 715             {
 716                 if (bDisRe)
 717                 {
 718                     /*******************************************
 719                      * D I S T A N C E   R E S T R A I N T S
 720                      *******************************************/
 721                     if (ndisre > 0)
 722                     {
 723                         GMX_RELEASE_ASSERT(blk_disre != nullptr,
 724                                            "Trying to dereference NULL blk_disre pointer");
 725 #if !GMX_DOUBLE
 726                         float* disre_rt     = blk_disre->sub[0].fval;
 727                         float* disre_rm3tav = blk_disre->sub[1].fval;
 728 #else
 729                         double* disre_rt     = blk_disre->sub[0].dval;
 730                         double* disre_rm3tav = blk_disre->sub[1].dval;
 731 #endif
 732
 733                         print_time(out_disre, fr.t);
 734                         if (violaver == nullptr)
 735                         {
 736                             snew(violaver, ndisre);
 737                         }
 738
 739                         /* Subtract bounds from distances, to calculate violations */
 740                         calc_violations(
 741                                 disre_rt, disre_rm3tav, nbounds, pair, bounds, violaver, &sumt, &sumaver);
 742
 743                         fprintf(out_disre, "  %8.4f  %8.4f\n", sumaver, sumt);
 744                         if (bDRAll)
 745                         {
 746                             print_time(fp_pairs, fr.t);
 747                             for (i = 0; (i < nset); i++)
 748                             {
 749                                 sss = set[i];
 750                                 fprintf(fp_pairs, "  %8.4f", mypow(disre_rm3tav[sss], minthird));
 751                                 fprintf(fp_pairs, "  %8.4f", disre_rt[sss]);
 752                             }
 753                             fprintf(fp_pairs, "\n");
 754                         }
 755                         teller_disre++;
 756                     }
 757                 }
 758
 759                 /*******************************************
 760                  * O R I R E S
 761                  *******************************************/
 762                 else
 763                 {
 764                     t_enxblock* blk = find_block_id_enxframe(&fr, enx_i, nullptr);
 765                     if (bORIRE && blk)
 766                     {
 767                         if (blk->nsub != 1)
 768                         {
 769                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly.");
 770                         }
 771
 772                         if (blk->sub[0].nr != nor)
 773                         {
 774                             gmx_fatal(FARGS,
 775                                       "Number of orientation restraints in energy file (%d) does "
 776                                       "not match with the topology (%d)",
 777                                       blk->sub[0].nr,
 778                                       nor);
 779                         }
 780                         if (bORA || bODA)
 781                         {
 782                             for (i = 0; i < nor; i++)
 783                             {
 784                                 orient[i] += blk_value(blk, 0, i);
 785                             }
 786                         }
 787                         if (bODR)
 788                         {
 789                             for (i = 0; i < nor; i++)
 790                             {
 791                                 real v = blk_value(blk, 0, i);
 792                                 odrms[i] += gmx::square(v - oobs[i]);
 793                             }
 794                         }
 795                         if (bORT)
 796                         {
 797                             fprintf(fort, "  %10f", fr.t);
 798                             for (i = 0; i < norsel; i++)
 799                             {
 800                                 fprintf(fort, " %g", blk_value(blk, 0, orsel[i]));
 801                             }
 802                             fprintf(fort, "\n");
 803                         }
 804                         if (bODT)
 805                         {
 806                             fprintf(fodt, "  %10f", fr.t);
 807                             for (i = 0; i < norsel; i++)
 808                             {
 809                                 fprintf(fodt, " %g", blk_value(blk, 0, orsel[i]) - oobs[orsel[i]]);
 810                             }
 811                             fprintf(fodt, "\n");
 812                         }
 813                         norfr++;
 814                     }
 815                     blk = find_block_id_enxframe(&fr, enxORT, nullptr);
 816                     if (bOTEN && blk)
 817                     {
 818                         if (blk->nsub != 1)
 819                         {
 820                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
 821                         }
 822
 823                         if (blk->sub[0].nr != nex * 12)
 824                         {
 825                             gmx_fatal(FARGS,
 826                                       "Number of orientation experiments in energy file (%d) does "
 827                                       "not match with the topology (%d)",
 828                                       blk->sub[0].nr / 12,
 829                                       nex);
 830                         }
 831                         fprintf(foten, "  %10f", fr.t);
 832                         for (i = 0; i < nex; i++)
 833                         {
 834                             for (j = 0; j < (bOvec ? 12 : 3); j++)
 835                             {
 836                                 fprintf(foten, " %g", blk_value(blk, 0, i * 12 + j));
 837                             }
 838                         }
 839                         fprintf(foten, "\n");
 840                     }
 841                 }
 842             }
 843             teller++;
 844         }
 845     } while (bCont && (timecheck == 0));
 846     free_enxframe(&fr);
 847
 848     fprintf(stderr, "\n");
 849     done_ener_file(fp);
 850     if (out_disre)
 851     {
 852         xvgrclose(out_disre);
 853     }
 854
 855     if (bDRAll)
 856     {
 857         xvgrclose(fp_pairs);
 858     }
 859
 860     if (bORT)
 861     {
 862         xvgrclose(fort);
 863     }
 864     if (bODT)
 865     {
 866         xvgrclose(fodt);
 867     }
 868     if (bORA)
 869     {
 870         FILE* out = xvgropen(
 871                 opt2fn("-ora", NFILE, fnm), "Average calculated orientations", "Restraint label", "", oenv);
 872         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 873         {
 874             fprintf(out, "%s", orinst_sub);
 875         }
 876         for (i = 0; i < nor; i++)
 877         {
 878             fprintf(out, "%5d  %g\n", or_label[i], orient[i] / real(norfr));
 879         }
 880         xvgrclose(out);
 881     }
 882     if (bODA)
 883     {
 884         FILE* out = xvgropen(
 885                 opt2fn("-oda", NFILE, fnm), "Average restraint deviation", "Restraint label", "", oenv);
 886         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 887         {
 888             fprintf(out, "%s", orinst_sub);
 889         }
 890         for (i = 0; i < nor; i++)
 891         {
 892             fprintf(out, "%5d  %g\n", or_label[i], orient[i] / real(norfr) - oobs[i]);
 893         }
 894         xvgrclose(out);
 895     }
 896     if (bODR)
 897     {
 898         FILE* out = xvgropen(opt2fn("-odr", NFILE, fnm),
 899                              "RMS orientation restraint deviations",
 900                              "Restraint label",
 901                              "",
 902                              oenv);
 903         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 904         {
 905             fprintf(out, "%s", orinst_sub);
 906         }
 907         for (i = 0; i < nor; i++)
 908         {
 909             fprintf(out, "%5d  %g\n", or_label[i], std::sqrt(odrms[i] / real(norfr)));
 910         }
 911         xvgrclose(out);
 912     }
 913     // Clean up orires variables.
 914     sfree(or_label);
 915     sfree(oobs);
 916     sfree(orient);
 917     sfree(odrms);
 918     sfree(orsel);
 919     if (bOTEN)
 920     {
 921         xvgrclose(foten);
 922     }
 923
 924     if (bDisRe)
 925     {
 926         analyse_disre(opt2fn("-viol", NFILE, fnm), teller_disre, violaver, bounds, index, pair, nbounds, oenv);
 927     }
 928     {
 929         const char* nxy = "-nxy";
 930
 931         do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
 932         do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
 933         do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
 934         do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
 935         do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
 936         do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
 937     }
 938     output_env_done(oenv);
 939
 940     return 0;
 941 }