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