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