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, gmx_localtop_t* top)
 202 {
 203     t_functype* functype;
 204     t_iparams*  ip;
 205     int         i, j, k, type, ftype, natom;
 206     t_ilist*    disres;
 207     t_iatom*    iatom;
 208     real*       b;
 209     int *       ind, *pair;
 210     int         nb, label1;
 211
 212     functype = top->idef.functype;
 213     ip       = top->idef.iparams;
 214
 215     /* Count how many distance restraint there are... */
 216     nb = top->idef.il[F_DISRES].nr;
 217     if (nb == 0)
 218     {
 219         gmx_fatal(FARGS, "No distance restraints in topology!\n");
 220     }
 221
 222     /* Allocate memory */
 223     snew(b, nb);
 224     snew(ind, nb);
 225     snew(pair, nb + 1);
 226
 227     /* Fill the bound array */
 228     nb = 0;
 229     for (i = 0; (i < top->idef.ntypes); i++)
 230     {
 231         ftype = functype[i];
 232         if (ftype == F_DISRES)
 233         {
 234
 235             label1  = ip[i].disres.label;
 236             b[nb]   = ip[i].disres.up1;
 237             ind[nb] = label1;
 238             nb++;
 239         }
 240     }
 241     *bounds = b;
 242
 243     /* Fill the index array */
 244     label1 = -1;
 245     disres = &(top->idef.il[F_DISRES]);
 246     iatom  = disres->iatoms;
 247     for (i = j = k = 0; (i < disres->nr);)
 248     {
 249         type  = iatom[i];
 250         ftype = top->idef.functype[type];
 251         natom = interaction_function[ftype].nratoms + 1;
 252         if (label1 != top->idef.iparams[type].disres.label)
 253         {
 254             pair[j] = k;
 255             label1  = top->idef.iparams[type].disres.label;
 256             j++;
 257         }
 258         k++;
 259         i += natom;
 260     }
 261     pair[j] = k;
 262     *npairs = k;
 263     if (j != nb)
 264     {
 265         gmx_incons("get_bounds for distance restraints");
 266     }
 267
 268     *index   = ind;
 269     *dr_pair = pair;
 270
 271     return nb;
 272 }
 273
 274 static void
 275 calc_violations(real rt[], real rav3[], int nb, const int index[], real bounds[], real* viol, double* st, double* sa)
 276 {
 277     const real sixth = 1.0 / 6.0;
 278     int        i, j;
 279     double     rsum, rav, sumaver, sumt;
 280
 281     sumaver = 0;
 282     sumt    = 0;
 283     for (i = 0; (i < nb); i++)
 284     {
 285         rsum = 0.0;
 286         rav  = 0.0;
 287         for (j = index[i]; (j < index[i + 1]); j++)
 288         {
 289             if (viol)
 290             {
 291                 viol[j] += mypow(rt[j], -3.0);
 292             }
 293             rav += gmx::square(rav3[j]);
 294             rsum += mypow(rt[j], -6);
 295         }
 296         rsum = std::max(0.0, mypow(rsum, -sixth) - bounds[i]);
 297         rav  = std::max(0.0, mypow(rav, -sixth) - bounds[i]);
 298
 299         sumt += rsum;
 300         sumaver += rav;
 301     }
 302     *st = sumt;
 303     *sa = sumaver;
 304 }
 305
 306 static void analyse_disre(const char*             voutfn,
 307                           int                     nframes,
 308                           real                    violaver[],
 309                           real                    bounds[],
 310                           int                     index[],
 311                           int                     pair[],
 312                           int                     nbounds,
 313                           const gmx_output_env_t* oenv)
 314 {
 315     FILE*  vout;
 316     double sum, sumt, sumaver;
 317     int    i, j;
 318
 319     /* Subtract bounds from distances, to calculate violations */
 320     calc_violations(violaver, violaver, nbounds, pair, bounds, nullptr, &sumt, &sumaver);
 321
 322 #ifdef DEBUG
 323     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n", sumaver);
 324     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sumt);
 325 #endif
 326     vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm", oenv);
 327     sum  = 0.0;
 328     sumt = 0.0;
 329     for (i = 0; (i < nbounds); i++)
 330     {
 331         /* Do ensemble averaging */
 332         sumaver = 0;
 333         for (j = pair[i]; (j < pair[i + 1]); j++)
 334         {
 335             sumaver += gmx::square(violaver[j] / real(nframes));
 336         }
 337         sumaver = std::max(0.0, mypow(sumaver, minsixth) - bounds[i]);
 338
 339         sumt += sumaver;
 340         sum = std::max(sum, sumaver);
 341         fprintf(vout, "%10d  %10.5e\n", index[i], sumaver);
 342     }
 343 #ifdef DEBUG
 344     for (j = 0; (j < dr.ndr); j++)
 345     {
 346         fprintf(vout, "%10d  %10.5e\n", j, mypow(violaver[j] / real(nframes), minthird));
 347     }
 348 #endif
 349     xvgrclose(vout);
 350
 351     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n", sumt);
 352     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
 353
 354     do_view(oenv, voutfn, "-graphtype bar");
 355 }
 356
 357 static void print_time(FILE* fp, double t)
 358 {
 359     fprintf(fp, "%12.6f", t);
 360 }
 361
 362 int gmx_nmr(int argc, char* argv[])
 363 {
 364     const char* desc[] = {
 365         "[THISMODULE] extracts distance or orientation restraint",
 366         "data from an energy file. The user is prompted to interactively",
 367         "select the desired terms.[PAR]",
 368
 369         "When the [TT]-viol[tt] option is set, the time averaged",
 370         "violations are plotted and the running time-averaged and",
 371         "instantaneous sum of violations are recalculated. Additionally",
 372         "running time-averaged and instantaneous distances between",
 373         "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
 374
 375         "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
 376         "[TT]-odt[tt] are used for analyzing orientation restraint data.",
 377         "The first two options plot the orientation, the last three the",
 378         "deviations of the orientations from the experimental values.",
 379         "The options that end on an 'a' plot the average over time",
 380         "as a function of restraint. The options that end on a 't'",
 381         "prompt the user for restraint label numbers and plot the data",
 382         "as a function of time. Option [TT]-odr[tt] plots the RMS",
 383         "deviation as a function of restraint.",
 384         "When the run used time or ensemble averaged orientation restraints,",
 385         "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
 386         "not ensemble-averaged orientations and deviations instead of",
 387         "the time and ensemble averages.[PAR]",
 388
 389         "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
 390         "tensor for each orientation restraint experiment. With option",
 391         "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
 392
 393
 394     };
 395     static gmx_bool bPrAll = FALSE;
 396     static gmx_bool bDp = FALSE, bOrinst = FALSE, bOvec = FALSE;
 397     static int      skip = 0;
 398     t_pargs         pa[] = {
 399         { "-dp", FALSE, etBOOL, { &bDp }, "Print energies in high precision" },
 400         { "-skip", FALSE, etINT, { &skip }, "Skip number of frames between data points" },
 401         { "-aver",
 402           FALSE,
 403           etBOOL,
 404           { &bPrAll },
 405           "Also print the exact average and rmsd stored in the energy frames (only when 1 term is "
 406           "requested)" },
 407         { "-orinst", FALSE, etBOOL, { &bOrinst }, "Analyse instantaneous orientation data" },
 408         { "-ovec", FALSE, etBOOL, { &bOvec }, "Also plot the eigenvectors with [TT]-oten[tt]" }
 409     };
 410     const char* drleg[] = { "Running average", "Instantaneous" };
 411
 412     FILE /* *out     = NULL,*/ *out_disre = nullptr, *fp_pairs = nullptr, *fort = nullptr,
 413                                *fodt = nullptr, *foten = nullptr;
 414     ener_file_t    fp;
 415     int            timecheck = 0;
 416     gmx_localtop_t top;
 417     gmx_enxnm_t*   enm = nullptr;
 418     t_enxframe     fr;
 419     int            nre, teller, teller_disre;
 420     int            nor = 0, nex = 0, norfr = 0, enx_i = 0;
 421     real *bounds = nullptr, *violaver = nullptr, *oobs = nullptr, *orient = nullptr, *odrms = nullptr;
 422     int *       index = nullptr, *pair = nullptr, norsel = 0, *orsel = nullptr, *or_label = nullptr;
 423     int         nbounds = 0, npairs;
 424     gmx_bool    bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN;
 425     gmx_bool    bCont;
 426     double      sumaver, sumt;
 427     int *       set = nullptr, i, j, k, nset, sss;
 428     char **     pairleg, **odtleg, **otenleg;
 429     char**      leg = nullptr;
 430     const char *anm_j, *anm_k, *resnm_j, *resnm_k;
 431     int         resnr_j, resnr_k;
 432     const char* orinst_sub = "@ subtitle \"instantaneous\"\n";
 433     char        buf[256];
 434     gmx_output_env_t* oenv;
 435     t_enxblock*       blk_disre = nullptr;
 436     int               ndisre    = 0;
 437
 438     t_filenm fnm[] = { { efEDR, "-f", nullptr, ffREAD },
 439                        { efEDR, "-f2", nullptr, ffOPTRD },
 440                        { efTPR, "-s", nullptr, ffOPTRD },
 441                        // { efXVG, "-o",    "energy",  ffWRITE },
 442                        { efXVG, "-viol", "violaver", ffOPTWR },
 443                        { efXVG, "-pairs", "pairs", ffOPTWR },
 444                        { efXVG, "-ora", "orienta", ffOPTWR },
 445                        { efXVG, "-ort", "orientt", ffOPTWR },
 446                        { efXVG, "-oda", "orideva", ffOPTWR },
 447                        { efXVG, "-odr", "oridevr", ffOPTWR },
 448                        { efXVG, "-odt", "oridevt", ffOPTWR },
 449                        { efXVG, "-oten", "oriten", ffOPTWR } };
 450 #define NFILE asize(fnm)
 451     int npargs;
 452
 453     npargs = asize(pa);
 454     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END, NFILE, fnm,
 455                            npargs, pa, asize(desc), desc, 0, nullptr, &oenv))
 456     {
 457         return 0;
 458     }
 459
 460     bDRAll = opt2bSet("-pairs", NFILE, fnm);
 461     bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
 462     bORA   = opt2bSet("-ora", NFILE, fnm);
 463     bORT   = opt2bSet("-ort", NFILE, fnm);
 464     bODA   = opt2bSet("-oda", NFILE, fnm);
 465     bODR   = opt2bSet("-odr", NFILE, fnm);
 466     bODT   = opt2bSet("-odt", NFILE, fnm);
 467     bORIRE = bORA || bORT || bODA || bODR || bODT;
 468     bOTEN  = opt2bSet("-oten", NFILE, fnm);
 469     if (!(bDRAll || bDisRe || bORA || bORT || bODA || bODR || bODT || bORIRE || bOTEN))
 470     {
 471         printf("No output selected. Run with -h to see options. Terminating program.\n");
 472         return 0;
 473     }
 474     nset = 0;
 475
 476     fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
 477     do_enxnms(fp, &nre, &enm);
 478     free_enxnms(nre, enm);
 479
 480     t_inputrec  irInstance;
 481     t_inputrec* ir = &irInstance;
 482     init_enxframe(&fr);
 483     gmx::TopologyInformation topInfo;
 484     if (!bDisRe)
 485     {
 486         if (bORIRE || bOTEN)
 487         {
 488             get_orires_parms(ftp2fn(efTPR, NFILE, fnm), ir, &nor, &nex, &or_label, &oobs);
 489         }
 490
 491         if (bORIRE)
 492         {
 493             if (bOrinst)
 494             {
 495                 enx_i = enxORI;
 496             }
 497             else
 498             {
 499                 enx_i = enxOR;
 500             }
 501
 502             if (bORA || bODA)
 503             {
 504                 snew(orient, nor);
 505             }
 506             if (bODR)
 507             {
 508                 snew(odrms, nor);
 509             }
 510             if (bORT || bODT)
 511             {
 512                 fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n");
 513                 fprintf(stderr, "End your selection with 0\n");
 514                 j     = -1;
 515                 orsel = nullptr;
 516                 do
 517                 {
 518                     j++;
 519                     srenew(orsel, j + 1);
 520                     if (1 != scanf("%d", &(orsel[j])))
 521                     {
 522                         gmx_fatal(FARGS, "Error reading user input");
 523                     }
 524                 } while (orsel[j] > 0);
 525                 if (orsel[0] == -1)
 526                 {
 527                     fprintf(stderr, "Selecting all %d orientation restraints\n", nor);
 528                     norsel = nor;
 529                     srenew(orsel, nor);
 530                     for (i = 0; i < nor; i++)
 531                     {
 532                         orsel[i] = i;
 533                     }
 534                 }
 535                 else
 536                 {
 537                     /* Build the selection */
 538                     norsel = 0;
 539                     for (i = 0; i < j; i++)
 540                     {
 541                         for (k = 0; k < nor; k++)
 542                         {
 543                             if (or_label[k] == orsel[i])
 544                             {
 545                                 orsel[norsel] = k;
 546                                 norsel++;
 547                                 break;
 548                             }
 549                         }
 550                         if (k == nor)
 551                         {
 552                             fprintf(stderr, "Orientation restraint label %d not found\n", orsel[i]);
 553                         }
 554                     }
 555                 }
 556                 snew(odtleg, norsel);
 557                 for (i = 0; i < norsel; i++)
 558                 {
 559                     snew(odtleg[i], 256);
 560                     sprintf(odtleg[i], "%d", or_label[orsel[i]]);
 561                 }
 562                 if (bORT)
 563                 {
 564                     fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
 565                                     "Time (ps)", "", oenv);
 566                     if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 567                     {
 568                         fprintf(fort, "%s", orinst_sub);
 569                     }
 570                     xvgr_legend(fort, norsel, odtleg, oenv);
 571                 }
 572                 if (bODT)
 573                 {
 574                     fodt = xvgropen(opt2fn("-odt", NFILE, fnm), "Orientation restraint deviation",
 575                                     "Time (ps)", "", oenv);
 576                     if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 577                     {
 578                         fprintf(fodt, "%s", orinst_sub);
 579                     }
 580                     xvgr_legend(fodt, norsel, odtleg, oenv);
 581                 }
 582                 for (i = 0; i < norsel; i++)
 583                 {
 584                     sfree(odtleg[i]);
 585                 }
 586                 sfree(odtleg);
 587             }
 588         }
 589         if (bOTEN)
 590         {
 591             foten = xvgropen(opt2fn("-oten", NFILE, fnm), "Order tensor", "Time (ps)", "", oenv);
 592             snew(otenleg, bOvec ? nex * 12 : nex * 3);
 593             for (i = 0; i < nex; i++)
 594             {
 595                 for (j = 0; j < 3; j++)
 596                 {
 597                     sprintf(buf, "eig%d", j + 1);
 598                     otenleg[(bOvec ? 12 : 3) * i + j] = gmx_strdup(buf);
 599                 }
 600                 if (bOvec)
 601                 {
 602                     for (j = 0; j < 9; j++)
 603                     {
 604                         sprintf(buf, "vec%d%s", j / 3 + 1, j % 3 == 0 ? "x" : (j % 3 == 1 ? "y" : "z"));
 605                         otenleg[12 * i + 3 + j] = gmx_strdup(buf);
 606                     }
 607                 }
 608             }
 609             xvgr_legend(foten, bOvec ? nex * 12 : nex * 3, otenleg, oenv);
 610             for (j = 0; j < 3; j++)
 611             {
 612                 sfree(otenleg[j]);
 613             }
 614             sfree(otenleg);
 615         }
 616     }
 617     else
 618     {
 619         {
 620             topInfo.fillFromInputFile(ftp2fn(efTPR, NFILE, fnm));
 621             gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO);
 622         }
 623         nbounds = get_bounds(&bounds, &index, &pair, &npairs, &top);
 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(opt2fn("-pairs", NFILE, fnm), "Pair Distances", "Time (ps)",
 630                                 "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                 t_iatom*   fa;
 665                 t_iparams* ip;
 666
 667                 fa = top.idef.il[F_DISRES].iatoms;
 668                 ip = top.idef.iparams;
 669                 if (blk_disre->nsub != 2 || (blk_disre->sub[0].nr != blk_disre->sub[1].nr))
 670                 {
 671                     gmx_incons("Number of disre sub-blocks not equal to 2");
 672                 }
 673
 674                 ndisre = blk_disre->sub[0].nr;
 675                 if (ndisre != top.idef.il[F_DISRES].nr / 3)
 676                 {
 677                     gmx_fatal(FARGS,
 678                               "Number of disre pairs in the energy file (%d) does not match the "
 679                               "number in the run input file (%d)\n",
 680                               ndisre, top.idef.il[F_DISRES].nr / 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], "%d %s %d %s (%d)", resnr_j, anm_j, resnr_k, anm_k,
 692                             ip[fa[3 * i]].disres.label);
 693                 }
 694                 set = select_it(ndisre, pairleg, &nset);
 695                 snew(leg, 2 * nset);
 696                 for (i = 0; (i < nset); i++)
 697                 {
 698                     snew(leg[2 * i], 32);
 699                     sprintf(leg[2 * i], "a %s", pairleg[set[i]]);
 700                     snew(leg[2 * i + 1], 32);
 701                     sprintf(leg[2 * i + 1], "i %s", pairleg[set[i]]);
 702                 }
 703                 xvgr_legend(fp_pairs, 2 * nset, leg, oenv);
 704             }
 705
 706             /*
 707              * Printing time, only when we do not want to skip frames
 708              */
 709             if (!skip || teller % skip == 0)
 710             {
 711                 if (bDisRe)
 712                 {
 713                     /*******************************************
 714                      * D I S T A N C E   R E S T R A I N T S
 715                      *******************************************/
 716                     if (ndisre > 0)
 717                     {
 718                         GMX_RELEASE_ASSERT(blk_disre != nullptr,
 719                                            "Trying to dereference NULL blk_disre pointer");
 720 #if !GMX_DOUBLE
 721                         float* disre_rt     = blk_disre->sub[0].fval;
 722                         float* disre_rm3tav = blk_disre->sub[1].fval;
 723 #else
 724                         double* disre_rt     = blk_disre->sub[0].dval;
 725                         double* disre_rm3tav = blk_disre->sub[1].dval;
 726 #endif
 727
 728                         print_time(out_disre, fr.t);
 729                         if (violaver == nullptr)
 730                         {
 731                             snew(violaver, ndisre);
 732                         }
 733
 734                         /* Subtract bounds from distances, to calculate violations */
 735                         calc_violations(disre_rt, disre_rm3tav, nbounds, pair, bounds, violaver,
 736                                         &sumt, &sumaver);
 737
 738                         fprintf(out_disre, "  %8.4f  %8.4f\n", sumaver, sumt);
 739                         if (bDRAll)
 740                         {
 741                             print_time(fp_pairs, fr.t);
 742                             for (i = 0; (i < nset); i++)
 743                             {
 744                                 sss = set[i];
 745                                 fprintf(fp_pairs, "  %8.4f", mypow(disre_rm3tav[sss], minthird));
 746                                 fprintf(fp_pairs, "  %8.4f", disre_rt[sss]);
 747                             }
 748                             fprintf(fp_pairs, "\n");
 749                         }
 750                         teller_disre++;
 751                     }
 752                 }
 753
 754                 /*******************************************
 755                  * O R I R E S
 756                  *******************************************/
 757                 else
 758                 {
 759                     t_enxblock* blk = find_block_id_enxframe(&fr, enx_i, nullptr);
 760                     if (bORIRE && blk)
 761                     {
 762                         if (blk->nsub != 1)
 763                         {
 764                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly.");
 765                         }
 766
 767                         if (blk->sub[0].nr != nor)
 768                         {
 769                             gmx_fatal(FARGS,
 770                                       "Number of orientation restraints in energy file (%d) does "
 771                                       "not match with the topology (%d)",
 772                                       blk->sub[0].nr, nor);
 773                         }
 774                         if (bORA || bODA)
 775                         {
 776                             for (i = 0; i < nor; i++)
 777                             {
 778                                 orient[i] += blk_value(blk, 0, i);
 779                             }
 780                         }
 781                         if (bODR)
 782                         {
 783                             for (i = 0; i < nor; i++)
 784                             {
 785                                 real v = blk_value(blk, 0, i);
 786                                 odrms[i] += gmx::square(v - oobs[i]);
 787                             }
 788                         }
 789                         if (bORT)
 790                         {
 791                             fprintf(fort, "  %10f", fr.t);
 792                             for (i = 0; i < norsel; i++)
 793                             {
 794                                 fprintf(fort, " %g", blk_value(blk, 0, orsel[i]));
 795                             }
 796                             fprintf(fort, "\n");
 797                         }
 798                         if (bODT)
 799                         {
 800                             fprintf(fodt, "  %10f", fr.t);
 801                             for (i = 0; i < norsel; i++)
 802                             {
 803                                 fprintf(fodt, " %g", blk_value(blk, 0, orsel[i]) - oobs[orsel[i]]);
 804                             }
 805                             fprintf(fodt, "\n");
 806                         }
 807                         norfr++;
 808                     }
 809                     blk = find_block_id_enxframe(&fr, enxORT, nullptr);
 810                     if (bOTEN && blk)
 811                     {
 812                         if (blk->nsub != 1)
 813                         {
 814                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
 815                         }
 816
 817                         if (blk->sub[0].nr != nex * 12)
 818                         {
 819                             gmx_fatal(FARGS,
 820                                       "Number of orientation experiments in energy file (%d) does "
 821                                       "not match with the topology (%d)",
 822                                       blk->sub[0].nr / 12, nex);
 823                         }
 824                         fprintf(foten, "  %10f", fr.t);
 825                         for (i = 0; i < nex; i++)
 826                         {
 827                             for (j = 0; j < (bOvec ? 12 : 3); j++)
 828                             {
 829                                 fprintf(foten, " %g", blk_value(blk, 0, i * 12 + j));
 830                             }
 831                         }
 832                         fprintf(foten, "\n");
 833                     }
 834                 }
 835             }
 836             teller++;
 837         }
 838     } while (bCont && (timecheck == 0));
 839     free_enxframe(&fr);
 840
 841     fprintf(stderr, "\n");
 842     done_ener_file(fp);
 843     if (out_disre)
 844     {
 845         xvgrclose(out_disre);
 846     }
 847
 848     if (bDRAll)
 849     {
 850         xvgrclose(fp_pairs);
 851     }
 852
 853     if (bORT)
 854     {
 855         xvgrclose(fort);
 856     }
 857     if (bODT)
 858     {
 859         xvgrclose(fodt);
 860     }
 861     if (bORA)
 862     {
 863         FILE* out = xvgropen(opt2fn("-ora", NFILE, fnm), "Average calculated orientations",
 864                              "Restraint label", "", oenv);
 865         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 866         {
 867             fprintf(out, "%s", orinst_sub);
 868         }
 869         for (i = 0; i < nor; i++)
 870         {
 871             fprintf(out, "%5d  %g\n", or_label[i], orient[i] / real(norfr));
 872         }
 873         xvgrclose(out);
 874     }
 875     if (bODA)
 876     {
 877         FILE* out = xvgropen(opt2fn("-oda", NFILE, fnm), "Average restraint deviation",
 878                              "Restraint label", "", oenv);
 879         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 880         {
 881             fprintf(out, "%s", orinst_sub);
 882         }
 883         for (i = 0; i < nor; i++)
 884         {
 885             fprintf(out, "%5d  %g\n", or_label[i], orient[i] / real(norfr) - oobs[i]);
 886         }
 887         xvgrclose(out);
 888     }
 889     if (bODR)
 890     {
 891         FILE* out = xvgropen(opt2fn("-odr", NFILE, fnm), "RMS orientation restraint deviations",
 892                              "Restraint label", "", oenv);
 893         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
 894         {
 895             fprintf(out, "%s", orinst_sub);
 896         }
 897         for (i = 0; i < nor; i++)
 898         {
 899             fprintf(out, "%5d  %g\n", or_label[i], std::sqrt(odrms[i] / real(norfr)));
 900         }
 901         xvgrclose(out);
 902     }
 903     // Clean up orires variables.
 904     sfree(or_label);
 905     sfree(oobs);
 906     sfree(orient);
 907     sfree(odrms);
 908     sfree(orsel);
 909     if (bOTEN)
 910     {
 911         xvgrclose(foten);
 912     }
 913
 914     if (bDisRe)
 915     {
 916         analyse_disre(opt2fn("-viol", NFILE, fnm), teller_disre, violaver, bounds, index, pair,
 917                       nbounds, oenv);
 918     }
 919     {
 920         const char* nxy = "-nxy";
 921
 922         do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
 923         do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
 924         do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
 925         do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
 926         do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
 927         do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
 928     }
 929     output_env_done(oenv);
 930
 931     return 0;
 932 }