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