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