src/gromacs/gmxana/gmx_energy.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/utility/arraysize.h"
  68 #include "gromacs/utility/cstringutil.h"
  69 #include "gromacs/utility/fatalerror.h"
  70 #include "gromacs/utility/gmxassert.h"
  71 #include "gromacs/utility/pleasecite.h"
  72 #include "gromacs/utility/smalloc.h"
  73 #include "gromacs/utility/strconvert.h"
  74
  75 static const int NOTSET = -23451;
  76
  77 typedef struct
  78 {
  79     real sum;
  80     real sum2;
  81 } exactsum_t;
  82
  83 typedef struct
  84 {
  85     real*       ener;
  86     exactsum_t* es;
  87     gmx_bool    bExactStat;
  88     double      av;
  89     double      rmsd;
  90     double      ee;
  91     double      slope;
  92 } enerdat_t;
  93
  94 typedef struct
  95 {
  96     int64_t    nsteps;
  97     int64_t    npoints;
  98     int        nframes;
  99     int*       step;
 100     int*       steps;
 101     int*       points;
 102     enerdat_t* s;
 103     gmx_bool   bHaveSums;
 104 } enerdata_t;
 105
 106 static void done_enerdata_t(int nset, enerdata_t* edat)
 107 {
 108     sfree(edat->step);
 109     sfree(edat->steps);
 110     sfree(edat->points);
 111     for (int i = 0; i < nset; i++)
 112     {
 113         sfree(edat->s[i].ener);
 114         sfree(edat->s[i].es);
 115     }
 116     sfree(edat->s);
 117 }
 118
 119 static void chomp(char* buf)
 120 {
 121     int len = std::strlen(buf);
 122
 123     while ((len > 0) && (buf[len - 1] == '\n'))
 124     {
 125         buf[len - 1] = '\0';
 126         len--;
 127     }
 128 }
 129
 130 static int* select_by_name(int nre, gmx_enxnm_t* nm, int* nset)
 131 {
 132     gmx_bool*   bE;
 133     int         k, kk, j, i, nmatch, nind, nss;
 134     int*        set;
 135     gmx_bool    bEOF, bVerbose = TRUE, bLong = FALSE;
 136     char *      ptr, buf[STRLEN];
 137     const char* fm4   = "%3d  %-14s";
 138     const char* fm2   = "%3d  %-34s";
 139     char**      newnm = nullptr;
 140
 141     if ((getenv("GMX_ENER_VERBOSE")) != nullptr)
 142     {
 143         bVerbose = FALSE;
 144     }
 145
 146     fprintf(stderr, "\n");
 147     fprintf(stderr, "Select the terms you want from the following list by\n");
 148     fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
 149     fprintf(stderr, "End your selection with an empty line or a zero.\n");
 150     fprintf(stderr, "-------------------------------------------------------------------\n");
 151
 152     snew(newnm, nre);
 153     j = 0;
 154     for (k = 0; k < nre; k++)
 155     {
 156         newnm[k] = gmx_strdup(nm[k].name);
 157         /* Insert dashes in all the names */
 158         while ((ptr = std::strchr(newnm[k], ' ')) != nullptr)
 159         {
 160             *ptr = '-';
 161         }
 162         if (bVerbose)
 163         {
 164             if (j == 0)
 165             {
 166                 if (k > 0)
 167                 {
 168                     fprintf(stderr, "\n");
 169                 }
 170                 bLong = FALSE;
 171                 for (kk = k; kk < k + 4; kk++)
 172                 {
 173                     if (kk < nre && std::strlen(nm[kk].name) > 14)
 174                     {
 175                         bLong = TRUE;
 176                     }
 177                 }
 178             }
 179             else
 180             {
 181                 fprintf(stderr, " ");
 182             }
 183             if (!bLong)
 184             {
 185                 fprintf(stderr, fm4, k + 1, newnm[k]);
 186                 j++;
 187                 if (j == 4)
 188                 {
 189                     j = 0;
 190                 }
 191             }
 192             else
 193             {
 194                 fprintf(stderr, fm2, k + 1, newnm[k]);
 195                 j++;
 196                 if (j == 2)
 197                 {
 198                     j = 0;
 199                 }
 200             }
 201         }
 202     }
 203     if (bVerbose)
 204     {
 205         fprintf(stderr, "\n\n");
 206     }
 207
 208     snew(bE, nre);
 209
 210     bEOF = FALSE;
 211     while (!bEOF && (fgets2(buf, STRLEN - 1, stdin)))
 212     {
 213         /* Remove newlines */
 214         chomp(buf);
 215
 216         /* Remove spaces */
 217         trim(buf);
 218
 219         /* Empty line means end of input */
 220         bEOF = (std::strlen(buf) == 0);
 221         if (!bEOF)
 222         {
 223             ptr = buf;
 224             do
 225             {
 226                 if (!bEOF)
 227                 {
 228                     /* First try to read an integer */
 229                     nss = sscanf(ptr, "%d", &nind);
 230                     if (nss == 1)
 231                     {
 232                         /* Zero means end of input */
 233                         if (nind == 0)
 234                         {
 235                             bEOF = TRUE;
 236                         }
 237                         else if ((1 <= nind) && (nind <= nre))
 238                         {
 239                             bE[nind - 1] = TRUE;
 240                         }
 241                         else
 242                         {
 243                             fprintf(stderr, "number %d is out of range\n", nind);
 244                         }
 245                     }
 246                     else
 247                     {
 248                         /* Now try to read a string */
 249                         nmatch = 0;
 250                         for (nind = 0; nind < nre; nind++)
 251                         {
 252                             if (gmx_strcasecmp(newnm[nind], ptr) == 0)
 253                             {
 254                                 bE[nind] = TRUE;
 255                                 nmatch++;
 256                             }
 257                         }
 258                         if (nmatch == 0)
 259                         {
 260                             i      = std::strlen(ptr);
 261                             nmatch = 0;
 262                             for (nind = 0; nind < nre; nind++)
 263                             {
 264                                 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
 265                                 {
 266                                     bE[nind] = TRUE;
 267                                     nmatch++;
 268                                 }
 269                             }
 270                             if (nmatch == 0)
 271                             {
 272                                 fprintf(stderr, "String '%s' does not match anything\n", ptr);
 273                             }
 274                         }
 275                     }
 276                 }
 277                 /* Look for the first space, and remove spaces from there */
 278                 if ((ptr = std::strchr(ptr, ' ')) != nullptr)
 279                 {
 280                     trim(ptr);
 281                 }
 282             } while (!bEOF && ((ptr != nullptr) && (std::strlen(ptr) > 0)));
 283         }
 284     }
 285
 286     snew(set, nre);
 287     for (i = (*nset) = 0; (i < nre); i++)
 288     {
 289         if (bE[i])
 290         {
 291             set[(*nset)++] = i;
 292         }
 293     }
 294
 295     sfree(bE);
 296
 297     if (*nset == 0)
 298     {
 299         gmx_fatal(FARGS, "No energy terms selected");
 300     }
 301
 302     for (i = 0; (i < nre); i++)
 303     {
 304         sfree(newnm[i]);
 305     }
 306     sfree(newnm);
 307
 308     return set;
 309 }
 310
 311 static void get_dhdl_parms(const char* topnm, t_inputrec* ir)
 312 {
 313     gmx_mtop_t mtop;
 314     int        natoms;
 315     matrix     box;
 316
 317     /* all we need is the ir to be able to write the label */
 318     read_tpx(topnm, ir, box, &natoms, nullptr, nullptr, &mtop);
 319 }
 320
 321 static void einstein_visco(const char*             fn,
 322                            const char*             fni,
 323                            int                     nsets,
 324                            int                     nint,
 325                            real**                  eneint,
 326                            real                    V,
 327                            real                    T,
 328                            double                  dt,
 329                            const gmx_output_env_t* oenv)
 330 {
 331     FILE * fp0, *fp1;
 332     double av[4], avold[4];
 333     double fac, di;
 334     int    i, j, m, nf4;
 335
 336     nf4 = nint / 4 + 1;
 337
 338     for (i = 0; i <= nsets; i++)
 339     {
 340         avold[i] = 0;
 341     }
 342     fp0 = xvgropen(fni, "Shear viscosity integral", "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
 343     fp1 = xvgropen(
 344             fn, "Shear viscosity using Einstein relation", "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
 345     for (i = 0; i < nf4; i++)
 346     {
 347         for (m = 0; m <= nsets; m++)
 348         {
 349             av[m] = 0;
 350         }
 351         for (j = 0; j < nint - i; j++)
 352         {
 353             for (m = 0; m < nsets; m++)
 354             {
 355                 di = gmx::square(eneint[m][j + i] - eneint[m][j]);
 356
 357                 av[m] += di;
 358                 av[nsets] += di / nsets;
 359             }
 360         }
 361         /* Convert to SI for the viscosity */
 362         fac = (V * gmx::c_nano * gmx::c_nano * gmx::c_nano * gmx::c_pico * 1e10)
 363               / (2 * gmx::c_boltzmann * T) / (nint - i);
 364         fprintf(fp0, "%10g", i * dt);
 365         for (m = 0; (m <= nsets); m++)
 366         {
 367             av[m] = fac * av[m];
 368             fprintf(fp0, "  %10g", av[m]);
 369         }
 370         fprintf(fp0, "\n");
 371         fprintf(fp1, "%10g", (i + 0.5) * dt);
 372         for (m = 0; (m <= nsets); m++)
 373         {
 374             fprintf(fp1, "  %10g", (av[m] - avold[m]) / dt);
 375             avold[m] = av[m];
 376         }
 377         fprintf(fp1, "\n");
 378     }
 379     xvgrclose(fp0);
 380     xvgrclose(fp1);
 381 }
 382
 383 typedef struct
 384 {
 385     int64_t np;
 386     double  sum;
 387     double  sav;
 388     double  sav2;
 389 } ee_sum_t;
 390
 391 typedef struct
 392 {
 393     int      b;
 394     ee_sum_t sum;
 395     int64_t  nst;
 396     int64_t  nst_min;
 397 } ener_ee_t;
 398
 399 static void clear_ee_sum(ee_sum_t* ees)
 400 {
 401     ees->sav  = 0;
 402     ees->sav2 = 0;
 403     ees->np   = 0;
 404     ees->sum  = 0;
 405 }
 406
 407 static void add_ee_sum(ee_sum_t* ees, double sum, int np)
 408 {
 409     ees->np += np;
 410     ees->sum += sum;
 411 }
 412
 413 static void add_ee_av(ee_sum_t* ees)
 414 {
 415     double av;
 416
 417     av = ees->sum / ees->np;
 418     ees->sav += av;
 419     ees->sav2 += av * av;
 420     ees->np  = 0;
 421     ees->sum = 0;
 422 }
 423
 424 static double calc_ee2(int nb, ee_sum_t* ees)
 425 {
 426     return (ees->sav2 / nb - gmx::square(ees->sav / nb)) / (nb - 1);
 427 }
 428
 429 static void set_ee_av(ener_ee_t* eee)
 430 {
 431     if (debug)
 432     {
 433         char buf[STEPSTRSIZE];
 434         fprintf(debug, "Storing average for err.est.: %s steps\n", gmx_step_str(eee->nst, buf));
 435     }
 436     add_ee_av(&eee->sum);
 437     eee->b++;
 438     if (eee->b == 1 || eee->nst < eee->nst_min)
 439     {
 440         eee->nst_min = eee->nst;
 441     }
 442     eee->nst = 0;
 443 }
 444
 445 static void calc_averages(int nset, enerdata_t* edat, int nbmin, int nbmax)
 446 {
 447     int         nb, i, f, nee;
 448     double      sum, sum2, sump, see2;
 449     int64_t     np, p, bound_nb;
 450     enerdat_t*  ed;
 451     exactsum_t* es;
 452     gmx_bool    bAllZero;
 453     double      x, sx, sy, sxx, sxy;
 454     ener_ee_t*  eee;
 455
 456     /* Check if we have exact statistics over all points */
 457     for (i = 0; i < nset; i++)
 458     {
 459         ed             = &edat->s[i];
 460         ed->bExactStat = FALSE;
 461         if (edat->bHaveSums)
 462         {
 463             /* All energy file sum entries 0 signals no exact sums.
 464              * But if all energy values are 0, we still have exact sums.
 465              */
 466             bAllZero = TRUE;
 467             for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
 468             {
 469                 if (ed->ener[i] != 0)
 470                 {
 471                     bAllZero = FALSE;
 472                 }
 473                 ed->bExactStat = (ed->es[f].sum != 0);
 474             }
 475             if (bAllZero)
 476             {
 477                 ed->bExactStat = TRUE;
 478             }
 479         }
 480     }
 481
 482     snew(eee, nbmax + 1);
 483     for (i = 0; i < nset; i++)
 484     {
 485         ed = &edat->s[i];
 486
 487         sum  = 0;
 488         sum2 = 0;
 489         np   = 0;
 490         sx   = 0;
 491         sy   = 0;
 492         sxx  = 0;
 493         sxy  = 0;
 494         for (nb = nbmin; nb <= nbmax; nb++)
 495         {
 496             eee[nb].b = 0;
 497             clear_ee_sum(&eee[nb].sum);
 498             eee[nb].nst     = 0;
 499             eee[nb].nst_min = 0;
 500         }
 501         for (f = 0; f < edat->nframes; f++)
 502         {
 503             es = &ed->es[f];
 504
 505             if (ed->bExactStat)
 506             {
 507                 /* Add the sum and the sum of variances to the totals. */
 508                 p    = edat->points[f];
 509                 sump = es->sum;
 510                 sum2 += es->sum2;
 511                 if (np > 0)
 512                 {
 513                     sum2 += gmx::square(sum / np - (sum + es->sum) / (np + p)) * np * (np + p) / p;
 514                 }
 515             }
 516             else
 517             {
 518                 /* Add a single value to the sum and sum of squares. */
 519                 p    = 1;
 520                 sump = ed->ener[f];
 521                 sum2 += gmx::square(sump);
 522             }
 523
 524             /* sum has to be increased after sum2 */
 525             np += p;
 526             sum += sump;
 527
 528             /* For the linear regression use variance 1/p.
 529              * Note that sump is the sum, not the average, so we don't need p*.
 530              */
 531             x = edat->step[f] - 0.5 * (edat->steps[f] - 1);
 532             sx += p * x;
 533             sy += sump;
 534             sxx += p * x * x;
 535             sxy += x * sump;
 536
 537             for (nb = nbmin; nb <= nbmax; nb++)
 538             {
 539                 /* Check if the current end step is closer to the desired
 540                  * block boundary than the next end step.
 541                  */
 542                 bound_nb = (edat->step[0] - 1) * nb + edat->nsteps * (eee[nb].b + 1);
 543                 if (eee[nb].nst > 0 && bound_nb - edat->step[f - 1] * nb < edat->step[f] * nb - bound_nb)
 544                 {
 545                     set_ee_av(&eee[nb]);
 546                 }
 547                 if (f == 0)
 548                 {
 549                     eee[nb].nst = 1;
 550                 }
 551                 else
 552                 {
 553                     eee[nb].nst += edat->step[f] - edat->step[f - 1];
 554                 }
 555                 if (ed->bExactStat)
 556                 {
 557                     add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
 558                 }
 559                 else
 560                 {
 561                     add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
 562                 }
 563                 bound_nb = (edat->step[0] - 1) * nb + edat->nsteps * (eee[nb].b + 1);
 564                 if (edat->step[f] * nb >= bound_nb)
 565                 {
 566                     set_ee_av(&eee[nb]);
 567                 }
 568             }
 569         }
 570
 571         edat->s[i].av = sum / np;
 572         if (ed->bExactStat)
 573         {
 574             edat->s[i].rmsd = std::sqrt(sum2 / np);
 575         }
 576         else
 577         {
 578             edat->s[i].rmsd = std::sqrt(sum2 / np - gmx::square(edat->s[i].av));
 579         }
 580
 581         if (edat->nframes > 1)
 582         {
 583             edat->s[i].slope = (np * sxy - sx * sy) / (np * sxx - sx * sx);
 584         }
 585         else
 586         {
 587             edat->s[i].slope = 0;
 588         }
 589
 590         nee  = 0;
 591         see2 = 0;
 592         for (nb = nbmin; nb <= nbmax; nb++)
 593         {
 594             /* Check if we actually got nb blocks and if the smallest
 595              * block is not shorter than 80% of the average.
 596              */
 597             if (debug)
 598             {
 599                 char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
 600                 fprintf(debug,
 601                         "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
 602                         nb,
 603                         eee[nb].b,
 604                         gmx_step_str(eee[nb].nst_min, buf1),
 605                         gmx_step_str(edat->nsteps, buf2));
 606             }
 607             if (eee[nb].b == nb && 5 * nb * eee[nb].nst_min >= 4 * edat->nsteps)
 608             {
 609                 see2 += calc_ee2(nb, &eee[nb].sum);
 610                 nee++;
 611             }
 612         }
 613         if (nee > 0)
 614         {
 615             edat->s[i].ee = std::sqrt(see2 / nee);
 616         }
 617         else
 618         {
 619             edat->s[i].ee = -1;
 620         }
 621     }
 622     sfree(eee);
 623 }
 624
 625 static enerdata_t* calc_sum(int nset, enerdata_t* edat, int nbmin, int nbmax)
 626 {
 627     enerdata_t* esum;
 628     enerdat_t*  s;
 629     int         f, i;
 630     double      sum;
 631
 632     snew(esum, 1);
 633     *esum = *edat;
 634     snew(esum->s, 1);
 635     s = &esum->s[0];
 636     snew(s->ener, esum->nframes);
 637     snew(s->es, esum->nframes);
 638
 639     s->bExactStat = TRUE;
 640     s->slope      = 0;
 641     for (i = 0; i < nset; i++)
 642     {
 643         if (!edat->s[i].bExactStat)
 644         {
 645             s->bExactStat = FALSE;
 646         }
 647         s->slope += edat->s[i].slope;
 648     }
 649
 650     for (f = 0; f < edat->nframes; f++)
 651     {
 652         sum = 0;
 653         for (i = 0; i < nset; i++)
 654         {
 655             sum += edat->s[i].ener[f];
 656         }
 657         s->ener[f] = sum;
 658         sum        = 0;
 659         for (i = 0; i < nset; i++)
 660         {
 661             sum += edat->s[i].es[f].sum;
 662         }
 663         s->es[f].sum  = sum;
 664         s->es[f].sum2 = 0;
 665     }
 666
 667     calc_averages(1, esum, nbmin, nbmax);
 668
 669     return esum;
 670 }
 671
 672 static void ee_pr(double ee, int buflen, char* buf)
 673 {
 674     snprintf(buf, buflen, "%s", "--");
 675     if (ee >= 0)
 676     {
 677         /* Round to two decimals by printing. */
 678         char tmp[100];
 679         snprintf(tmp, sizeof(tmp), "%.1e", ee);
 680         double rnd = gmx::doubleFromString(tmp);
 681         snprintf(buf, buflen, "%g", rnd);
 682     }
 683 }
 684
 685 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t* edat)
 686 {
 687     /* Remove the drift by performing a fit to y = ax+b.
 688        Uses 5 iterations. */
 689     int    i, j, k;
 690     double delta;
 691
 692     edat->npoints = edat->nframes;
 693     edat->nsteps  = edat->nframes;
 694
 695     for (k = 0; (k < 5); k++)
 696     {
 697         for (i = 0; (i < nset); i++)
 698         {
 699             delta = edat->s[i].slope * dt;
 700
 701             if (nullptr != debug)
 702             {
 703                 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
 704             }
 705
 706             for (j = 0; (j < edat->nframes); j++)
 707             {
 708                 edat->s[i].ener[j] -= j * delta;
 709                 edat->s[i].es[j].sum  = 0;
 710                 edat->s[i].es[j].sum2 = 0;
 711             }
 712         }
 713         calc_averages(nset, edat, nbmin, nbmax);
 714     }
 715 }
 716
 717 static void calc_fluctuation_props(FILE*       fp,
 718                                    gmx_bool    bDriftCorr,
 719                                    real        dt,
 720                                    int         nset,
 721                                    int         nmol,
 722                                    char**      leg,
 723                                    enerdata_t* edat,
 724                                    int         nbmin,
 725                                    int         nbmax)
 726 {
 727     int    i, j;
 728     double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, varet;
 729     double NANO3;
 730     enum
 731     {
 732         eVol,
 733         eEnth,
 734         eTemp,
 735         eEtot,
 736         eNR
 737     };
 738     const char* my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
 739     int         ii[eNR];
 740
 741     NANO3 = gmx::c_nano * gmx::c_nano * gmx::c_nano;
 742     if (!bDriftCorr)
 743     {
 744         fprintf(fp,
 745                 "\nYou may want to use the -driftcorr flag in order to correct\n"
 746                 "for spurious drift in the graphs. Note that this is not\n"
 747                 "a substitute for proper equilibration and sampling!\n");
 748     }
 749     else
 750     {
 751         remove_drift(nset, nbmin, nbmax, dt, edat);
 752     }
 753     for (i = 0; (i < eNR); i++)
 754     {
 755         for (ii[i] = 0; (ii[i] < nset && (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++) {}
 756 /*        if (ii[i] < nset)
 757             fprintf(fp,"Found %s data.\n",my_ener[i]);
 758  */ }
 759 /* Compute it all! */
 760 alpha = kappa = cp = dcp = cv = NOTSET;
 761
 762 /* Temperature */
 763 tt = NOTSET;
 764 if (ii[eTemp] < nset)
 765 {
 766     tt = edat->s[ii[eTemp]].av;
 767 }
 768 /* Volume */
 769 vv = varv = NOTSET;
 770 if ((ii[eVol] < nset) && (ii[eTemp] < nset))
 771 {
 772     vv    = edat->s[ii[eVol]].av * NANO3;
 773     varv  = gmx::square(edat->s[ii[eVol]].rmsd * NANO3);
 774     kappa = (varv / vv) / (gmx::c_boltzmann * tt);
 775 }
 776 /* Enthalpy */
 777 hh = varh = NOTSET;
 778 if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
 779 {
 780     hh   = gmx::c_kilo * edat->s[ii[eEnth]].av / gmx::c_avogadro;
 781     varh = gmx::square(gmx::c_kilo * edat->s[ii[eEnth]].rmsd / gmx::c_avogadro);
 782     cp   = gmx::c_avogadro * ((varh / nmol) / (gmx::c_boltzmann * tt * tt));
 783 }
 784 /* Total energy */
 785 if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
 786 {
 787     /* Only compute cv in constant volume runs, which we can test
 788        by checking whether the enthalpy was computed.
 789      */
 790     varet = gmx::square(edat->s[ii[eEtot]].rmsd);
 791     cv    = gmx::c_kilo * ((varet / nmol) / (gmx::c_boltz * tt * tt));
 792 }
 793 /* Alpha, dcp */
 794 if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
 795 {
 796     double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
 797     vh_sum = v_sum = h_sum = 0;
 798     for (j = 0; (j < edat->nframes); j++)
 799     {
 800         v = edat->s[ii[eVol]].ener[j] * NANO3;
 801         h = gmx::c_kilo * edat->s[ii[eEnth]].ener[j] / gmx::c_avogadro;
 802         v_sum += v;
 803         h_sum += h;
 804         vh_sum += (v * h);
 805     }
 806     vh_aver = vh_sum / edat->nframes;
 807     v_aver  = v_sum / edat->nframes;
 808     h_aver  = h_sum / edat->nframes;
 809     alpha   = (vh_aver - v_aver * h_aver) / (v_aver * gmx::c_boltzmann * tt * tt);
 810     dcp     = (v_aver * gmx::c_avogadro / nmol) * tt * gmx::square(alpha) / (kappa);
 811 }
 812
 813 if (tt != NOTSET)
 814 {
 815     if (nmol < 2)
 816     {
 817         fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n", nmol);
 818     }
 819     fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
 820     fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
 821     fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
 822     fprintf(fp, "please use the g_dos program.\n\n");
 823     fprintf(fp,
 824             "WARNING: Please verify that your simulations are converged and perform\n"
 825             "a block-averaging error analysis (not implemented in g_energy yet)\n");
 826
 827     if (debug != nullptr)
 828     {
 829         if (varv != NOTSET)
 830         {
 831             fprintf(fp, "varv  =  %10g (m^6)\n", varv * gmx::c_avogadro / nmol);
 832         }
 833     }
 834     if (vv != NOTSET)
 835     {
 836         fprintf(fp, "Volume                                   = %10g m^3/mol\n", vv * gmx::c_avogadro / nmol);
 837     }
 838     if (varh != NOTSET)
 839     {
 840         fprintf(fp,
 841                 "Enthalpy                                 = %10g kJ/mol\n",
 842                 hh * gmx::c_avogadro / (gmx::c_kilo * nmol));
 843     }
 844     if (alpha != NOTSET)
 845     {
 846         fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n", alpha);
 847     }
 848     if (kappa != NOTSET)
 849     {
 850         fprintf(fp, "Isothermal Compressibility Kappa         = %10g (m^3/J)\n", kappa);
 851         fprintf(fp, "Adiabatic bulk modulus                   = %10g (J/m^3)\n", 1.0 / kappa);
 852     }
 853     if (cp != NOTSET)
 854     {
 855         fprintf(fp, "Heat capacity at constant pressure Cp    = %10g J/(mol K)\n", cp);
 856     }
 857     if (cv != NOTSET)
 858     {
 859         fprintf(fp, "Heat capacity at constant volume Cv      = %10g J/(mol K)\n", cv);
 860     }
 861     if (dcp != NOTSET)
 862     {
 863         fprintf(fp, "Cp-Cv                                    =  %10g J/(mol K)\n", dcp);
 864     }
 865     please_cite(fp, "Allen1987a");
 866 }
 867 else
 868 {
 869     fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
 870 }
 871 }
 872
 873 static void analyse_ener(gmx_bool                bCorr,
 874                          const char*             corrfn,
 875                          const char*             eviscofn,
 876                          const char*             eviscoifn,
 877                          gmx_bool                bFee,
 878                          gmx_bool                bSum,
 879                          gmx_bool                bFluct,
 880                          gmx_bool                bVisco,
 881                          const char*             visfn,
 882                          int                     nmol,
 883                          int64_t                 start_step,
 884                          double                  start_t,
 885                          int64_t                 step,
 886                          double                  t,
 887                          real                    reftemp,
 888                          enerdata_t*             edat,
 889                          int                     nset,
 890                          const int               set[],
 891                          const gmx_bool*         bIsEner,
 892                          char**                  leg,
 893                          gmx_enxnm_t*            enm,
 894                          real                    Vaver,
 895                          real                    ezero,
 896                          int                     nbmin,
 897                          int                     nbmax,
 898                          const gmx_output_env_t* oenv)
 899 {
 900     FILE* fp;
 901     /* Check out the printed manual for equations! */
 902     double      Dt, aver, stddev, errest, delta_t, totaldrift;
 903     enerdata_t* esum = nullptr;
 904     real        integral, intBulk, Temp = 0, Pres = 0;
 905     real        pr_aver, pr_stddev, pr_errest;
 906     double      beta = 0, expE, expEtot, *fee = nullptr;
 907     int64_t     nsteps;
 908     int         nexact, nnotexact;
 909     int         i, j, nout;
 910     char        buf[256], eebuf[100];
 911
 912     nsteps = step - start_step + 1;
 913     if (nsteps < 1)
 914     {
 915         fprintf(stdout, "Not enough steps (%s) for statistics\n", gmx_step_str(nsteps, buf));
 916     }
 917     else
 918     {
 919         /* Calculate the time difference */
 920         delta_t = t - start_t;
 921
 922         fprintf(stdout,
 923                 "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
 924                 gmx_step_str(nsteps, buf),
 925                 start_t,
 926                 t,
 927                 nset);
 928
 929         calc_averages(nset, edat, nbmin, nbmax);
 930
 931         if (bSum)
 932         {
 933             esum = calc_sum(nset, edat, nbmin, nbmax);
 934         }
 935
 936         if (!edat->bHaveSums)
 937         {
 938             nexact    = 0;
 939             nnotexact = nset;
 940         }
 941         else
 942         {
 943             nexact    = 0;
 944             nnotexact = 0;
 945             for (i = 0; (i < nset); i++)
 946             {
 947                 if (edat->s[i].bExactStat)
 948                 {
 949                     nexact++;
 950                 }
 951                 else
 952                 {
 953                     nnotexact++;
 954                 }
 955             }
 956         }
 957
 958         if (nnotexact == 0)
 959         {
 960             fprintf(stdout, "All statistics are over %s points\n", gmx_step_str(edat->npoints, buf));
 961         }
 962         else if (nexact == 0 || edat->npoints == edat->nframes)
 963         {
 964             fprintf(stdout, "All statistics are over %d points (frames)\n", edat->nframes);
 965         }
 966         else
 967         {
 968             fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
 969             for (i = 0; (i < nset); i++)
 970             {
 971                 if (!edat->s[i].bExactStat)
 972                 {
 973                     fprintf(stdout, " '%s'", leg[i]);
 974                 }
 975             }
 976             fprintf(stdout,
 977                     " %s has statistics over %d points (frames)\n",
 978                     nnotexact == 1 ? "is" : "are",
 979                     edat->nframes);
 980             fprintf(stdout, "All other statistics are over %s points\n", gmx_step_str(edat->npoints, buf));
 981         }
 982         fprintf(stdout, "\n");
 983
 984         fprintf(stdout,
 985                 "%-24s %10s %10s %10s %10s",
 986                 "Energy",
 987                 "Average",
 988                 "Err.Est.",
 989                 "RMSD",
 990                 "Tot-Drift");
 991         if (bFee)
 992         {
 993             fprintf(stdout, "  %10s\n", "-kT ln<e^(E/kT)>");
 994         }
 995         else
 996         {
 997             fprintf(stdout, "\n");
 998         }
 999         fprintf(stdout,
1000                 "-------------------------------------------------------------------------------"
1001                 "\n");
1002
1003         /* Initiate locals, only used with -sum */
1004         expEtot = 0;
1005         if (bFee)
1006         {
1007             beta = 1.0 / (gmx::c_boltz * reftemp);
1008             snew(fee, nset);
1009         }
1010         for (i = 0; (i < nset); i++)
1011         {
1012             aver   = edat->s[i].av;
1013             stddev = edat->s[i].rmsd;
1014             errest = edat->s[i].ee;
1015
1016             if (bFee)
1017             {
1018                 expE = 0;
1019                 for (j = 0; (j < edat->nframes); j++)
1020                 {
1021                     expE += std::exp(beta * (edat->s[i].ener[j] - aver) / nmol);
1022                 }
1023                 if (bSum)
1024                 {
1025                     expEtot += expE / edat->nframes;
1026                 }
1027
1028                 fee[i] = std::log(expE / edat->nframes) / beta + aver / nmol;
1029             }
1030             if (std::strstr(leg[i], "empera") != nullptr)
1031             {
1032                 Temp = aver;
1033             }
1034             else if (std::strstr(leg[i], "olum") != nullptr)
1035             {
1036                 Vaver = aver;
1037             }
1038             else if (std::strstr(leg[i], "essure") != nullptr)
1039             {
1040                 Pres = aver;
1041             }
1042             if (bIsEner[i])
1043             {
1044                 pr_aver   = aver / nmol - ezero;
1045                 pr_stddev = stddev / nmol;
1046                 pr_errest = errest / nmol;
1047             }
1048             else
1049             {
1050                 pr_aver   = aver;
1051                 pr_stddev = stddev;
1052                 pr_errest = errest;
1053             }
1054
1055             /* Multiply the slope in steps with the number of steps taken */
1056             totaldrift = (edat->nsteps - 1) * edat->s[i].slope;
1057             if (bIsEner[i])
1058             {
1059                 totaldrift /= nmol;
1060             }
1061
1062             ee_pr(pr_errest, sizeof(eebuf), eebuf);
1063             fprintf(stdout, "%-24s %10g %10s %10g %10g", leg[i], pr_aver, eebuf, pr_stddev, totaldrift);
1064             if (bFee)
1065             {
1066                 fprintf(stdout, "  %10g", fee[i]);
1067             }
1068
1069             fprintf(stdout, "  (%s)\n", enm[set[i]].unit);
1070
1071             if (bFluct)
1072             {
1073                 for (j = 0; (j < edat->nframes); j++)
1074                 {
1075                     edat->s[i].ener[j] -= aver;
1076                 }
1077             }
1078         }
1079         if (bSum)
1080         {
1081             totaldrift = (edat->nsteps - 1) * esum->s[0].slope;
1082             ee_pr(esum->s[0].ee / nmol, sizeof(eebuf), eebuf);
1083             fprintf(stdout,
1084                     "%-24s %10g %10s %10s %10g  (%s)",
1085                     "Total",
1086                     esum->s[0].av / nmol,
1087                     eebuf,
1088                     "--",
1089                     totaldrift / nmol,
1090                     enm[set[0]].unit);
1091             /* pr_aver,pr_stddev,a,totaldrift */
1092             if (bFee)
1093             {
1094                 fprintf(stdout,
1095                         "  %10g  %10g\n",
1096                         std::log(expEtot) / beta + esum->s[0].av / nmol,
1097                         std::log(expEtot) / beta);
1098             }
1099             else
1100             {
1101                 fprintf(stdout, "\n");
1102             }
1103         }
1104
1105         /* Do correlation function */
1106         if (edat->nframes > 1)
1107         {
1108             Dt = delta_t / (edat->nframes - 1);
1109         }
1110         else
1111         {
1112             Dt = 0;
1113         }
1114         if (bVisco)
1115         {
1116             const char* leg[] = { "Shear", "Bulk" };
1117             real        factor;
1118             real**      eneset;
1119             real**      eneint;
1120
1121             /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1122
1123             /* Symmetrise tensor! (and store in first three elements)
1124              * And subtract average pressure!
1125              */
1126             snew(eneset, 12);
1127             for (i = 0; i < 12; i++)
1128             {
1129                 snew(eneset[i], edat->nframes);
1130             }
1131             for (i = 0; (i < edat->nframes); i++)
1132             {
1133                 eneset[0][i] = 0.5 * (edat->s[1].ener[i] + edat->s[3].ener[i]);
1134                 eneset[1][i] = 0.5 * (edat->s[2].ener[i] + edat->s[6].ener[i]);
1135                 eneset[2][i] = 0.5 * (edat->s[5].ener[i] + edat->s[7].ener[i]);
1136                 for (j = 3; j <= 11; j++)
1137                 {
1138                     eneset[j][i] = edat->s[j].ener[i];
1139                 }
1140                 eneset[11][i] -= Pres;
1141             }
1142
1143             /* Determine integrals of the off-diagonal pressure elements */
1144             snew(eneint, 3);
1145             for (i = 0; i < 3; i++)
1146             {
1147                 snew(eneint[i], edat->nframes + 1);
1148             }
1149             eneint[0][0] = 0;
1150             eneint[1][0] = 0;
1151             eneint[2][0] = 0;
1152             for (i = 0; i < edat->nframes; i++)
1153             {
1154                 eneint[0][i + 1] =
1155                         eneint[0][i]
1156                         + 0.5 * (edat->s[1].es[i].sum + edat->s[3].es[i].sum) * Dt / edat->points[i];
1157                 eneint[1][i + 1] =
1158                         eneint[1][i]
1159                         + 0.5 * (edat->s[2].es[i].sum + edat->s[6].es[i].sum) * Dt / edat->points[i];
1160                 eneint[2][i + 1] =
1161                         eneint[2][i]
1162                         + 0.5 * (edat->s[5].es[i].sum + edat->s[7].es[i].sum) * Dt / edat->points[i];
1163             }
1164
1165             einstein_visco(eviscofn, eviscoifn, 3, edat->nframes + 1, eneint, Vaver, Temp, Dt, oenv);
1166
1167             for (i = 0; i < 3; i++)
1168             {
1169                 sfree(eneint[i]);
1170             }
1171             sfree(eneint);
1172
1173             /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1174             /* Do it for shear viscosity */
1175             std::strcpy(buf, "Shear Viscosity");
1176             low_do_autocorr(corrfn,
1177                             oenv,
1178                             buf,
1179                             edat->nframes,
1180                             3,
1181                             (edat->nframes + 1) / 2,
1182                             eneset,
1183                             Dt,
1184                             eacNormal,
1185                             1,
1186                             TRUE,
1187                             FALSE,
1188                             FALSE,
1189                             0.0,
1190                             0.0,
1191                             0);
1192
1193             /* Now for bulk viscosity */
1194             std::strcpy(buf, "Bulk Viscosity");
1195             low_do_autocorr(corrfn,
1196                             oenv,
1197                             buf,
1198                             edat->nframes,
1199                             1,
1200                             (edat->nframes + 1) / 2,
1201                             &(eneset[11]),
1202                             Dt,
1203                             eacNormal,
1204                             1,
1205                             TRUE,
1206                             FALSE,
1207                             FALSE,
1208                             0.0,
1209                             0.0,
1210                             0);
1211
1212             factor = (Vaver * 1e-26 / (gmx::c_boltzmann * Temp)) * Dt;
1213             fp     = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1214             xvgr_legend(fp, asize(leg), leg, oenv);
1215
1216             /* Use trapezium rule for integration */
1217             integral = 0;
1218             intBulk  = 0;
1219             nout     = get_acfnout();
1220             if ((nout < 2) || (nout >= edat->nframes / 2))
1221             {
1222                 nout = edat->nframes / 2;
1223             }
1224             for (i = 1; (i < nout); i++)
1225             {
1226                 integral += 0.5 * (eneset[0][i - 1] + eneset[0][i]) * factor;
1227                 intBulk += 0.5 * (eneset[11][i - 1] + eneset[11][i]) * factor;
1228                 fprintf(fp, "%10g  %10g  %10g\n", (i * Dt), integral, intBulk);
1229             }
1230             xvgrclose(fp);
1231
1232             for (i = 0; i < 12; i++)
1233             {
1234                 sfree(eneset[i]);
1235             }
1236             sfree(eneset);
1237         }
1238         else if (bCorr)
1239         {
1240             if (bFluct)
1241             {
1242                 std::strcpy(buf, "Autocorrelation of Energy Fluctuations");
1243             }
1244             else
1245             {
1246                 std::strcpy(buf, "Energy Autocorrelation");
1247             }
1248 #if 0
1249             do_autocorr(corrfn, oenv, buf, edat->nframes,
1250                         bSum ? 1                 : nset,
1251                         bSum ? &edat->s[nset-1].energyGroupPairTerms : eneset,
1252                         (delta_t/edat->nframes), eacNormal, FALSE);
1253 #endif
1254         }
1255     }
1256 }
1257
1258 static void print_time(FILE* fp, double t)
1259 {
1260     fprintf(fp, "%12.6f", t);
1261 }
1262
1263 static void print1(FILE* fp, gmx_bool bDp, real e)
1264 {
1265     if (bDp)
1266     {
1267         fprintf(fp, "  %16.12f", e);
1268     }
1269     else
1270     {
1271         fprintf(fp, "  %10.6f", e);
1272     }
1273 }
1274
1275 static void fec(const char*             ene2fn,
1276                 const char*             runavgfn,
1277                 real                    reftemp,
1278                 int                     nset,
1279                 const int               set[],
1280                 char*                   leg[],
1281                 enerdata_t*             edat,
1282                 double                  time[],
1283                 const gmx_output_env_t* oenv)
1284 {
1285     const char*  ravgleg[] = { "\\8D\\4E = E\\sB\\N-E\\sA\\N", "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N" };
1286     FILE*        fp;
1287     ener_file_t  enx;
1288     int          timecheck, nenergy, nenergy2, maxenergy;
1289     int          i, j;
1290     gmx_bool     bCont;
1291     real         aver, beta;
1292     real**       eneset2;
1293     double       dE, sum;
1294     gmx_enxnm_t* enm = nullptr;
1295     t_enxframe*  fr;
1296     char         buf[22];
1297
1298     /* read second energy file */
1299     snew(fr, 1);
1300     enm = nullptr;
1301     enx = open_enx(ene2fn, "r");
1302     do_enxnms(enx, &(fr->nre), &enm);
1303
1304     snew(eneset2, nset + 1);
1305     nenergy2  = 0;
1306     maxenergy = 0;
1307     timecheck = 0;
1308     do
1309     {
1310         /* This loop searches for the first frame (when -b option is given),
1311          * or when this has been found it reads just one energy frame
1312          */
1313         do
1314         {
1315             bCont = do_enx(enx, fr);
1316
1317             if (bCont)
1318             {
1319                 timecheck = check_times(fr->t);
1320             }
1321
1322         } while (bCont && (timecheck < 0));
1323
1324         /* Store energies for analysis afterwards... */
1325         if ((timecheck == 0) && bCont)
1326         {
1327             if (fr->nre > 0)
1328             {
1329                 if (nenergy2 >= maxenergy)
1330                 {
1331                     maxenergy += 1000;
1332                     for (i = 0; i <= nset; i++)
1333                     {
1334                         srenew(eneset2[i], maxenergy);
1335                     }
1336                 }
1337                 GMX_RELEASE_ASSERT(time != nullptr, "trying to dereference NULL time pointer");
1338
1339                 if (fr->t != time[nenergy2])
1340                 {
1341                     fprintf(stderr,
1342                             "\nWARNING time mismatch %g!=%g at frame %s\n",
1343                             fr->t,
1344                             time[nenergy2],
1345                             gmx_step_str(fr->step, buf));
1346                 }
1347                 for (i = 0; i < nset; i++)
1348                 {
1349                     eneset2[i][nenergy2] = fr->ener[set[i]].e;
1350                 }
1351                 nenergy2++;
1352             }
1353         }
1354     } while (bCont && (timecheck == 0));
1355
1356     /* check */
1357     if (edat->nframes != nenergy2)
1358     {
1359         fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n", edat->nframes, nenergy2);
1360     }
1361     nenergy = std::min(edat->nframes, nenergy2);
1362
1363     /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1364     fp = nullptr;
1365     if (runavgfn)
1366     {
1367         fp = xvgropen(runavgfn,
1368                       "Running average free energy difference",
1369                       "Time (" unit_time ")",
1370                       "\\8D\\4E (" unit_energy ")",
1371                       oenv);
1372         xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
1373     }
1374     fprintf(stdout, "\n%-24s %10s\n", "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
1375     sum  = 0;
1376     beta = 1.0 / (gmx::c_boltz * reftemp);
1377     for (i = 0; i < nset; i++)
1378     {
1379         if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
1380         {
1381             fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n", leg[i], enm[set[i]].name);
1382         }
1383         for (j = 0; j < nenergy; j++)
1384         {
1385             dE = eneset2[i][j] - edat->s[i].ener[j];
1386             sum += std::exp(-dE * beta);
1387             if (fp)
1388             {
1389                 fprintf(fp, "%10g %10g %10g\n", time[j], dE, -gmx::c_boltz * reftemp * std::log(sum / (j + 1)));
1390             }
1391         }
1392         aver = -gmx::c_boltz * reftemp * std::log(sum / nenergy);
1393         fprintf(stdout, "%-24s %10g\n", leg[i], aver);
1394     }
1395     if (fp)
1396     {
1397         xvgrclose(fp);
1398     }
1399     sfree(fr);
1400 }
1401
1402
1403 static void do_dhdl(t_enxframe*             fr,
1404                     const t_inputrec*       ir,
1405                     FILE**                  fp_dhdl,
1406                     const char*             filename,
1407                     gmx_bool                bDp,
1408                     int*                    blocks,
1409                     int*                    hists,
1410                     int*                    samples,
1411                     int*                    nlambdas,
1412                     const gmx_output_env_t* oenv)
1413 {
1414     const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1415     char        title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
1416     char        buf[STRLEN];
1417     int         nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
1418     int         i, j, k;
1419     /* coll data */
1420     double       temp = 0, start_time = 0, delta_time = 0, start_lambda = 0;
1421     static int   setnr             = 0;
1422     double*      native_lambda_vec = nullptr;
1423     const char** lambda_components = nullptr;
1424     int          n_lambda_vec      = 0;
1425     bool         firstPass         = true;
1426
1427     /* now count the blocks & handle the global dh data */
1428     for (i = 0; i < fr->nblock; i++)
1429     {
1430         if (fr->block[i].id == enxDHHIST)
1431         {
1432             nblock_hist++;
1433         }
1434         else if (fr->block[i].id == enxDH)
1435         {
1436             nblock_dh++;
1437         }
1438         else if (fr->block[i].id == enxDHCOLL)
1439         {
1440             nblock_dhcoll++;
1441             if ((fr->block[i].nsub < 1) || (fr->block[i].sub[0].type != xdr_datatype_double)
1442                 || (fr->block[i].sub[0].nr < 5))
1443             {
1444                 gmx_fatal(FARGS, "Unexpected block data");
1445             }
1446
1447             /* read the data from the DHCOLL block */
1448             temp         = fr->block[i].sub[0].dval[0];
1449             start_time   = fr->block[i].sub[0].dval[1];
1450             delta_time   = fr->block[i].sub[0].dval[2];
1451             start_lambda = fr->block[i].sub[0].dval[3];
1452             if (fr->block[i].nsub > 1)
1453             {
1454                 if (firstPass)
1455                 {
1456                     n_lambda_vec = fr->block[i].sub[1].ival[1];
1457                     snew(lambda_components, n_lambda_vec);
1458                     snew(native_lambda_vec, n_lambda_vec);
1459                     firstPass = false;
1460                 }
1461                 else
1462                 {
1463                     if (n_lambda_vec != fr->block[i].sub[1].ival[1])
1464                     {
1465                         gmx_fatal(FARGS, "Unexpected change of basis set in lambda");
1466                     }
1467                 }
1468                 for (j = 0; j < n_lambda_vec; j++)
1469                 {
1470                     native_lambda_vec[j] = fr->block[i].sub[0].dval[5 + j];
1471                     lambda_components[j] =
1472                             enumValueToStringSingular(static_cast<FreeEnergyPerturbationCouplingType>(
1473                                     fr->block[i].sub[1].ival[2 + j]));
1474                 }
1475             }
1476         }
1477     }
1478     // Clean up!
1479     sfree(native_lambda_vec);
1480     sfree(lambda_components);
1481     if (nblock_hist == 0 && nblock_dh == 0)
1482     {
1483         /* don't do anything */
1484         return;
1485     }
1486     if (nblock_hist > 0 && nblock_dh > 0)
1487     {
1488         gmx_fatal(FARGS,
1489                   "This energy file contains both histogram dhdl data and non-histogram dhdl data. "
1490                   "Don't know what to do.");
1491     }
1492     if (!*fp_dhdl)
1493     {
1494         if (nblock_dh > 0)
1495         {
1496             /* we have standard, non-histogram data --
1497                call open_dhdl to open the file */
1498             /* TODO this is an ugly hack that needs to be fixed: this will only
1499                work if the order of data is always the same and if we're
1500                only using the g_energy compiled with the mdrun that produced
1501                the ener.edr. */
1502             *fp_dhdl = open_dhdl(filename, ir, oenv);
1503         }
1504         else
1505         {
1506             sprintf(title, "N(%s)", deltag);
1507             sprintf(label_x, "%s (%s)", deltag, unit_energy);
1508             sprintf(label_y, "Samples");
1509             *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1510             sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
1511             xvgr_subtitle(*fp_dhdl, buf, oenv);
1512         }
1513     }
1514
1515     (*hists) += nblock_hist;
1516     (*blocks) += nblock_dh;
1517     (*nlambdas) = nblock_hist + nblock_dh;
1518
1519     /* write the data */
1520     if (nblock_hist > 0)
1521     {
1522         int64_t sum = 0;
1523         /* histograms */
1524         for (i = 0; i < fr->nblock; i++)
1525         {
1526             t_enxblock* blk = &(fr->block[i]);
1527             if (blk->id == enxDHHIST)
1528             {
1529                 double  foreign_lambda, dx;
1530                 int64_t x0;
1531                 int     nhist, derivative;
1532
1533                 /* check the block types etc. */
1534                 if ((blk->nsub < 2) || (blk->sub[0].type != xdr_datatype_double)
1535                     || (blk->sub[1].type != xdr_datatype_int64) || (blk->sub[0].nr < 2)
1536                     || (blk->sub[1].nr < 2))
1537                 {
1538                     gmx_fatal(FARGS, "Unexpected block data in file");
1539                 }
1540                 foreign_lambda = blk->sub[0].dval[0];
1541                 dx             = blk->sub[0].dval[1];
1542                 nhist          = blk->sub[1].lval[0];
1543                 derivative     = blk->sub[1].lval[1];
1544                 for (j = 0; j < nhist; j++)
1545                 {
1546                     const char* lg[1];
1547                     x0 = blk->sub[1].lval[2 + j];
1548
1549                     if (!derivative)
1550                     {
1551                         sprintf(legend, "N(%s(%s=%g) | %s=%g)", deltag, lambda, foreign_lambda, lambda, start_lambda);
1552                     }
1553                     else
1554                     {
1555                         sprintf(legend, "N(%s | %s=%g)", dhdl, lambda, start_lambda);
1556                     }
1557
1558                     lg[0] = legend;
1559                     xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1560                     setnr++;
1561                     for (k = 0; k < blk->sub[j + 2].nr; k++)
1562                     {
1563                         int    hist;
1564                         double xmin, xmax;
1565
1566                         hist = blk->sub[j + 2].ival[k];
1567                         xmin = (x0 + k) * dx;
1568                         xmax = (x0 + k + 1) * dx;
1569                         fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist, xmax, hist);
1570                         sum += hist;
1571                     }
1572                     /* multiple histogram data blocks in one histogram
1573                        mean that the second one is the reverse of the first one:
1574                        for dhdl derivatives, it's important to know both the
1575                        maximum and minimum values */
1576                     dx = -dx;
1577                 }
1578             }
1579         }
1580         (*samples) += static_cast<int>(sum / nblock_hist);
1581     }
1582     else
1583     {
1584         /* raw dh */
1585         int len = 0;
1586
1587         for (i = 0; i < fr->nblock; i++)
1588         {
1589             t_enxblock* blk = &(fr->block[i]);
1590             if (blk->id == enxDH)
1591             {
1592                 if (len == 0)
1593                 {
1594                     len = blk->sub[2].nr;
1595                 }
1596                 else
1597                 {
1598                     if (len != blk->sub[2].nr)
1599                     {
1600                         gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1601                     }
1602                 }
1603             }
1604         }
1605         (*samples) += len;
1606
1607         for (i = 0; i < len; i++)
1608         {
1609             double time = start_time + delta_time * i;
1610
1611             fprintf(*fp_dhdl, "%.4f ", time);
1612
1613             for (j = 0; j < fr->nblock; j++)
1614             {
1615                 t_enxblock* blk = &(fr->block[j]);
1616                 if (blk->id == enxDH)
1617                 {
1618                     double value;
1619                     if (blk->sub[2].type == xdr_datatype_float)
1620                     {
1621                         value = blk->sub[2].fval[i];
1622                     }
1623                     else
1624                     {
1625                         value = blk->sub[2].dval[i];
1626                     }
1627                     /* we need to decide which data type it is based on the count*/
1628
1629                     if (j == 1 && ir->bExpanded)
1630                     {
1631                         fprintf(*fp_dhdl, "%4d", static_cast<int>(value)); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */
1632                     }
1633                     else
1634                     {
1635                         if (bDp)
1636                         {
1637                             fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */
1638                         }
1639                         else
1640                         {
1641                             fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */
1642                         }
1643                     }
1644                 }
1645             }
1646             fprintf(*fp_dhdl, "\n");
1647         }
1648     }
1649 }
1650
1651
1652 int gmx_energy(int argc, char* argv[])
1653 {
1654     const char* desc[] = {
1655         "[THISMODULE] extracts energy components",
1656         "from an energy file. The user is prompted to interactively",
1657         "select the desired energy terms.[PAR]",
1658
1659         "Average, RMSD, and drift are calculated with full precision from the",
1660         "simulation (see printed manual). Drift is calculated by performing",
1661         "a least-squares fit of the data to a straight line. The reported total drift",
1662         "is the difference of the fit at the first and last point.",
1663         "An error estimate of the average is given based on a block averages",
1664         "over 5 blocks using the full-precision averages. The error estimate",
1665         "can be performed over multiple block lengths with the options",
1666         "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1667         "[BB]Note[bb] that in most cases the energy files contains averages over all",
1668         "MD steps, or over many more points than the number of frames in",
1669         "energy file. This makes the [THISMODULE] statistics output more accurate",
1670         "than the [REF].xvg[ref] output. When exact averages are not present in the energy",
1671         "file, the statistics mentioned above are simply over the single, per-frame",
1672         "energy values.[PAR]",
1673
1674         "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1675
1676         "Some fluctuation-dependent properties can be calculated provided",
1677         "the correct energy terms are selected, and that the command line option",
1678         "[TT]-fluct_props[tt] is given. The following properties",
1679         "will be computed:",
1680         "",
1681         "===============================  ===================",
1682         "Property                         Energy terms needed",
1683         "===============================  ===================",
1684         "Heat capacity C[SUB]p[sub] (NPT sims):    Enthalpy, Temp",
1685         "Heat capacity C[SUB]v[sub] (NVT sims):    Etot, Temp",
1686         "Thermal expansion coeff. (NPT):  Enthalpy, Vol, Temp",
1687         "Isothermal compressibility:      Vol, Temp",
1688         "Adiabatic bulk modulus:          Vol, Temp",
1689         "===============================  ===================",
1690         "",
1691         "You always need to set the number of molecules [TT]-nmol[tt].",
1692         "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1693         "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]",
1694
1695         "Option [TT]-odh[tt] extracts and plots the free energy data",
1696         "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1697         "from the [TT]ener.edr[tt] file.[PAR]",
1698
1699         "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1700         "difference with an ideal gas state::",
1701         "",
1702         "  [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT ",
1703         "  [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln]",
1704         "  [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT ",
1705         "  [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln]",
1706         "",
1707         "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1708         "the average is over the ensemble (or time in a trajectory).",
1709         "Note that this is in principle",
1710         "only correct when averaging over the whole (Boltzmann) ensemble",
1711         "and using the potential energy. This also allows for an entropy",
1712         "estimate using::",
1713         "",
1714         "  [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ",
1715         "  ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T",
1716         "  [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ",
1717         "  ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T",
1718         "",
1719
1720         "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1721         "difference is calculated::",
1722         "",
1723         "  dF = -kT ",
1724         "  [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub]) / ",
1725         "  kT[exp][chevron][SUB]A[sub][ln],",
1726         "",
1727         "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1728         "files, and the average is over the ensemble A. The running average",
1729         "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1730         "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1731
1732     };
1733     static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1734     static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1735     static int      nmol = 1, nbmin = 5, nbmax = 5;
1736     static real     reftemp = 300.0, ezero = 0;
1737     t_pargs         pa[] = {
1738         { "-fee", FALSE, etBOOL, { &bFee }, "Do a free energy estimate" },
1739         { "-fetemp",
1740           FALSE,
1741           etREAL,
1742           { &reftemp },
1743           "Reference temperature for free energy calculation" },
1744         { "-zero", FALSE, etREAL, { &ezero }, "Subtract a zero-point energy" },
1745         { "-sum",
1746           FALSE,
1747           etBOOL,
1748           { &bSum },
1749           "Sum the energy terms selected rather than display them all" },
1750         { "-dp", FALSE, etBOOL, { &bDp }, "Print energies in high precision" },
1751         { "-nbmin", FALSE, etINT, { &nbmin }, "Minimum number of blocks for error estimate" },
1752         { "-nbmax", FALSE, etINT, { &nbmax }, "Maximum number of blocks for error estimate" },
1753         { "-mutot",
1754           FALSE,
1755           etBOOL,
1756           { &bMutot },
1757           "Compute the total dipole moment from the components" },
1758         { "-aver",
1759           FALSE,
1760           etBOOL,
1761           { &bPrAll },
1762           "Also print the exact average and rmsd stored in the energy frames (only when 1 term is "
1763           "requested)" },
1764         { "-nmol",
1765           FALSE,
1766           etINT,
1767           { &nmol },
1768           "Number of molecules in your sample: the energies are divided by this number" },
1769         { "-fluct_props",
1770           FALSE,
1771           etBOOL,
1772           { &bFluctProps },
1773           "Compute properties based on energy fluctuations, like heat capacity" },
1774         { "-driftcorr",
1775           FALSE,
1776           etBOOL,
1777           { &bDriftCorr },
1778           "Useful only for calculations of fluctuation properties. The drift in the observables "
1779           "will be subtracted before computing the fluctuation properties." },
1780         { "-fluc",
1781           FALSE,
1782           etBOOL,
1783           { &bFluct },
1784           "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1785         { "-orinst", FALSE, etBOOL, { &bOrinst }, "Analyse instantaneous orientation data" },
1786         { "-ovec", FALSE, etBOOL, { &bOvec }, "Also plot the eigenvectors with [TT]-oten[tt]" }
1787     };
1788     static const char* setnm[] = { "Pres-XX", "Pres-XY",     "Pres-XZ", "Pres-YX",
1789                                    "Pres-YY", "Pres-YZ",     "Pres-ZX", "Pres-ZY",
1790                                    "Pres-ZZ", "Temperature", "Volume",  "Pressure" };
1791
1792     FILE*        out     = nullptr;
1793     FILE*        fp_dhdl = nullptr;
1794     ener_file_t  fp;
1795     int          timecheck = 0;
1796     enerdata_t   edat;
1797     gmx_enxnm_t* enm = nullptr;
1798     t_enxframe * frame, *fr = nullptr;
1799     int          cur = 0;
1800 #define NEXT (1 - cur)
1801     int               nre, nfr;
1802     int64_t           start_step;
1803     real              start_t;
1804     gmx_bool          bDHDL;
1805     gmx_bool          bFoundStart, bCont, bVisco;
1806     double            sum, dbl;
1807     double*           time = nullptr;
1808     real              Vaver;
1809     int *             set     = nullptr, i, j, nset, sss;
1810     gmx_bool*         bIsEner = nullptr;
1811     char**            leg     = nullptr;
1812     char              buf[256];
1813     gmx_output_env_t* oenv;
1814     int               dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
1815
1816     t_filenm fnm[] = {
1817         { efEDR, "-f", nullptr, ffREAD },        { efEDR, "-f2", nullptr, ffOPTRD },
1818         { efTPR, "-s", nullptr, ffOPTRD },       { efXVG, "-o", "energy", ffWRITE },
1819         { efXVG, "-viol", "violaver", ffOPTWR }, { efXVG, "-pairs", "pairs", ffOPTWR },
1820         { efXVG, "-corr", "enecorr", ffOPTWR },  { efXVG, "-vis", "visco", ffOPTWR },
1821         { efXVG, "-evisco", "evisco", ffOPTWR }, { efXVG, "-eviscoi", "eviscoi", ffOPTWR },
1822         { efXVG, "-ravg", "runavgdf", ffOPTWR }, { efXVG, "-odh", "dhdl", ffOPTWR }
1823     };
1824 #define NFILE asize(fnm)
1825     int      npargs;
1826     t_pargs* ppa;
1827
1828     npargs = asize(pa);
1829     ppa    = add_acf_pargs(&npargs, pa);
1830     if (!parse_common_args(
1831                 &argc, argv, PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END, NFILE, fnm, npargs, ppa, asize(desc), desc, 0, nullptr, &oenv))
1832     {
1833         sfree(ppa);
1834         return 0;
1835     }
1836
1837     bDHDL = opt2bSet("-odh", NFILE, fnm);
1838
1839     nset = 0;
1840
1841     snew(frame, 2);
1842     fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
1843     do_enxnms(fp, &nre, &enm);
1844
1845     Vaver = -1;
1846
1847     bVisco = opt2bSet("-vis", NFILE, fnm);
1848
1849     t_inputrec  irInstance;
1850     t_inputrec* ir = &irInstance;
1851
1852     if (!bDHDL)
1853     {
1854         if (bVisco)
1855         {
1856             nset = asize(setnm);
1857             snew(set, nset);
1858             /* This is nasty code... To extract Pres tensor, Volume and Temperature */
1859             for (j = 0; j < nset; j++)
1860             {
1861                 for (i = 0; i < nre; i++)
1862                 {
1863                     if (std::strstr(enm[i].name, setnm[j]))
1864                     {
1865                         set[j] = i;
1866                         break;
1867                     }
1868                 }
1869                 if (i == nre)
1870                 {
1871                     if (gmx_strcasecmp(setnm[j], "Volume") == 0)
1872                     {
1873                         printf("Enter the box volume (" unit_volume "): ");
1874                         if (1 != scanf("%lf", &dbl))
1875                         {
1876                             gmx_fatal(FARGS, "Error reading user input");
1877                         }
1878                         Vaver = dbl;
1879                     }
1880                     else
1881                     {
1882                         gmx_fatal(FARGS, "Could not find term %s for viscosity calculation", setnm[j]);
1883                     }
1884                 }
1885             }
1886         }
1887         else
1888         {
1889             set = select_by_name(nre, enm, &nset);
1890         }
1891         /* Print all the different units once */
1892         sprintf(buf, "(%s)", enm[set[0]].unit);
1893         for (i = 1; i < nset; i++)
1894         {
1895             for (j = 0; j < i; j++)
1896             {
1897                 if (std::strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
1898                 {
1899                     break;
1900                 }
1901             }
1902             if (j == i)
1903             {
1904                 std::strcat(buf, ", (");
1905                 std::strcat(buf, enm[set[i]].unit);
1906                 std::strcat(buf, ")");
1907             }
1908         }
1909         out = xvgropen(opt2fn("-o", NFILE, fnm), "GROMACS Energies", "Time (ps)", buf, oenv);
1910
1911         snew(leg, nset + 1);
1912         for (i = 0; (i < nset); i++)
1913         {
1914             leg[i] = enm[set[i]].name;
1915         }
1916         if (bSum)
1917         {
1918             leg[nset] = gmx_strdup("Sum");
1919             xvgr_legend(out, nset + 1, leg, oenv);
1920         }
1921         else
1922         {
1923             xvgr_legend(out, nset, leg, oenv);
1924         }
1925
1926         snew(bIsEner, nset);
1927         for (i = 0; (i < nset); i++)
1928         {
1929             bIsEner[i] = FALSE;
1930             for (j = 0; (j <= F_ETOT); j++)
1931             {
1932                 bIsEner[i] = bIsEner[i]
1933                              || (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
1934             }
1935         }
1936         if (bPrAll && nset > 1)
1937         {
1938             gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
1939         }
1940     }
1941     else if (bDHDL)
1942     {
1943         get_dhdl_parms(ftp2fn(efTPR, NFILE, fnm), ir);
1944     }
1945
1946     /* Initiate energies and set them to zero */
1947     edat.nsteps    = 0;
1948     edat.npoints   = 0;
1949     edat.nframes   = 0;
1950     edat.step      = nullptr;
1951     edat.steps     = nullptr;
1952     edat.points    = nullptr;
1953     edat.bHaveSums = TRUE;
1954     snew(edat.s, nset);
1955
1956     /* Initiate counters */
1957     bFoundStart = FALSE;
1958     start_step  = 0;
1959     start_t     = 0;
1960     do
1961     {
1962         /* This loop searches for the first frame (when -b option is given),
1963          * or when this has been found it reads just one energy frame
1964          */
1965         do
1966         {
1967             bCont = do_enx(fp, &(frame[NEXT]));
1968             if (bCont)
1969             {
1970                 timecheck = check_times(frame[NEXT].t);
1971             }
1972         } while (bCont && (timecheck < 0));
1973
1974         if ((timecheck == 0) && bCont)
1975         {
1976             /* We read a valid frame, so we can use it */
1977             fr = &(frame[NEXT]);
1978
1979             if (fr->nre > 0)
1980             {
1981                 /* The frame contains energies, so update cur */
1982                 cur = NEXT;
1983
1984                 if (edat.nframes % 1000 == 0)
1985                 {
1986                     srenew(edat.step, edat.nframes + 1000);
1987                     std::memset(&(edat.step[edat.nframes]), 0, 1000 * sizeof(edat.step[0]));
1988                     srenew(edat.steps, edat.nframes + 1000);
1989                     std::memset(&(edat.steps[edat.nframes]), 0, 1000 * sizeof(edat.steps[0]));
1990                     srenew(edat.points, edat.nframes + 1000);
1991                     std::memset(&(edat.points[edat.nframes]), 0, 1000 * sizeof(edat.points[0]));
1992
1993                     for (i = 0; i < nset; i++)
1994                     {
1995                         srenew(edat.s[i].ener, edat.nframes + 1000);
1996                         std::memset(&(edat.s[i].ener[edat.nframes]), 0, 1000 * sizeof(edat.s[i].ener[0]));
1997                         srenew(edat.s[i].es, edat.nframes + 1000);
1998                         std::memset(&(edat.s[i].es[edat.nframes]), 0, 1000 * sizeof(edat.s[i].es[0]));
1999                     }
2000                 }
2001
2002                 nfr            = edat.nframes;
2003                 edat.step[nfr] = fr->step;
2004
2005                 if (!bFoundStart)
2006                 {
2007                     bFoundStart = TRUE;
2008                     /* Initiate the previous step data */
2009                     start_step = fr->step;
2010                     start_t    = fr->t;
2011                     /* Initiate the energy sums */
2012                     edat.steps[nfr]  = 1;
2013                     edat.points[nfr] = 1;
2014                     for (i = 0; i < nset; i++)
2015                     {
2016                         sss                    = set[i];
2017                         edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2018                         edat.s[i].es[nfr].sum2 = 0;
2019                     }
2020                     edat.nsteps  = 1;
2021                     edat.npoints = 1;
2022                 }
2023                 else
2024                 {
2025                     edat.steps[nfr] = fr->nsteps;
2026
2027                     if (fr->nsum <= 1)
2028                     {
2029                         /* mdrun only calculated the energy at energy output
2030                          * steps. We don't need to check step intervals.
2031                          */
2032                         edat.points[nfr] = 1;
2033                         for (i = 0; i < nset; i++)
2034                         {
2035                             sss                    = set[i];
2036                             edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2037                             edat.s[i].es[nfr].sum2 = 0;
2038                         }
2039                         edat.npoints += 1;
2040                         edat.bHaveSums = FALSE;
2041                     }
2042                     else if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2043                     {
2044                         /* We have statistics  to the previous frame */
2045                         edat.points[nfr] = fr->nsum;
2046                         for (i = 0; i < nset; i++)
2047                         {
2048                             sss                    = set[i];
2049                             edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
2050                             edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2051                         }
2052                         edat.npoints += fr->nsum;
2053                     }
2054                     else
2055                     {
2056                         /* The interval does not match fr->nsteps:
2057                          * can not do exact averages.
2058                          */
2059                         edat.bHaveSums = FALSE;
2060                     }
2061
2062                     edat.nsteps = fr->step - start_step + 1;
2063                 }
2064                 for (i = 0; i < nset; i++)
2065                 {
2066                     edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2067                 }
2068             }
2069             /*
2070              * Store energies for analysis afterwards...
2071              */
2072             if (!bDHDL && (fr->nre > 0))
2073             {
2074                 if (edat.nframes % 1000 == 0)
2075                 {
2076                     srenew(time, edat.nframes + 1000);
2077                 }
2078                 time[edat.nframes] = fr->t;
2079                 edat.nframes++;
2080             }
2081             if (bDHDL)
2082             {
2083                 do_dhdl(fr, ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2084             }
2085
2086             /*******************************************
2087              * E N E R G I E S
2088              *******************************************/
2089             else
2090             {
2091                 if (fr->nre > 0)
2092                 {
2093                     if (bPrAll)
2094                     {
2095                         /* We skip frames with single points (usually only the first frame),
2096                          * since they would result in an average plot with outliers.
2097                          */
2098                         if (fr->nsum > 1)
2099                         {
2100                             print_time(out, fr->t);
2101                             print1(out, bDp, fr->ener[set[0]].e);
2102                             print1(out, bDp, fr->ener[set[0]].esum / fr->nsum);
2103                             print1(out, bDp, std::sqrt(fr->ener[set[0]].eav / fr->nsum));
2104                             fprintf(out, "\n");
2105                         }
2106                     }
2107                     else
2108                     {
2109                         print_time(out, fr->t);
2110                         if (bSum)
2111                         {
2112                             sum = 0;
2113                             for (i = 0; i < nset; i++)
2114                             {
2115                                 sum += fr->ener[set[i]].e;
2116                             }
2117                             print1(out, bDp, sum / nmol - ezero);
2118                         }
2119                         else
2120                         {
2121                             for (i = 0; (i < nset); i++)
2122                             {
2123                                 if (bIsEner[i])
2124                                 {
2125                                     print1(out, bDp, (fr->ener[set[i]].e) / nmol - ezero);
2126                                 }
2127                                 else
2128                                 {
2129                                     print1(out, bDp, fr->ener[set[i]].e);
2130                                 }
2131                             }
2132                         }
2133                         fprintf(out, "\n");
2134                     }
2135                 }
2136             }
2137         }
2138     } while (bCont && (timecheck == 0));
2139
2140     fprintf(stderr, "\n");
2141     done_ener_file(fp);
2142     if (out)
2143     {
2144         xvgrclose(out);
2145     }
2146
2147     if (bDHDL)
2148     {
2149         if (fp_dhdl)
2150         {
2151             gmx_fio_fclose(fp_dhdl);
2152             printf("\n\nWrote %d lambda values with %d samples as ", dh_lambdas, dh_samples);
2153             if (dh_hists > 0)
2154             {
2155                 printf("%d dH histograms ", dh_hists);
2156             }
2157             if (dh_blocks > 0)
2158             {
2159                 printf("%d dH data blocks ", dh_blocks);
2160             }
2161             printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2162         }
2163         else
2164         {
2165             gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2166         }
2167     }
2168     else
2169     {
2170         double dt = (frame[cur].t - start_t) / (edat.nframes - 1);
2171         analyse_ener(opt2bSet("-corr", NFILE, fnm),
2172                      opt2fn("-corr", NFILE, fnm),
2173                      opt2fn("-evisco", NFILE, fnm),
2174                      opt2fn("-eviscoi", NFILE, fnm),
2175                      bFee,
2176                      bSum,
2177                      bFluct,
2178                      bVisco,
2179                      opt2fn("-vis", NFILE, fnm),
2180                      nmol,
2181                      start_step,
2182                      start_t,
2183                      frame[cur].step,
2184                      frame[cur].t,
2185                      reftemp,
2186                      &edat,
2187                      nset,
2188                      set,
2189                      bIsEner,
2190                      leg,
2191                      enm,
2192                      Vaver,
2193                      ezero,
2194                      nbmin,
2195                      nbmax,
2196                      oenv);
2197         if (bFluctProps)
2198         {
2199             calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat, nbmin, nbmax);
2200         }
2201     }
2202     if (opt2bSet("-f2", NFILE, fnm))
2203     {
2204         fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm), reftemp, nset, set, leg, &edat, time, oenv);
2205     }
2206     // Clean up!
2207     done_enerdata_t(nset, &edat);
2208     sfree(time);
2209     free_enxframe(&frame[0]);
2210     free_enxframe(&frame[1]);
2211     sfree(frame);
2212     free_enxnms(nre, enm);
2213     sfree(ppa);
2214     sfree(set);
2215     sfree(leg);
2216     sfree(bIsEner);
2217     {
2218         const char* nxy = "-nxy";
2219
2220         do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2221         do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2222         do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);
2223     }
2224     output_env_done(oenv);
2225
2226     return 0;
2227 }