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