src/gromacs/gmxana/gmx_nmr.cpp

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