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                     GMX_ASSERT(topInfo.hasTopology(), "Need to have a valid topology");
 690                     mtopGetAtomAndResidueName(*topInfo.mtop(), j, &molb, &anm_j, &resnr_j, &resnm_j, nullptr);
 691                     mtopGetAtomAndResidueName(*topInfo.mtop(), k, &molb, &anm_k, &resnr_k, &resnm_k, nullptr);
 692                     sprintf(pairleg[i],
 693                             "%d %s %d %s (%d)",
 694                             resnr_j,
 695                             anm_j,
 696                             resnr_k,
 697                             anm_k,
 698                             ip[fa[3 * i]].disres.label);
 699                 }
 700                 set = select_it(ndisre, pairleg, &nset);
 701                 snew(leg, 2 * nset);
 702                 for (i = 0; (i < nset); i++)
 703                 {
 704                     snew(leg[2 * i], 32);
 705                     sprintf(leg[2 * i], "a %s", pairleg[set[i]]);
 706                     snew(leg[2 * i + 1], 32);
 707                     sprintf(leg[2 * i + 1], "i %s", pairleg[set[i]]);
 708                 }
 709                 xvgr_legend(fp_pairs, 2 * nset, leg, oenv);
 710             }
 711
 712             /*
 713              * Printing time, only when we do not want to skip frames
 714              */
 715             if (!skip || teller % skip == 0)
 716             {
 717                 if (bDisRe)
 718                 {
 719                     /*******************************************
 720                      * D I S T A N C E   R E S T R A I N T S
 721                      *******************************************/
 722                     if (ndisre > 0)
 723                     {
 724                         GMX_RELEASE_ASSERT(blk_disre != nullptr,
 725                                            "Trying to dereference NULL blk_disre pointer");
 726 #if !GMX_DOUBLE
 727                         float* disre_rt     = blk_disre->sub[0].fval;
 728                         float* disre_rm3tav = blk_disre->sub[1].fval;
 729 #else
 730                         double* disre_rt     = blk_disre->sub[0].dval;
 731                         double* disre_rm3tav = blk_disre->sub[1].dval;
 732 #endif
 733
 734                         print_time(out_disre, fr.t);
 735                         if (violaver == nullptr)
 736                         {
 737                             snew(violaver, ndisre);
 738                         }
 739
 740                         /* Subtract bounds from distances, to calculate violations */
 741                         calc_violations(
 742                                 disre_rt, disre_rm3tav, nbounds, pair, bounds, violaver, &sumt, &sumaver);
 743
 744                         fprintf(out_disre, "  %8.4f  %8.4f\n", sumaver, sumt);
 745                         if (bDRAll)
 746                         {
 747                             print_time(fp_pairs, fr.t);
 748                             for (i = 0; (i < nset); i++)
 749                             {
 750                                 sss = set[i];
 751                                 fprintf(fp_pairs, "  %8.4f", mypow(disre_rm3tav[sss], minthird));
 752                                 fprintf(fp_pairs, "  %8.4f", disre_rt[sss]);
 753                             }
 754                             fprintf(fp_pairs, "\n");
 755                         }
 756                         teller_disre++;
 757                     }
 758                 }
 759
 760                 /*******************************************
 761                  * O R I R E S
 762                  *******************************************/
 763                 else
 764                 {
 765                     t_enxblock* blk = find_block_id_enxframe(&fr, enx_i, nullptr);
 766                     if (bORIRE && blk)
 767                     {
 768                         if (blk->nsub != 1)
 769                         {
 770                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly.");
 771                         }
 772
 773                         if (blk->sub[0].nr != nor)
 774                         {
 775                             gmx_fatal(FARGS,
 776                                       "Number of orientation restraints in energy file (%d) does "
 777                                       "not match with the topology (%d)",
 778                                       blk->sub[0].nr,
 779                                       nor);
 780                         }
 781                         if (bORA || bODA)
 782                         {
 783                             for (i = 0; i < nor; i++)
 784                             {
 785                                 orient[i] += blk_value(blk, 0, i);
 786                             }
 787                         }
 788                         if (bODR)
 789                         {
 790                             for (i = 0; i < nor; i++)
 791                             {
 792                                 real v = blk_value(blk, 0, i);
 793                                 odrms[i] += gmx::square(v - oobs[i]);
 794                             }
 795                         }
 796                         if (bORT)
 797                         {
 798                             fprintf(fort, "  %10f", fr.t);
 799                             for (i = 0; i < norsel; i++)
 800                             {
 801                                 fprintf(fort, " %g", blk_value(blk, 0, orsel[i]));
 802                             }
 803                             fprintf(fort, "\n");
 804                         }
 805                         if (bODT)
 806                         {
 807                             fprintf(fodt, "  %10f", fr.t);
 808                             for (i = 0; i < norsel; i++)
 809                             {
 810                                 fprintf(fodt, " %g", blk_value(blk, 0, orsel[i]) - oobs[orsel[i]]);
 811                             }
 812                             fprintf(fodt, "\n");
 813                         }
 814                         norfr++;
 815                     }
 816                     blk = find_block_id_enxframe(&fr, enxORT, nullptr);
 817                     if (bOTEN && blk)
 818                     {
 819                         if (blk->nsub != 1)
 820                         {
 821                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
 822                         }
 823
 824                         if (blk->sub[0].nr != nex * 12)
 825                         {
 826                             gmx_fatal(FARGS,
 827                                       "Number of orientation experiments in energy file (%d) does "
 828                                       "not match with the topology (%d)",
 829                                       blk->sub[0].nr / 12,
 830                                       nex);
 831                         }
 832                         fprintf(foten, "  %10f", fr.t);
 833                         for (i = 0; i < nex; i++)
 834                         {
 835                             for (j = 0; j < (bOvec ? 12 : 3); j++)
 836                             {
 837                                 fprintf(foten, " %g", blk_value(blk, 0, i * 12 + j));
 838                             }
 839                         }
 840                         fprintf(foten, "\n");
 841                     }
 842                 }
 843             }
 844             teller++;
 845         }
 846     } while (bCont && (timecheck == 0));
 847     free_enxframe(&fr);
 848
 849     fprintf(stderr, "\n");
 850     done_ener_file(fp);
 851     if (out_disre)
 852     {
 853         xvgrclose(out_disre);
 854     }
 855
 856     if (bDRAll)
 857     {
 858         xvgrclose(fp_pairs);
 859     }
 860
 861     if (bORT)
 862     {
 863         xvgrclose(fort);
 864     }
 865     if (bODT)
 866     {
 867         xvgrclose(fodt);
 868     }
 869     if (bORA)
 870     {
 871         FILE* out = xvgropen(
 872                 opt2fn("-ora", NFILE, fnm), "Average calculated orientations", "Restraint label", "", oenv);
 873         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 874         {
 875             fprintf(out, "%s", orinst_sub);
 876         }
 877         for (i = 0; i < nor; i++)
 878         {
 879             fprintf(out, "%5d  %g\n", or_label[i], orient[i] / real(norfr));
 880         }
 881         xvgrclose(out);
 882     }
 883     if (bODA)
 884     {
 885         FILE* out = xvgropen(
 886                 opt2fn("-oda", NFILE, fnm), "Average restraint deviation", "Restraint label", "", oenv);
 887         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 888         {
 889             fprintf(out, "%s", orinst_sub);
 890         }
 891         for (i = 0; i < nor; i++)
 892         {
 893             fprintf(out, "%5d  %g\n", or_label[i], orient[i] / real(norfr) - oobs[i]);
 894         }
 895         xvgrclose(out);
 896     }
 897     if (bODR)
 898     {
 899         FILE* out = xvgropen(opt2fn("-odr", NFILE, fnm),
 900                              "RMS orientation restraint deviations",
 901                              "Restraint label",
 902                              "",
 903                              oenv);
 904         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 905         {
 906             fprintf(out, "%s", orinst_sub);
 907         }
 908         for (i = 0; i < nor; i++)
 909         {
 910             fprintf(out, "%5d  %g\n", or_label[i], std::sqrt(odrms[i] / real(norfr)));
 911         }
 912         xvgrclose(out);
 913     }
 914     // Clean up orires variables.
 915     sfree(or_label);
 916     sfree(oobs);
 917     sfree(orient);
 918     sfree(odrms);
 919     sfree(orsel);
 920     if (bOTEN)
 921     {
 922         xvgrclose(foten);
 923     }
 924
 925     if (bDisRe)
 926     {
 927         analyse_disre(opt2fn("-viol", NFILE, fnm), teller_disre, violaver, bounds, index, pair, nbounds, oenv);
 928     }
 929     {
 930         const char* nxy = "-nxy";
 931
 932         do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
 933         do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
 934         do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
 935         do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
 936         do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
 937         do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
 938     }
 939     output_env_done(oenv);
 940
 941     return 0;
 942 }