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