2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
46 #include "gromacs/commandline/pargs.h"
47 #include "gromacs/commandline/viewit.h"
48 #include "gromacs/correlationfunctions/autocorr.h"
49 #include "gromacs/fileio/enxio.h"
50 #include "gromacs/fileio/gmxfio.h"
51 #include "gromacs/fileio/tpxio.h"
52 #include "gromacs/fileio/trxio.h"
53 #include "gromacs/fileio/xvgr.h"
54 #include "gromacs/gmxana/gmx_ana.h"
55 #include "gromacs/gmxana/gstat.h"
56 #include "gromacs/math/functions.h"
57 #include "gromacs/math/units.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/mdlib/energyoutput.h"
60 #include "gromacs/mdtypes/inputrec.h"
61 #include "gromacs/mdtypes/md_enums.h"
62 #include "gromacs/topology/ifunc.h"
63 #include "gromacs/topology/mtop_lookup.h"
64 #include "gromacs/topology/mtop_util.h"
65 #include "gromacs/topology/topology.h"
66 #include "gromacs/trajectory/energyframe.h"
67 #include "gromacs/utility/arraysize.h"
68 #include "gromacs/utility/cstringutil.h"
69 #include "gromacs/utility/fatalerror.h"
70 #include "gromacs/utility/gmxassert.h"
71 #include "gromacs/utility/pleasecite.h"
72 #include "gromacs/utility/smalloc.h"
73 #include "gromacs/utility/strconvert.h"
75 static const int NOTSET = -23451;
106 static void done_enerdata_t(int nset, enerdata_t* edat)
111 for (int i = 0; i < nset; i++)
113 sfree(edat->s[i].ener);
114 sfree(edat->s[i].es);
119 static void chomp(char* buf)
121 int len = std::strlen(buf);
123 while ((len > 0) && (buf[len - 1] == '\n'))
130 static int* select_by_name(int nre, gmx_enxnm_t* nm, int* nset)
133 int k, kk, j, i, nmatch, nind, nss;
135 gmx_bool bEOF, bVerbose = TRUE, bLong = FALSE;
136 char * ptr, buf[STRLEN];
137 const char* fm4 = "%3d %-14s";
138 const char* fm2 = "%3d %-34s";
139 char** newnm = nullptr;
141 if ((getenv("GMX_ENER_VERBOSE")) != nullptr)
146 fprintf(stderr, "\n");
147 fprintf(stderr, "Select the terms you want from the following list by\n");
148 fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
149 fprintf(stderr, "End your selection with an empty line or a zero.\n");
150 fprintf(stderr, "-------------------------------------------------------------------\n");
154 for (k = 0; k < nre; k++)
156 newnm[k] = gmx_strdup(nm[k].name);
157 /* Insert dashes in all the names */
158 while ((ptr = std::strchr(newnm[k], ' ')) != nullptr)
168 fprintf(stderr, "\n");
171 for (kk = k; kk < k + 4; kk++)
173 if (kk < nre && std::strlen(nm[kk].name) > 14)
181 fprintf(stderr, " ");
185 fprintf(stderr, fm4, k + 1, newnm[k]);
194 fprintf(stderr, fm2, k + 1, newnm[k]);
205 fprintf(stderr, "\n\n");
211 while (!bEOF && (fgets2(buf, STRLEN - 1, stdin)))
213 /* Remove newlines */
219 /* Empty line means end of input */
220 bEOF = (std::strlen(buf) == 0);
228 /* First try to read an integer */
229 nss = sscanf(ptr, "%d", &nind);
232 /* Zero means end of input */
237 else if ((1 <= nind) && (nind <= nre))
243 fprintf(stderr, "number %d is out of range\n", nind);
248 /* Now try to read a string */
250 for (nind = 0; nind < nre; nind++)
252 if (gmx_strcasecmp(newnm[nind], ptr) == 0)
260 i = std::strlen(ptr);
262 for (nind = 0; nind < nre; nind++)
264 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
272 fprintf(stderr, "String '%s' does not match anything\n", ptr);
277 /* Look for the first space, and remove spaces from there */
278 if ((ptr = std::strchr(ptr, ' ')) != nullptr)
282 } while (!bEOF && ((ptr != nullptr) && (std::strlen(ptr) > 0)));
287 for (i = (*nset) = 0; (i < nre); i++)
299 gmx_fatal(FARGS, "No energy terms selected");
302 for (i = 0; (i < nre); i++)
311 static void get_dhdl_parms(const char* topnm, t_inputrec* ir)
317 /* all we need is the ir to be able to write the label */
318 read_tpx(topnm, ir, box, &natoms, nullptr, nullptr, &mtop);
321 static void einstein_visco(const char* fn,
329 const gmx_output_env_t* oenv)
332 double av[4], avold[4];
338 for (i = 0; i <= nsets; i++)
342 fp0 = xvgropen(fni, "Shear viscosity integral", "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
343 fp1 = xvgropen(fn, "Shear viscosity using Einstein relation", "Time (ps)",
344 "(kg m\\S-1\\N s\\S-1\\N)", oenv);
345 for (i = 0; i < nf4; i++)
347 for (m = 0; m <= nsets; m++)
351 for (j = 0; j < nint - i; j++)
353 for (m = 0; m < nsets; m++)
355 di = gmx::square(eneint[m][j + i] - eneint[m][j]);
358 av[nsets] += di / nsets;
361 /* Convert to SI for the viscosity */
362 fac = (V * NANO * NANO * NANO * PICO * 1e10) / (2 * BOLTZMANN * T) / (nint - i);
363 fprintf(fp0, "%10g", i * dt);
364 for (m = 0; (m <= nsets); m++)
367 fprintf(fp0, " %10g", av[m]);
370 fprintf(fp1, "%10g", (i + 0.5) * dt);
371 for (m = 0; (m <= nsets); m++)
373 fprintf(fp1, " %10g", (av[m] - avold[m]) / dt);
398 static void clear_ee_sum(ee_sum_t* ees)
406 static void add_ee_sum(ee_sum_t* ees, double sum, int np)
412 static void add_ee_av(ee_sum_t* ees)
416 av = ees->sum / ees->np;
418 ees->sav2 += av * av;
423 static double calc_ee2(int nb, ee_sum_t* ees)
425 return (ees->sav2 / nb - gmx::square(ees->sav / nb)) / (nb - 1);
428 static void set_ee_av(ener_ee_t* eee)
432 char buf[STEPSTRSIZE];
433 fprintf(debug, "Storing average for err.est.: %s steps\n", gmx_step_str(eee->nst, buf));
435 add_ee_av(&eee->sum);
437 if (eee->b == 1 || eee->nst < eee->nst_min)
439 eee->nst_min = eee->nst;
444 static void calc_averages(int nset, enerdata_t* edat, int nbmin, int nbmax)
447 double sum, sum2, sump, see2;
448 int64_t np, p, bound_nb;
452 double x, sx, sy, sxx, sxy;
455 /* Check if we have exact statistics over all points */
456 for (i = 0; i < nset; i++)
459 ed->bExactStat = FALSE;
462 /* All energy file sum entries 0 signals no exact sums.
463 * But if all energy values are 0, we still have exact sums.
466 for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
468 if (ed->ener[i] != 0)
472 ed->bExactStat = (ed->es[f].sum != 0);
476 ed->bExactStat = TRUE;
481 snew(eee, nbmax + 1);
482 for (i = 0; i < nset; i++)
493 for (nb = nbmin; nb <= nbmax; nb++)
496 clear_ee_sum(&eee[nb].sum);
500 for (f = 0; f < edat->nframes; f++)
506 /* Add the sum and the sum of variances to the totals. */
512 sum2 += gmx::square(sum / np - (sum + es->sum) / (np + p)) * np * (np + p) / p;
517 /* Add a single value to the sum and sum of squares. */
520 sum2 += gmx::square(sump);
523 /* sum has to be increased after sum2 */
527 /* For the linear regression use variance 1/p.
528 * Note that sump is the sum, not the average, so we don't need p*.
530 x = edat->step[f] - 0.5 * (edat->steps[f] - 1);
536 for (nb = nbmin; nb <= nbmax; nb++)
538 /* Check if the current end step is closer to the desired
539 * block boundary than the next end step.
541 bound_nb = (edat->step[0] - 1) * nb + edat->nsteps * (eee[nb].b + 1);
542 if (eee[nb].nst > 0 && bound_nb - edat->step[f - 1] * nb < edat->step[f] * nb - bound_nb)
552 eee[nb].nst += edat->step[f] - edat->step[f - 1];
556 add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
560 add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
562 bound_nb = (edat->step[0] - 1) * nb + edat->nsteps * (eee[nb].b + 1);
563 if (edat->step[f] * nb >= bound_nb)
570 edat->s[i].av = sum / np;
573 edat->s[i].rmsd = std::sqrt(sum2 / np);
577 edat->s[i].rmsd = std::sqrt(sum2 / np - gmx::square(edat->s[i].av));
580 if (edat->nframes > 1)
582 edat->s[i].slope = (np * sxy - sx * sy) / (np * sxx - sx * sx);
586 edat->s[i].slope = 0;
591 for (nb = nbmin; nb <= nbmax; nb++)
593 /* Check if we actually got nb blocks and if the smallest
594 * block is not shorter than 80% of the average.
598 char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
599 fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n", nb,
600 eee[nb].b, gmx_step_str(eee[nb].nst_min, buf1), gmx_step_str(edat->nsteps, buf2));
602 if (eee[nb].b == nb && 5 * nb * eee[nb].nst_min >= 4 * edat->nsteps)
604 see2 += calc_ee2(nb, &eee[nb].sum);
610 edat->s[i].ee = std::sqrt(see2 / nee);
620 static enerdata_t* calc_sum(int nset, enerdata_t* edat, int nbmin, int nbmax)
631 snew(s->ener, esum->nframes);
632 snew(s->es, esum->nframes);
634 s->bExactStat = TRUE;
636 for (i = 0; i < nset; i++)
638 if (!edat->s[i].bExactStat)
640 s->bExactStat = FALSE;
642 s->slope += edat->s[i].slope;
645 for (f = 0; f < edat->nframes; f++)
648 for (i = 0; i < nset; i++)
650 sum += edat->s[i].ener[f];
654 for (i = 0; i < nset; i++)
656 sum += edat->s[i].es[f].sum;
662 calc_averages(1, esum, nbmin, nbmax);
667 static void ee_pr(double ee, int buflen, char* buf)
669 snprintf(buf, buflen, "%s", "--");
672 /* Round to two decimals by printing. */
674 snprintf(tmp, sizeof(tmp), "%.1e", ee);
675 double rnd = gmx::doubleFromString(tmp);
676 snprintf(buf, buflen, "%g", rnd);
680 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t* edat)
682 /* Remove the drift by performing a fit to y = ax+b.
683 Uses 5 iterations. */
687 edat->npoints = edat->nframes;
688 edat->nsteps = edat->nframes;
690 for (k = 0; (k < 5); k++)
692 for (i = 0; (i < nset); i++)
694 delta = edat->s[i].slope * dt;
696 if (nullptr != debug)
698 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
701 for (j = 0; (j < edat->nframes); j++)
703 edat->s[i].ener[j] -= j * delta;
704 edat->s[i].es[j].sum = 0;
705 edat->s[i].es[j].sum2 = 0;
708 calc_averages(nset, edat, nbmin, nbmax);
712 static void calc_fluctuation_props(FILE* fp,
723 double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, varet;
733 const char* my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
736 NANO3 = NANO * NANO * NANO;
740 "\nYou may want to use the -driftcorr flag in order to correct\n"
741 "for spurious drift in the graphs. Note that this is not\n"
742 "a substitute for proper equilibration and sampling!\n");
746 remove_drift(nset, nbmin, nbmax, dt, edat);
748 for (i = 0; (i < eNR); i++)
750 for (ii[i] = 0; (ii[i] < nset && (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++) {}
752 fprintf(fp,"Found %s data.\n",my_ener[i]);
754 /* Compute it all! */
755 alpha = kappa = cp = dcp = cv = NOTSET;
759 if (ii[eTemp] < nset)
761 tt = edat->s[ii[eTemp]].av;
765 if ((ii[eVol] < nset) && (ii[eTemp] < nset))
767 vv = edat->s[ii[eVol]].av * NANO3;
768 varv = gmx::square(edat->s[ii[eVol]].rmsd * NANO3);
769 kappa = (varv / vv) / (BOLTZMANN * tt);
773 if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
775 hh = KILO * edat->s[ii[eEnth]].av / AVOGADRO;
776 varh = gmx::square(KILO * edat->s[ii[eEnth]].rmsd / AVOGADRO);
777 cp = AVOGADRO * ((varh / nmol) / (BOLTZMANN * tt * tt));
780 if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
782 /* Only compute cv in constant volume runs, which we can test
783 by checking whether the enthalpy was computed.
785 varet = gmx::square(edat->s[ii[eEtot]].rmsd);
786 cv = KILO * ((varet / nmol) / (BOLTZ * tt * tt));
789 if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
791 double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
792 vh_sum = v_sum = h_sum = 0;
793 for (j = 0; (j < edat->nframes); j++)
795 v = edat->s[ii[eVol]].ener[j] * NANO3;
796 h = KILO * edat->s[ii[eEnth]].ener[j] / AVOGADRO;
801 vh_aver = vh_sum / edat->nframes;
802 v_aver = v_sum / edat->nframes;
803 h_aver = h_sum / edat->nframes;
804 alpha = (vh_aver - v_aver * h_aver) / (v_aver * BOLTZMANN * tt * tt);
805 dcp = (v_aver * AVOGADRO / nmol) * tt * gmx::square(alpha) / (kappa);
812 fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n", nmol);
814 fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
815 fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
816 fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
817 fprintf(fp, "please use the g_dos program.\n\n");
819 "WARNING: Please verify that your simulations are converged and perform\n"
820 "a block-averaging error analysis (not implemented in g_energy yet)\n");
822 if (debug != nullptr)
826 fprintf(fp, "varv = %10g (m^6)\n", varv * AVOGADRO / nmol);
831 fprintf(fp, "Volume = %10g m^3/mol\n", vv * AVOGADRO / nmol);
835 fprintf(fp, "Enthalpy = %10g kJ/mol\n",
836 hh * AVOGADRO / (KILO * nmol));
840 fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n", alpha);
844 fprintf(fp, "Isothermal Compressibility Kappa = %10g (m^3/J)\n", kappa);
845 fprintf(fp, "Adiabatic bulk modulus = %10g (J/m^3)\n", 1.0 / kappa);
849 fprintf(fp, "Heat capacity at constant pressure Cp = %10g J/(mol K)\n", cp);
853 fprintf(fp, "Heat capacity at constant volume Cv = %10g J/(mol K)\n", cv);
857 fprintf(fp, "Cp-Cv = %10g J/(mol K)\n", dcp);
859 please_cite(fp, "Allen1987a");
863 fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
867 static void analyse_ener(gmx_bool bCorr,
869 const char* eviscofn,
870 const char* eviscoifn,
885 const gmx_bool* bIsEner,
892 const gmx_output_env_t* oenv)
895 /* Check out the printed manual for equations! */
896 double Dt, aver, stddev, errest, delta_t, totaldrift;
897 enerdata_t* esum = nullptr;
898 real integral, intBulk, Temp = 0, Pres = 0;
899 real pr_aver, pr_stddev, pr_errest;
900 double beta = 0, expE, expEtot, *fee = nullptr;
902 int nexact, nnotexact;
904 char buf[256], eebuf[100];
906 nsteps = step - start_step + 1;
909 fprintf(stdout, "Not enough steps (%s) for statistics\n", gmx_step_str(nsteps, buf));
913 /* Calculate the time difference */
914 delta_t = t - start_t;
916 fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
917 gmx_step_str(nsteps, buf), start_t, t, nset);
919 calc_averages(nset, edat, nbmin, nbmax);
923 esum = calc_sum(nset, edat, nbmin, nbmax);
926 if (!edat->bHaveSums)
935 for (i = 0; (i < nset); i++)
937 if (edat->s[i].bExactStat)
950 fprintf(stdout, "All statistics are over %s points\n", gmx_step_str(edat->npoints, buf));
952 else if (nexact == 0 || edat->npoints == edat->nframes)
954 fprintf(stdout, "All statistics are over %d points (frames)\n", edat->nframes);
958 fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
959 for (i = 0; (i < nset); i++)
961 if (!edat->s[i].bExactStat)
963 fprintf(stdout, " '%s'", leg[i]);
966 fprintf(stdout, " %s has statistics over %d points (frames)\n",
967 nnotexact == 1 ? "is" : "are", edat->nframes);
968 fprintf(stdout, "All other statistics are over %s points\n", gmx_step_str(edat->npoints, buf));
970 fprintf(stdout, "\n");
972 fprintf(stdout, "%-24s %10s %10s %10s %10s", "Energy", "Average", "Err.Est.", "RMSD",
976 fprintf(stdout, " %10s\n", "-kT ln<e^(E/kT)>");
980 fprintf(stdout, "\n");
983 "-------------------------------------------------------------------------------"
986 /* Initiate locals, only used with -sum */
990 beta = 1.0 / (BOLTZ * reftemp);
993 for (i = 0; (i < nset); i++)
995 aver = edat->s[i].av;
996 stddev = edat->s[i].rmsd;
997 errest = edat->s[i].ee;
1002 for (j = 0; (j < edat->nframes); j++)
1004 expE += std::exp(beta * (edat->s[i].ener[j] - aver) / nmol);
1008 expEtot += expE / edat->nframes;
1011 fee[i] = std::log(expE / edat->nframes) / beta + aver / nmol;
1013 if (std::strstr(leg[i], "empera") != nullptr)
1017 else if (std::strstr(leg[i], "olum") != nullptr)
1021 else if (std::strstr(leg[i], "essure") != nullptr)
1027 pr_aver = aver / nmol - ezero;
1028 pr_stddev = stddev / nmol;
1029 pr_errest = errest / nmol;
1038 /* Multiply the slope in steps with the number of steps taken */
1039 totaldrift = (edat->nsteps - 1) * edat->s[i].slope;
1045 ee_pr(pr_errest, sizeof(eebuf), eebuf);
1046 fprintf(stdout, "%-24s %10g %10s %10g %10g", leg[i], pr_aver, eebuf, pr_stddev, totaldrift);
1049 fprintf(stdout, " %10g", fee[i]);
1052 fprintf(stdout, " (%s)\n", enm[set[i]].unit);
1056 for (j = 0; (j < edat->nframes); j++)
1058 edat->s[i].ener[j] -= aver;
1064 totaldrift = (edat->nsteps - 1) * esum->s[0].slope;
1065 ee_pr(esum->s[0].ee / nmol, sizeof(eebuf), eebuf);
1066 fprintf(stdout, "%-24s %10g %10s %10s %10g (%s)", "Total", esum->s[0].av / nmol, eebuf,
1067 "--", totaldrift / nmol, enm[set[0]].unit);
1068 /* pr_aver,pr_stddev,a,totaldrift */
1071 fprintf(stdout, " %10g %10g\n", std::log(expEtot) / beta + esum->s[0].av / nmol,
1072 std::log(expEtot) / beta);
1076 fprintf(stdout, "\n");
1080 /* Do correlation function */
1081 if (edat->nframes > 1)
1083 Dt = delta_t / (edat->nframes - 1);
1091 const char* leg[] = { "Shear", "Bulk" };
1096 /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1098 /* Symmetrise tensor! (and store in first three elements)
1099 * And subtract average pressure!
1102 for (i = 0; i < 12; i++)
1104 snew(eneset[i], edat->nframes);
1106 for (i = 0; (i < edat->nframes); i++)
1108 eneset[0][i] = 0.5 * (edat->s[1].ener[i] + edat->s[3].ener[i]);
1109 eneset[1][i] = 0.5 * (edat->s[2].ener[i] + edat->s[6].ener[i]);
1110 eneset[2][i] = 0.5 * (edat->s[5].ener[i] + edat->s[7].ener[i]);
1111 for (j = 3; j <= 11; j++)
1113 eneset[j][i] = edat->s[j].ener[i];
1115 eneset[11][i] -= Pres;
1118 /* Determine integrals of the off-diagonal pressure elements */
1120 for (i = 0; i < 3; i++)
1122 snew(eneint[i], edat->nframes + 1);
1127 for (i = 0; i < edat->nframes; i++)
1131 + 0.5 * (edat->s[1].es[i].sum + edat->s[3].es[i].sum) * Dt / edat->points[i];
1134 + 0.5 * (edat->s[2].es[i].sum + edat->s[6].es[i].sum) * Dt / edat->points[i];
1137 + 0.5 * (edat->s[5].es[i].sum + edat->s[7].es[i].sum) * Dt / edat->points[i];
1140 einstein_visco(eviscofn, eviscoifn, 3, edat->nframes + 1, eneint, Vaver, Temp, Dt, oenv);
1142 for (i = 0; i < 3; i++)
1148 /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1149 /* Do it for shear viscosity */
1150 std::strcpy(buf, "Shear Viscosity");
1151 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3, (edat->nframes + 1) / 2, eneset,
1152 Dt, eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1154 /* Now for bulk viscosity */
1155 std::strcpy(buf, "Bulk Viscosity");
1156 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1, (edat->nframes + 1) / 2,
1157 &(eneset[11]), Dt, eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1159 factor = (Vaver * 1e-26 / (BOLTZMANN * Temp)) * Dt;
1160 fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1161 xvgr_legend(fp, asize(leg), leg, oenv);
1163 /* Use trapezium rule for integration */
1166 nout = get_acfnout();
1167 if ((nout < 2) || (nout >= edat->nframes / 2))
1169 nout = edat->nframes / 2;
1171 for (i = 1; (i < nout); i++)
1173 integral += 0.5 * (eneset[0][i - 1] + eneset[0][i]) * factor;
1174 intBulk += 0.5 * (eneset[11][i - 1] + eneset[11][i]) * factor;
1175 fprintf(fp, "%10g %10g %10g\n", (i * Dt), integral, intBulk);
1179 for (i = 0; i < 12; i++)
1189 std::strcpy(buf, "Autocorrelation of Energy Fluctuations");
1193 std::strcpy(buf, "Energy Autocorrelation");
1196 do_autocorr(corrfn, oenv, buf, edat->nframes,
1198 bSum ? &edat->s[nset-1].ener : eneset,
1199 (delta_t/edat->nframes), eacNormal, FALSE);
1205 static void print_time(FILE* fp, double t)
1207 fprintf(fp, "%12.6f", t);
1210 static void print1(FILE* fp, gmx_bool bDp, real e)
1214 fprintf(fp, " %16.12f", e);
1218 fprintf(fp, " %10.6f", e);
1222 static void fec(const char* ene2fn,
1223 const char* runavgfn,
1230 const gmx_output_env_t* oenv)
1232 const char* ravgleg[] = { "\\8D\\4E = E\\sB\\N-E\\sA\\N", "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N" };
1235 int timecheck, nenergy, nenergy2, maxenergy;
1241 gmx_enxnm_t* enm = nullptr;
1245 /* read second energy file */
1248 enx = open_enx(ene2fn, "r");
1249 do_enxnms(enx, &(fr->nre), &enm);
1251 snew(eneset2, nset + 1);
1257 /* This loop searches for the first frame (when -b option is given),
1258 * or when this has been found it reads just one energy frame
1262 bCont = do_enx(enx, fr);
1266 timecheck = check_times(fr->t);
1269 } while (bCont && (timecheck < 0));
1271 /* Store energies for analysis afterwards... */
1272 if ((timecheck == 0) && bCont)
1276 if (nenergy2 >= maxenergy)
1279 for (i = 0; i <= nset; i++)
1281 srenew(eneset2[i], maxenergy);
1284 GMX_RELEASE_ASSERT(time != nullptr, "trying to dereference NULL time pointer");
1286 if (fr->t != time[nenergy2])
1288 fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n", fr->t,
1289 time[nenergy2], gmx_step_str(fr->step, buf));
1291 for (i = 0; i < nset; i++)
1293 eneset2[i][nenergy2] = fr->ener[set[i]].e;
1298 } while (bCont && (timecheck == 0));
1301 if (edat->nframes != nenergy2)
1303 fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n", edat->nframes, nenergy2);
1305 nenergy = std::min(edat->nframes, nenergy2);
1307 /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1311 fp = xvgropen(runavgfn, "Running average free energy difference", "Time (" unit_time ")",
1312 "\\8D\\4E (" unit_energy ")", oenv);
1313 xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
1315 fprintf(stdout, "\n%-24s %10s\n", "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
1317 beta = 1.0 / (BOLTZ * reftemp);
1318 for (i = 0; i < nset; i++)
1320 if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
1322 fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n", leg[i], enm[set[i]].name);
1324 for (j = 0; j < nenergy; j++)
1326 dE = eneset2[i][j] - edat->s[i].ener[j];
1327 sum += std::exp(-dE * beta);
1330 fprintf(fp, "%10g %10g %10g\n", time[j], dE, -BOLTZ * reftemp * std::log(sum / (j + 1)));
1333 aver = -BOLTZ * reftemp * std::log(sum / nenergy);
1334 fprintf(stdout, "%-24s %10g\n", leg[i], aver);
1344 static void do_dhdl(t_enxframe* fr,
1345 const t_inputrec* ir,
1347 const char* filename,
1353 const gmx_output_env_t* oenv)
1355 const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1356 char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
1358 int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
1361 double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0;
1362 static int setnr = 0;
1363 double* native_lambda_vec = nullptr;
1364 const char** lambda_components = nullptr;
1365 int n_lambda_vec = 0;
1366 bool firstPass = true;
1368 /* now count the blocks & handle the global dh data */
1369 for (i = 0; i < fr->nblock; i++)
1371 if (fr->block[i].id == enxDHHIST)
1375 else if (fr->block[i].id == enxDH)
1379 else if (fr->block[i].id == enxDHCOLL)
1382 if ((fr->block[i].nsub < 1) || (fr->block[i].sub[0].type != xdr_datatype_double)
1383 || (fr->block[i].sub[0].nr < 5))
1385 gmx_fatal(FARGS, "Unexpected block data");
1388 /* read the data from the DHCOLL block */
1389 temp = fr->block[i].sub[0].dval[0];
1390 start_time = fr->block[i].sub[0].dval[1];
1391 delta_time = fr->block[i].sub[0].dval[2];
1392 start_lambda = fr->block[i].sub[0].dval[3];
1393 if (fr->block[i].nsub > 1)
1397 n_lambda_vec = fr->block[i].sub[1].ival[1];
1398 snew(lambda_components, n_lambda_vec);
1399 snew(native_lambda_vec, n_lambda_vec);
1404 if (n_lambda_vec != fr->block[i].sub[1].ival[1])
1406 gmx_fatal(FARGS, "Unexpected change of basis set in lambda");
1409 for (j = 0; j < n_lambda_vec; j++)
1411 native_lambda_vec[j] = fr->block[i].sub[0].dval[5 + j];
1412 lambda_components[j] = efpt_singular_names[fr->block[i].sub[1].ival[2 + j]];
1418 sfree(native_lambda_vec);
1419 sfree(lambda_components);
1420 if (nblock_hist == 0 && nblock_dh == 0)
1422 /* don't do anything */
1425 if (nblock_hist > 0 && nblock_dh > 0)
1428 "This energy file contains both histogram dhdl data and non-histogram dhdl data. "
1429 "Don't know what to do.");
1435 /* we have standard, non-histogram data --
1436 call open_dhdl to open the file */
1437 /* TODO this is an ugly hack that needs to be fixed: this will only
1438 work if the order of data is always the same and if we're
1439 only using the g_energy compiled with the mdrun that produced
1441 *fp_dhdl = open_dhdl(filename, ir, oenv);
1445 sprintf(title, "N(%s)", deltag);
1446 sprintf(label_x, "%s (%s)", deltag, unit_energy);
1447 sprintf(label_y, "Samples");
1448 *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1449 sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
1450 xvgr_subtitle(*fp_dhdl, buf, oenv);
1454 (*hists) += nblock_hist;
1455 (*blocks) += nblock_dh;
1456 (*nlambdas) = nblock_hist + nblock_dh;
1458 /* write the data */
1459 if (nblock_hist > 0)
1463 for (i = 0; i < fr->nblock; i++)
1465 t_enxblock* blk = &(fr->block[i]);
1466 if (blk->id == enxDHHIST)
1468 double foreign_lambda, dx;
1470 int nhist, derivative;
1472 /* check the block types etc. */
1473 if ((blk->nsub < 2) || (blk->sub[0].type != xdr_datatype_double)
1474 || (blk->sub[1].type != xdr_datatype_int64) || (blk->sub[0].nr < 2)
1475 || (blk->sub[1].nr < 2))
1477 gmx_fatal(FARGS, "Unexpected block data in file");
1479 foreign_lambda = blk->sub[0].dval[0];
1480 dx = blk->sub[0].dval[1];
1481 nhist = blk->sub[1].lval[0];
1482 derivative = blk->sub[1].lval[1];
1483 for (j = 0; j < nhist; j++)
1486 x0 = blk->sub[1].lval[2 + j];
1490 sprintf(legend, "N(%s(%s=%g) | %s=%g)", deltag, lambda, foreign_lambda,
1491 lambda, start_lambda);
1495 sprintf(legend, "N(%s | %s=%g)", dhdl, lambda, start_lambda);
1499 xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1501 for (k = 0; k < blk->sub[j + 2].nr; k++)
1506 hist = blk->sub[j + 2].ival[k];
1507 xmin = (x0 + k) * dx;
1508 xmax = (x0 + k + 1) * dx;
1509 fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist, xmax, hist);
1512 /* multiple histogram data blocks in one histogram
1513 mean that the second one is the reverse of the first one:
1514 for dhdl derivatives, it's important to know both the
1515 maximum and minimum values */
1520 (*samples) += static_cast<int>(sum / nblock_hist);
1527 for (i = 0; i < fr->nblock; i++)
1529 t_enxblock* blk = &(fr->block[i]);
1530 if (blk->id == enxDH)
1534 len = blk->sub[2].nr;
1538 if (len != blk->sub[2].nr)
1540 gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1547 for (i = 0; i < len; i++)
1549 double time = start_time + delta_time * i;
1551 fprintf(*fp_dhdl, "%.4f ", time);
1553 for (j = 0; j < fr->nblock; j++)
1555 t_enxblock* blk = &(fr->block[j]);
1556 if (blk->id == enxDH)
1559 if (blk->sub[2].type == xdr_datatype_float)
1561 value = blk->sub[2].fval[i];
1565 value = blk->sub[2].dval[i];
1567 /* we need to decide which data type it is based on the count*/
1569 if (j == 1 && ir->bExpanded)
1571 fprintf(*fp_dhdl, "%4d",
1572 static_cast<int>(value)); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */
1578 fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */
1582 fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */
1587 fprintf(*fp_dhdl, "\n");
1593 int gmx_energy(int argc, char* argv[])
1595 const char* desc[] = {
1596 "[THISMODULE] extracts energy components",
1597 "from an energy file. The user is prompted to interactively",
1598 "select the desired energy terms.[PAR]",
1600 "Average, RMSD, and drift are calculated with full precision from the",
1601 "simulation (see printed manual). Drift is calculated by performing",
1602 "a least-squares fit of the data to a straight line. The reported total drift",
1603 "is the difference of the fit at the first and last point.",
1604 "An error estimate of the average is given based on a block averages",
1605 "over 5 blocks using the full-precision averages. The error estimate",
1606 "can be performed over multiple block lengths with the options",
1607 "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1608 "[BB]Note[bb] that in most cases the energy files contains averages over all",
1609 "MD steps, or over many more points than the number of frames in",
1610 "energy file. This makes the [THISMODULE] statistics output more accurate",
1611 "than the [REF].xvg[ref] output. When exact averages are not present in the energy",
1612 "file, the statistics mentioned above are simply over the single, per-frame",
1613 "energy values.[PAR]",
1615 "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1617 "Some fluctuation-dependent properties can be calculated provided",
1618 "the correct energy terms are selected, and that the command line option",
1619 "[TT]-fluct_props[tt] is given. The following properties",
1620 "will be computed:",
1622 "=============================== ===================",
1623 "Property Energy terms needed",
1624 "=============================== ===================",
1625 "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp",
1626 "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp",
1627 "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp",
1628 "Isothermal compressibility: Vol, Temp",
1629 "Adiabatic bulk modulus: Vol, Temp",
1630 "=============================== ===================",
1632 "You always need to set the number of molecules [TT]-nmol[tt].",
1633 "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1634 "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]",
1636 "Option [TT]-odh[tt] extracts and plots the free energy data",
1637 "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1638 "from the [TT]ener.edr[tt] file.[PAR]",
1640 "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1641 "difference with an ideal gas state::",
1643 " [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT ",
1644 " [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln]",
1645 " [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT ",
1646 " [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln]",
1648 "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1649 "the average is over the ensemble (or time in a trajectory).",
1650 "Note that this is in principle",
1651 "only correct when averaging over the whole (Boltzmann) ensemble",
1652 "and using the potential energy. This also allows for an entropy",
1655 " [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ",
1656 " ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T",
1657 " [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ",
1658 " ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T",
1661 "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1662 "difference is calculated::",
1665 " [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub]) / ",
1666 " kT[exp][chevron][SUB]A[sub][ln],",
1668 "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1669 "files, and the average is over the ensemble A. The running average",
1670 "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1671 "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1674 static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1675 static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1676 static int nmol = 1, nbmin = 5, nbmax = 5;
1677 static real reftemp = 300.0, ezero = 0;
1679 { "-fee", FALSE, etBOOL, { &bFee }, "Do a free energy estimate" },
1684 "Reference temperature for free energy calculation" },
1685 { "-zero", FALSE, etREAL, { &ezero }, "Subtract a zero-point energy" },
1690 "Sum the energy terms selected rather than display them all" },
1691 { "-dp", FALSE, etBOOL, { &bDp }, "Print energies in high precision" },
1692 { "-nbmin", FALSE, etINT, { &nbmin }, "Minimum number of blocks for error estimate" },
1693 { "-nbmax", FALSE, etINT, { &nbmax }, "Maximum number of blocks for error estimate" },
1698 "Compute the total dipole moment from the components" },
1703 "Also print the exact average and rmsd stored in the energy frames (only when 1 term is "
1709 "Number of molecules in your sample: the energies are divided by this number" },
1714 "Compute properties based on energy fluctuations, like heat capacity" },
1719 "Useful only for calculations of fluctuation properties. The drift in the observables "
1720 "will be subtracted before computing the fluctuation properties." },
1725 "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1726 { "-orinst", FALSE, etBOOL, { &bOrinst }, "Analyse instantaneous orientation data" },
1727 { "-ovec", FALSE, etBOOL, { &bOvec }, "Also plot the eigenvectors with [TT]-oten[tt]" }
1729 static const char* setnm[] = { "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX",
1730 "Pres-YY", "Pres-YZ", "Pres-ZX", "Pres-ZY",
1731 "Pres-ZZ", "Temperature", "Volume", "Pressure" };
1733 FILE* out = nullptr;
1734 FILE* fp_dhdl = nullptr;
1738 gmx_enxnm_t* enm = nullptr;
1739 t_enxframe * frame, *fr = nullptr;
1741 #define NEXT (1 - cur)
1746 gmx_bool bFoundStart, bCont, bVisco;
1748 double* time = nullptr;
1750 int * set = nullptr, i, j, nset, sss;
1751 gmx_bool* bIsEner = nullptr;
1752 char** leg = nullptr;
1754 gmx_output_env_t* oenv;
1755 int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
1758 { efEDR, "-f", nullptr, ffREAD }, { efEDR, "-f2", nullptr, ffOPTRD },
1759 { efTPR, "-s", nullptr, ffOPTRD }, { efXVG, "-o", "energy", ffWRITE },
1760 { efXVG, "-viol", "violaver", ffOPTWR }, { efXVG, "-pairs", "pairs", ffOPTWR },
1761 { efXVG, "-corr", "enecorr", ffOPTWR }, { efXVG, "-vis", "visco", ffOPTWR },
1762 { efXVG, "-evisco", "evisco", ffOPTWR }, { efXVG, "-eviscoi", "eviscoi", ffOPTWR },
1763 { efXVG, "-ravg", "runavgdf", ffOPTWR }, { efXVG, "-odh", "dhdl", ffOPTWR }
1765 #define NFILE asize(fnm)
1770 ppa = add_acf_pargs(&npargs, pa);
1771 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END, NFILE, fnm,
1772 npargs, ppa, asize(desc), desc, 0, nullptr, &oenv))
1778 bDHDL = opt2bSet("-odh", NFILE, fnm);
1783 fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
1784 do_enxnms(fp, &nre, &enm);
1788 bVisco = opt2bSet("-vis", NFILE, fnm);
1790 t_inputrec irInstance;
1791 t_inputrec* ir = &irInstance;
1797 nset = asize(setnm);
1799 /* This is nasty code... To extract Pres tensor, Volume and Temperature */
1800 for (j = 0; j < nset; j++)
1802 for (i = 0; i < nre; i++)
1804 if (std::strstr(enm[i].name, setnm[j]))
1812 if (gmx_strcasecmp(setnm[j], "Volume") == 0)
1814 printf("Enter the box volume (" unit_volume "): ");
1815 if (1 != scanf("%lf", &dbl))
1817 gmx_fatal(FARGS, "Error reading user input");
1823 gmx_fatal(FARGS, "Could not find term %s for viscosity calculation", setnm[j]);
1830 set = select_by_name(nre, enm, &nset);
1832 /* Print all the different units once */
1833 sprintf(buf, "(%s)", enm[set[0]].unit);
1834 for (i = 1; i < nset; i++)
1836 for (j = 0; j < i; j++)
1838 if (std::strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
1845 std::strcat(buf, ", (");
1846 std::strcat(buf, enm[set[i]].unit);
1847 std::strcat(buf, ")");
1850 out = xvgropen(opt2fn("-o", NFILE, fnm), "GROMACS Energies", "Time (ps)", buf, oenv);
1852 snew(leg, nset + 1);
1853 for (i = 0; (i < nset); i++)
1855 leg[i] = enm[set[i]].name;
1859 leg[nset] = gmx_strdup("Sum");
1860 xvgr_legend(out, nset + 1, leg, oenv);
1864 xvgr_legend(out, nset, leg, oenv);
1867 snew(bIsEner, nset);
1868 for (i = 0; (i < nset); i++)
1871 for (j = 0; (j <= F_ETOT); j++)
1873 bIsEner[i] = bIsEner[i]
1874 || (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
1877 if (bPrAll && nset > 1)
1879 gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
1884 get_dhdl_parms(ftp2fn(efTPR, NFILE, fnm), ir);
1887 /* Initiate energies and set them to zero */
1891 edat.step = nullptr;
1892 edat.steps = nullptr;
1893 edat.points = nullptr;
1894 edat.bHaveSums = TRUE;
1897 /* Initiate counters */
1898 bFoundStart = FALSE;
1903 /* This loop searches for the first frame (when -b option is given),
1904 * or when this has been found it reads just one energy frame
1908 bCont = do_enx(fp, &(frame[NEXT]));
1911 timecheck = check_times(frame[NEXT].t);
1913 } while (bCont && (timecheck < 0));
1915 if ((timecheck == 0) && bCont)
1917 /* We read a valid frame, so we can use it */
1918 fr = &(frame[NEXT]);
1922 /* The frame contains energies, so update cur */
1925 if (edat.nframes % 1000 == 0)
1927 srenew(edat.step, edat.nframes + 1000);
1928 std::memset(&(edat.step[edat.nframes]), 0, 1000 * sizeof(edat.step[0]));
1929 srenew(edat.steps, edat.nframes + 1000);
1930 std::memset(&(edat.steps[edat.nframes]), 0, 1000 * sizeof(edat.steps[0]));
1931 srenew(edat.points, edat.nframes + 1000);
1932 std::memset(&(edat.points[edat.nframes]), 0, 1000 * sizeof(edat.points[0]));
1934 for (i = 0; i < nset; i++)
1936 srenew(edat.s[i].ener, edat.nframes + 1000);
1937 std::memset(&(edat.s[i].ener[edat.nframes]), 0, 1000 * sizeof(edat.s[i].ener[0]));
1938 srenew(edat.s[i].es, edat.nframes + 1000);
1939 std::memset(&(edat.s[i].es[edat.nframes]), 0, 1000 * sizeof(edat.s[i].es[0]));
1944 edat.step[nfr] = fr->step;
1949 /* Initiate the previous step data */
1950 start_step = fr->step;
1952 /* Initiate the energy sums */
1953 edat.steps[nfr] = 1;
1954 edat.points[nfr] = 1;
1955 for (i = 0; i < nset; i++)
1958 edat.s[i].es[nfr].sum = fr->ener[sss].e;
1959 edat.s[i].es[nfr].sum2 = 0;
1966 edat.steps[nfr] = fr->nsteps;
1970 /* mdrun only calculated the energy at energy output
1971 * steps. We don't need to check step intervals.
1973 edat.points[nfr] = 1;
1974 for (i = 0; i < nset; i++)
1977 edat.s[i].es[nfr].sum = fr->ener[sss].e;
1978 edat.s[i].es[nfr].sum2 = 0;
1981 edat.bHaveSums = FALSE;
1983 else if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
1985 /* We have statistics to the previous frame */
1986 edat.points[nfr] = fr->nsum;
1987 for (i = 0; i < nset; i++)
1990 edat.s[i].es[nfr].sum = fr->ener[sss].esum;
1991 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
1993 edat.npoints += fr->nsum;
1997 /* The interval does not match fr->nsteps:
1998 * can not do exact averages.
2000 edat.bHaveSums = FALSE;
2003 edat.nsteps = fr->step - start_step + 1;
2005 for (i = 0; i < nset; i++)
2007 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2011 * Store energies for analysis afterwards...
2013 if (!bDHDL && (fr->nre > 0))
2015 if (edat.nframes % 1000 == 0)
2017 srenew(time, edat.nframes + 1000);
2019 time[edat.nframes] = fr->t;
2024 do_dhdl(fr, ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists,
2025 &dh_samples, &dh_lambdas, oenv);
2028 /*******************************************
2030 *******************************************/
2037 /* We skip frames with single points (usually only the first frame),
2038 * since they would result in an average plot with outliers.
2042 print_time(out, fr->t);
2043 print1(out, bDp, fr->ener[set[0]].e);
2044 print1(out, bDp, fr->ener[set[0]].esum / fr->nsum);
2045 print1(out, bDp, std::sqrt(fr->ener[set[0]].eav / fr->nsum));
2051 print_time(out, fr->t);
2055 for (i = 0; i < nset; i++)
2057 sum += fr->ener[set[i]].e;
2059 print1(out, bDp, sum / nmol - ezero);
2063 for (i = 0; (i < nset); i++)
2067 print1(out, bDp, (fr->ener[set[i]].e) / nmol - ezero);
2071 print1(out, bDp, fr->ener[set[i]].e);
2080 } while (bCont && (timecheck == 0));
2082 fprintf(stderr, "\n");
2093 gmx_fio_fclose(fp_dhdl);
2094 printf("\n\nWrote %d lambda values with %d samples as ", dh_lambdas, dh_samples);
2097 printf("%d dH histograms ", dh_hists);
2101 printf("%d dH data blocks ", dh_blocks);
2103 printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2107 gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2112 double dt = (frame[cur].t - start_t) / (edat.nframes - 1);
2113 analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
2114 opt2fn("-evisco", NFILE, fnm), opt2fn("-eviscoi", NFILE, fnm), bFee, bSum,
2115 bFluct, bVisco, opt2fn("-vis", NFILE, fnm), nmol, start_step, start_t,
2116 frame[cur].step, frame[cur].t, reftemp, &edat, nset, set, bIsEner, leg, enm,
2117 Vaver, ezero, nbmin, nbmax, oenv);
2120 calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat, nbmin, nbmax);
2123 if (opt2bSet("-f2", NFILE, fnm))
2125 fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm), reftemp, nset, set, leg, &edat,
2129 done_enerdata_t(nset, &edat);
2131 free_enxframe(&frame[0]);
2132 free_enxframe(&frame[1]);
2134 free_enxnms(nre, enm);
2140 const char* nxy = "-nxy";
2142 do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2143 do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2144 do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);
2146 output_env_done(oenv);