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