src/gromacs/gmxana/gmx_nmr.cpp

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