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,2021, 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);
344 fn, "Shear viscosity using Einstein relation", "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
345 for (i = 0; i < nf4; i++)
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];
600 "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
603 gmx_step_str(eee[nb].nst_min, buf1),
604 gmx_step_str(edat->nsteps, buf2));
606 if (eee[nb].b == nb && 5 * nb * eee[nb].nst_min >= 4 * edat->nsteps)
608 see2 += calc_ee2(nb, &eee[nb].sum);
614 edat->s[i].ee = std::sqrt(see2 / nee);
624 static enerdata_t* calc_sum(int nset, enerdata_t* edat, int nbmin, int nbmax)
635 snew(s->ener, esum->nframes);
636 snew(s->es, esum->nframes);
638 s->bExactStat = TRUE;
640 for (i = 0; i < nset; i++)
642 if (!edat->s[i].bExactStat)
644 s->bExactStat = FALSE;
646 s->slope += edat->s[i].slope;
649 for (f = 0; f < edat->nframes; f++)
652 for (i = 0; i < nset; i++)
654 sum += edat->s[i].ener[f];
658 for (i = 0; i < nset; i++)
660 sum += edat->s[i].es[f].sum;
666 calc_averages(1, esum, nbmin, nbmax);
671 static void ee_pr(double ee, int buflen, char* buf)
673 snprintf(buf, buflen, "%s", "--");
676 /* Round to two decimals by printing. */
678 snprintf(tmp, sizeof(tmp), "%.1e", ee);
679 double rnd = gmx::doubleFromString(tmp);
680 snprintf(buf, buflen, "%g", rnd);
684 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t* edat)
686 /* Remove the drift by performing a fit to y = ax+b.
687 Uses 5 iterations. */
691 edat->npoints = edat->nframes;
692 edat->nsteps = edat->nframes;
694 for (k = 0; (k < 5); k++)
696 for (i = 0; (i < nset); i++)
698 delta = edat->s[i].slope * dt;
700 if (nullptr != debug)
702 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
705 for (j = 0; (j < edat->nframes); j++)
707 edat->s[i].ener[j] -= j * delta;
708 edat->s[i].es[j].sum = 0;
709 edat->s[i].es[j].sum2 = 0;
712 calc_averages(nset, edat, nbmin, nbmax);
716 static void calc_fluctuation_props(FILE* fp,
727 double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, varet;
737 const char* my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
740 NANO3 = NANO * NANO * NANO;
744 "\nYou may want to use the -driftcorr flag in order to correct\n"
745 "for spurious drift in the graphs. Note that this is not\n"
746 "a substitute for proper equilibration and sampling!\n");
750 remove_drift(nset, nbmin, nbmax, dt, edat);
752 for (i = 0; (i < eNR); i++)
754 for (ii[i] = 0; (ii[i] < nset && (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++) {}
756 fprintf(fp,"Found %s data.\n",my_ener[i]);
758 /* Compute it all! */
759 alpha = kappa = cp = dcp = cv = NOTSET;
763 if (ii[eTemp] < nset)
765 tt = edat->s[ii[eTemp]].av;
769 if ((ii[eVol] < nset) && (ii[eTemp] < nset))
771 vv = edat->s[ii[eVol]].av * NANO3;
772 varv = gmx::square(edat->s[ii[eVol]].rmsd * NANO3);
773 kappa = (varv / vv) / (BOLTZMANN * tt);
777 if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
779 hh = KILO * edat->s[ii[eEnth]].av / AVOGADRO;
780 varh = gmx::square(KILO * edat->s[ii[eEnth]].rmsd / AVOGADRO);
781 cp = AVOGADRO * ((varh / nmol) / (BOLTZMANN * tt * tt));
784 if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
786 /* Only compute cv in constant volume runs, which we can test
787 by checking whether the enthalpy was computed.
789 varet = gmx::square(edat->s[ii[eEtot]].rmsd);
790 cv = KILO * ((varet / nmol) / (BOLTZ * tt * tt));
793 if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
795 double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
796 vh_sum = v_sum = h_sum = 0;
797 for (j = 0; (j < edat->nframes); j++)
799 v = edat->s[ii[eVol]].ener[j] * NANO3;
800 h = KILO * edat->s[ii[eEnth]].ener[j] / AVOGADRO;
805 vh_aver = vh_sum / edat->nframes;
806 v_aver = v_sum / edat->nframes;
807 h_aver = h_sum / edat->nframes;
808 alpha = (vh_aver - v_aver * h_aver) / (v_aver * BOLTZMANN * tt * tt);
809 dcp = (v_aver * AVOGADRO / nmol) * tt * gmx::square(alpha) / (kappa);
816 fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n", nmol);
818 fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
819 fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
820 fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
821 fprintf(fp, "please use the g_dos program.\n\n");
823 "WARNING: Please verify that your simulations are converged and perform\n"
824 "a block-averaging error analysis (not implemented in g_energy yet)\n");
826 if (debug != nullptr)
830 fprintf(fp, "varv = %10g (m^6)\n", varv * AVOGADRO / nmol);
835 fprintf(fp, "Volume = %10g m^3/mol\n", vv * AVOGADRO / nmol);
839 fprintf(fp, "Enthalpy = %10g kJ/mol\n", hh * AVOGADRO / (KILO * nmol));
843 fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n", alpha);
847 fprintf(fp, "Isothermal Compressibility Kappa = %10g (m^3/J)\n", kappa);
848 fprintf(fp, "Adiabatic bulk modulus = %10g (J/m^3)\n", 1.0 / kappa);
852 fprintf(fp, "Heat capacity at constant pressure Cp = %10g J/(mol K)\n", cp);
856 fprintf(fp, "Heat capacity at constant volume Cv = %10g J/(mol K)\n", cv);
860 fprintf(fp, "Cp-Cv = %10g J/(mol K)\n", dcp);
862 please_cite(fp, "Allen1987a");
866 fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
870 static void analyse_ener(gmx_bool bCorr,
872 const char* eviscofn,
873 const char* eviscoifn,
888 const gmx_bool* bIsEner,
895 const gmx_output_env_t* oenv)
898 /* Check out the printed manual for equations! */
899 double Dt, aver, stddev, errest, delta_t, totaldrift;
900 enerdata_t* esum = nullptr;
901 real integral, intBulk, Temp = 0, Pres = 0;
902 real pr_aver, pr_stddev, pr_errest;
903 double beta = 0, expE, expEtot, *fee = nullptr;
905 int nexact, nnotexact;
907 char buf[256], eebuf[100];
909 nsteps = step - start_step + 1;
912 fprintf(stdout, "Not enough steps (%s) for statistics\n", gmx_step_str(nsteps, buf));
916 /* Calculate the time difference */
917 delta_t = t - start_t;
920 "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
921 gmx_step_str(nsteps, buf),
926 calc_averages(nset, edat, nbmin, nbmax);
930 esum = calc_sum(nset, edat, nbmin, nbmax);
933 if (!edat->bHaveSums)
942 for (i = 0; (i < nset); i++)
944 if (edat->s[i].bExactStat)
957 fprintf(stdout, "All statistics are over %s points\n", gmx_step_str(edat->npoints, buf));
959 else if (nexact == 0 || edat->npoints == edat->nframes)
961 fprintf(stdout, "All statistics are over %d points (frames)\n", edat->nframes);
965 fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
966 for (i = 0; (i < nset); i++)
968 if (!edat->s[i].bExactStat)
970 fprintf(stdout, " '%s'", leg[i]);
974 " %s has statistics over %d points (frames)\n",
975 nnotexact == 1 ? "is" : "are",
977 fprintf(stdout, "All other statistics are over %s points\n", gmx_step_str(edat->npoints, buf));
979 fprintf(stdout, "\n");
982 "%-24s %10s %10s %10s %10s",
990 fprintf(stdout, " %10s\n", "-kT ln<e^(E/kT)>");
994 fprintf(stdout, "\n");
997 "-------------------------------------------------------------------------------"
1000 /* Initiate locals, only used with -sum */
1004 beta = 1.0 / (BOLTZ * reftemp);
1007 for (i = 0; (i < nset); i++)
1009 aver = edat->s[i].av;
1010 stddev = edat->s[i].rmsd;
1011 errest = edat->s[i].ee;
1016 for (j = 0; (j < edat->nframes); j++)
1018 expE += std::exp(beta * (edat->s[i].ener[j] - aver) / nmol);
1022 expEtot += expE / edat->nframes;
1025 fee[i] = std::log(expE / edat->nframes) / beta + aver / nmol;
1027 if (std::strstr(leg[i], "empera") != nullptr)
1031 else if (std::strstr(leg[i], "olum") != nullptr)
1035 else if (std::strstr(leg[i], "essure") != nullptr)
1041 pr_aver = aver / nmol - ezero;
1042 pr_stddev = stddev / nmol;
1043 pr_errest = errest / nmol;
1052 /* Multiply the slope in steps with the number of steps taken */
1053 totaldrift = (edat->nsteps - 1) * edat->s[i].slope;
1059 ee_pr(pr_errest, sizeof(eebuf), eebuf);
1060 fprintf(stdout, "%-24s %10g %10s %10g %10g", leg[i], pr_aver, eebuf, pr_stddev, totaldrift);
1063 fprintf(stdout, " %10g", fee[i]);
1066 fprintf(stdout, " (%s)\n", enm[set[i]].unit);
1070 for (j = 0; (j < edat->nframes); j++)
1072 edat->s[i].ener[j] -= aver;
1078 totaldrift = (edat->nsteps - 1) * esum->s[0].slope;
1079 ee_pr(esum->s[0].ee / nmol, sizeof(eebuf), eebuf);
1081 "%-24s %10g %10s %10s %10g (%s)",
1083 esum->s[0].av / nmol,
1088 /* pr_aver,pr_stddev,a,totaldrift */
1093 std::log(expEtot) / beta + esum->s[0].av / nmol,
1094 std::log(expEtot) / beta);
1098 fprintf(stdout, "\n");
1102 /* Do correlation function */
1103 if (edat->nframes > 1)
1105 Dt = delta_t / (edat->nframes - 1);
1113 const char* leg[] = { "Shear", "Bulk" };
1118 /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1120 /* Symmetrise tensor! (and store in first three elements)
1121 * And subtract average pressure!
1124 for (i = 0; i < 12; i++)
1126 snew(eneset[i], edat->nframes);
1128 for (i = 0; (i < edat->nframes); i++)
1130 eneset[0][i] = 0.5 * (edat->s[1].ener[i] + edat->s[3].ener[i]);
1131 eneset[1][i] = 0.5 * (edat->s[2].ener[i] + edat->s[6].ener[i]);
1132 eneset[2][i] = 0.5 * (edat->s[5].ener[i] + edat->s[7].ener[i]);
1133 for (j = 3; j <= 11; j++)
1135 eneset[j][i] = edat->s[j].ener[i];
1137 eneset[11][i] -= Pres;
1140 /* Determine integrals of the off-diagonal pressure elements */
1142 for (i = 0; i < 3; i++)
1144 snew(eneint[i], edat->nframes + 1);
1149 for (i = 0; i < edat->nframes; i++)
1153 + 0.5 * (edat->s[1].es[i].sum + edat->s[3].es[i].sum) * Dt / edat->points[i];
1156 + 0.5 * (edat->s[2].es[i].sum + edat->s[6].es[i].sum) * Dt / edat->points[i];
1159 + 0.5 * (edat->s[5].es[i].sum + edat->s[7].es[i].sum) * Dt / edat->points[i];
1162 einstein_visco(eviscofn, eviscoifn, 3, edat->nframes + 1, eneint, Vaver, Temp, Dt, oenv);
1164 for (i = 0; i < 3; i++)
1170 /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1171 /* Do it for shear viscosity */
1172 std::strcpy(buf, "Shear Viscosity");
1173 low_do_autocorr(corrfn,
1178 (edat->nframes + 1) / 2,
1190 /* Now for bulk viscosity */
1191 std::strcpy(buf, "Bulk Viscosity");
1192 low_do_autocorr(corrfn,
1197 (edat->nframes + 1) / 2,
1209 factor = (Vaver * 1e-26 / (BOLTZMANN * Temp)) * Dt;
1210 fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1211 xvgr_legend(fp, asize(leg), leg, oenv);
1213 /* Use trapezium rule for integration */
1216 nout = get_acfnout();
1217 if ((nout < 2) || (nout >= edat->nframes / 2))
1219 nout = edat->nframes / 2;
1221 for (i = 1; (i < nout); i++)
1223 integral += 0.5 * (eneset[0][i - 1] + eneset[0][i]) * factor;
1224 intBulk += 0.5 * (eneset[11][i - 1] + eneset[11][i]) * factor;
1225 fprintf(fp, "%10g %10g %10g\n", (i * Dt), integral, intBulk);
1229 for (i = 0; i < 12; i++)
1239 std::strcpy(buf, "Autocorrelation of Energy Fluctuations");
1243 std::strcpy(buf, "Energy Autocorrelation");
1246 do_autocorr(corrfn, oenv, buf, edat->nframes,
1248 bSum ? &edat->s[nset-1].ener : eneset,
1249 (delta_t/edat->nframes), eacNormal, FALSE);
1255 static void print_time(FILE* fp, double t)
1257 fprintf(fp, "%12.6f", t);
1260 static void print1(FILE* fp, gmx_bool bDp, real e)
1264 fprintf(fp, " %16.12f", e);
1268 fprintf(fp, " %10.6f", e);
1272 static void fec(const char* ene2fn,
1273 const char* runavgfn,
1280 const gmx_output_env_t* oenv)
1282 const char* ravgleg[] = { "\\8D\\4E = E\\sB\\N-E\\sA\\N", "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N" };
1285 int timecheck, nenergy, nenergy2, maxenergy;
1291 gmx_enxnm_t* enm = nullptr;
1295 /* read second energy file */
1298 enx = open_enx(ene2fn, "r");
1299 do_enxnms(enx, &(fr->nre), &enm);
1301 snew(eneset2, nset + 1);
1307 /* This loop searches for the first frame (when -b option is given),
1308 * or when this has been found it reads just one energy frame
1312 bCont = do_enx(enx, fr);
1316 timecheck = check_times(fr->t);
1319 } while (bCont && (timecheck < 0));
1321 /* Store energies for analysis afterwards... */
1322 if ((timecheck == 0) && bCont)
1326 if (nenergy2 >= maxenergy)
1329 for (i = 0; i <= nset; i++)
1331 srenew(eneset2[i], maxenergy);
1334 GMX_RELEASE_ASSERT(time != nullptr, "trying to dereference NULL time pointer");
1336 if (fr->t != time[nenergy2])
1339 "\nWARNING time mismatch %g!=%g at frame %s\n",
1342 gmx_step_str(fr->step, buf));
1344 for (i = 0; i < nset; i++)
1346 eneset2[i][nenergy2] = fr->ener[set[i]].e;
1351 } while (bCont && (timecheck == 0));
1354 if (edat->nframes != nenergy2)
1356 fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n", edat->nframes, nenergy2);
1358 nenergy = std::min(edat->nframes, nenergy2);
1360 /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1364 fp = xvgropen(runavgfn,
1365 "Running average free energy difference",
1366 "Time (" unit_time ")",
1367 "\\8D\\4E (" unit_energy ")",
1369 xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
1371 fprintf(stdout, "\n%-24s %10s\n", "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
1373 beta = 1.0 / (BOLTZ * reftemp);
1374 for (i = 0; i < nset; i++)
1376 if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
1378 fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n", leg[i], enm[set[i]].name);
1380 for (j = 0; j < nenergy; j++)
1382 dE = eneset2[i][j] - edat->s[i].ener[j];
1383 sum += std::exp(-dE * beta);
1386 fprintf(fp, "%10g %10g %10g\n", time[j], dE, -BOLTZ * reftemp * std::log(sum / (j + 1)));
1389 aver = -BOLTZ * reftemp * std::log(sum / nenergy);
1390 fprintf(stdout, "%-24s %10g\n", leg[i], aver);
1400 static void do_dhdl(t_enxframe* fr,
1401 const t_inputrec* ir,
1403 const char* filename,
1409 const gmx_output_env_t* oenv)
1411 const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1412 char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
1414 int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
1417 double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0;
1418 static int setnr = 0;
1419 double* native_lambda_vec = nullptr;
1420 const char** lambda_components = nullptr;
1421 int n_lambda_vec = 0;
1422 bool firstPass = true;
1424 /* now count the blocks & handle the global dh data */
1425 for (i = 0; i < fr->nblock; i++)
1427 if (fr->block[i].id == enxDHHIST)
1431 else if (fr->block[i].id == enxDH)
1435 else if (fr->block[i].id == enxDHCOLL)
1438 if ((fr->block[i].nsub < 1) || (fr->block[i].sub[0].type != xdr_datatype_double)
1439 || (fr->block[i].sub[0].nr < 5))
1441 gmx_fatal(FARGS, "Unexpected block data");
1444 /* read the data from the DHCOLL block */
1445 temp = fr->block[i].sub[0].dval[0];
1446 start_time = fr->block[i].sub[0].dval[1];
1447 delta_time = fr->block[i].sub[0].dval[2];
1448 start_lambda = fr->block[i].sub[0].dval[3];
1449 if (fr->block[i].nsub > 1)
1453 n_lambda_vec = fr->block[i].sub[1].ival[1];
1454 snew(lambda_components, n_lambda_vec);
1455 snew(native_lambda_vec, n_lambda_vec);
1460 if (n_lambda_vec != fr->block[i].sub[1].ival[1])
1462 gmx_fatal(FARGS, "Unexpected change of basis set in lambda");
1465 for (j = 0; j < n_lambda_vec; j++)
1467 native_lambda_vec[j] = fr->block[i].sub[0].dval[5 + j];
1468 lambda_components[j] =
1469 enumValueToStringSingular(static_cast<FreeEnergyPerturbationCouplingType>(
1470 fr->block[i].sub[1].ival[2 + j]));
1476 sfree(native_lambda_vec);
1477 sfree(lambda_components);
1478 if (nblock_hist == 0 && nblock_dh == 0)
1480 /* don't do anything */
1483 if (nblock_hist > 0 && nblock_dh > 0)
1486 "This energy file contains both histogram dhdl data and non-histogram dhdl data. "
1487 "Don't know what to do.");
1493 /* we have standard, non-histogram data --
1494 call open_dhdl to open the file */
1495 /* TODO this is an ugly hack that needs to be fixed: this will only
1496 work if the order of data is always the same and if we're
1497 only using the g_energy compiled with the mdrun that produced
1499 *fp_dhdl = open_dhdl(filename, ir, oenv);
1503 sprintf(title, "N(%s)", deltag);
1504 sprintf(label_x, "%s (%s)", deltag, unit_energy);
1505 sprintf(label_y, "Samples");
1506 *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1507 sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
1508 xvgr_subtitle(*fp_dhdl, buf, oenv);
1512 (*hists) += nblock_hist;
1513 (*blocks) += nblock_dh;
1514 (*nlambdas) = nblock_hist + nblock_dh;
1516 /* write the data */
1517 if (nblock_hist > 0)
1521 for (i = 0; i < fr->nblock; i++)
1523 t_enxblock* blk = &(fr->block[i]);
1524 if (blk->id == enxDHHIST)
1526 double foreign_lambda, dx;
1528 int nhist, derivative;
1530 /* check the block types etc. */
1531 if ((blk->nsub < 2) || (blk->sub[0].type != xdr_datatype_double)
1532 || (blk->sub[1].type != xdr_datatype_int64) || (blk->sub[0].nr < 2)
1533 || (blk->sub[1].nr < 2))
1535 gmx_fatal(FARGS, "Unexpected block data in file");
1537 foreign_lambda = blk->sub[0].dval[0];
1538 dx = blk->sub[0].dval[1];
1539 nhist = blk->sub[1].lval[0];
1540 derivative = blk->sub[1].lval[1];
1541 for (j = 0; j < nhist; j++)
1544 x0 = blk->sub[1].lval[2 + j];
1548 sprintf(legend, "N(%s(%s=%g) | %s=%g)", deltag, lambda, foreign_lambda, lambda, start_lambda);
1552 sprintf(legend, "N(%s | %s=%g)", dhdl, lambda, start_lambda);
1556 xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1558 for (k = 0; k < blk->sub[j + 2].nr; k++)
1563 hist = blk->sub[j + 2].ival[k];
1564 xmin = (x0 + k) * dx;
1565 xmax = (x0 + k + 1) * dx;
1566 fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist, xmax, hist);
1569 /* multiple histogram data blocks in one histogram
1570 mean that the second one is the reverse of the first one:
1571 for dhdl derivatives, it's important to know both the
1572 maximum and minimum values */
1577 (*samples) += static_cast<int>(sum / nblock_hist);
1584 for (i = 0; i < fr->nblock; i++)
1586 t_enxblock* blk = &(fr->block[i]);
1587 if (blk->id == enxDH)
1591 len = blk->sub[2].nr;
1595 if (len != blk->sub[2].nr)
1597 gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1604 for (i = 0; i < len; i++)
1606 double time = start_time + delta_time * i;
1608 fprintf(*fp_dhdl, "%.4f ", time);
1610 for (j = 0; j < fr->nblock; j++)
1612 t_enxblock* blk = &(fr->block[j]);
1613 if (blk->id == enxDH)
1616 if (blk->sub[2].type == xdr_datatype_float)
1618 value = blk->sub[2].fval[i];
1622 value = blk->sub[2].dval[i];
1624 /* we need to decide which data type it is based on the count*/
1626 if (j == 1 && ir->bExpanded)
1628 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! */
1634 fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */
1638 fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */
1643 fprintf(*fp_dhdl, "\n");
1649 int gmx_energy(int argc, char* argv[])
1651 const char* desc[] = {
1652 "[THISMODULE] extracts energy components",
1653 "from an energy file. The user is prompted to interactively",
1654 "select the desired energy terms.[PAR]",
1656 "Average, RMSD, and drift are calculated with full precision from the",
1657 "simulation (see printed manual). Drift is calculated by performing",
1658 "a least-squares fit of the data to a straight line. The reported total drift",
1659 "is the difference of the fit at the first and last point.",
1660 "An error estimate of the average is given based on a block averages",
1661 "over 5 blocks using the full-precision averages. The error estimate",
1662 "can be performed over multiple block lengths with the options",
1663 "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1664 "[BB]Note[bb] that in most cases the energy files contains averages over all",
1665 "MD steps, or over many more points than the number of frames in",
1666 "energy file. This makes the [THISMODULE] statistics output more accurate",
1667 "than the [REF].xvg[ref] output. When exact averages are not present in the energy",
1668 "file, the statistics mentioned above are simply over the single, per-frame",
1669 "energy values.[PAR]",
1671 "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1673 "Some fluctuation-dependent properties can be calculated provided",
1674 "the correct energy terms are selected, and that the command line option",
1675 "[TT]-fluct_props[tt] is given. The following properties",
1676 "will be computed:",
1678 "=============================== ===================",
1679 "Property Energy terms needed",
1680 "=============================== ===================",
1681 "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp",
1682 "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp",
1683 "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp",
1684 "Isothermal compressibility: Vol, Temp",
1685 "Adiabatic bulk modulus: Vol, Temp",
1686 "=============================== ===================",
1688 "You always need to set the number of molecules [TT]-nmol[tt].",
1689 "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1690 "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]",
1692 "Option [TT]-odh[tt] extracts and plots the free energy data",
1693 "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1694 "from the [TT]ener.edr[tt] file.[PAR]",
1696 "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1697 "difference with an ideal gas state::",
1699 " [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT ",
1700 " [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln]",
1701 " [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT ",
1702 " [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln]",
1704 "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1705 "the average is over the ensemble (or time in a trajectory).",
1706 "Note that this is in principle",
1707 "only correct when averaging over the whole (Boltzmann) ensemble",
1708 "and using the potential energy. This also allows for an entropy",
1711 " [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ",
1712 " ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T",
1713 " [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ",
1714 " ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T",
1717 "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1718 "difference is calculated::",
1721 " [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub]) / ",
1722 " kT[exp][chevron][SUB]A[sub][ln],",
1724 "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1725 "files, and the average is over the ensemble A. The running average",
1726 "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1727 "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1730 static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1731 static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1732 static int nmol = 1, nbmin = 5, nbmax = 5;
1733 static real reftemp = 300.0, ezero = 0;
1735 { "-fee", FALSE, etBOOL, { &bFee }, "Do a free energy estimate" },
1740 "Reference temperature for free energy calculation" },
1741 { "-zero", FALSE, etREAL, { &ezero }, "Subtract a zero-point energy" },
1746 "Sum the energy terms selected rather than display them all" },
1747 { "-dp", FALSE, etBOOL, { &bDp }, "Print energies in high precision" },
1748 { "-nbmin", FALSE, etINT, { &nbmin }, "Minimum number of blocks for error estimate" },
1749 { "-nbmax", FALSE, etINT, { &nbmax }, "Maximum number of blocks for error estimate" },
1754 "Compute the total dipole moment from the components" },
1759 "Also print the exact average and rmsd stored in the energy frames (only when 1 term is "
1765 "Number of molecules in your sample: the energies are divided by this number" },
1770 "Compute properties based on energy fluctuations, like heat capacity" },
1775 "Useful only for calculations of fluctuation properties. The drift in the observables "
1776 "will be subtracted before computing the fluctuation properties." },
1781 "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1782 { "-orinst", FALSE, etBOOL, { &bOrinst }, "Analyse instantaneous orientation data" },
1783 { "-ovec", FALSE, etBOOL, { &bOvec }, "Also plot the eigenvectors with [TT]-oten[tt]" }
1785 static const char* setnm[] = { "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX",
1786 "Pres-YY", "Pres-YZ", "Pres-ZX", "Pres-ZY",
1787 "Pres-ZZ", "Temperature", "Volume", "Pressure" };
1789 FILE* out = nullptr;
1790 FILE* fp_dhdl = nullptr;
1794 gmx_enxnm_t* enm = nullptr;
1795 t_enxframe * frame, *fr = nullptr;
1797 #define NEXT (1 - cur)
1802 gmx_bool bFoundStart, bCont, bVisco;
1804 double* time = nullptr;
1806 int * set = nullptr, i, j, nset, sss;
1807 gmx_bool* bIsEner = nullptr;
1808 char** leg = nullptr;
1810 gmx_output_env_t* oenv;
1811 int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
1814 { efEDR, "-f", nullptr, ffREAD }, { efEDR, "-f2", nullptr, ffOPTRD },
1815 { efTPR, "-s", nullptr, ffOPTRD }, { efXVG, "-o", "energy", ffWRITE },
1816 { efXVG, "-viol", "violaver", ffOPTWR }, { efXVG, "-pairs", "pairs", ffOPTWR },
1817 { efXVG, "-corr", "enecorr", ffOPTWR }, { efXVG, "-vis", "visco", ffOPTWR },
1818 { efXVG, "-evisco", "evisco", ffOPTWR }, { efXVG, "-eviscoi", "eviscoi", ffOPTWR },
1819 { efXVG, "-ravg", "runavgdf", ffOPTWR }, { efXVG, "-odh", "dhdl", ffOPTWR }
1821 #define NFILE asize(fnm)
1826 ppa = add_acf_pargs(&npargs, pa);
1827 if (!parse_common_args(
1828 &argc, argv, PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END, NFILE, fnm, npargs, ppa, asize(desc), desc, 0, nullptr, &oenv))
1834 bDHDL = opt2bSet("-odh", NFILE, fnm);
1839 fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
1840 do_enxnms(fp, &nre, &enm);
1844 bVisco = opt2bSet("-vis", NFILE, fnm);
1846 t_inputrec irInstance;
1847 t_inputrec* ir = &irInstance;
1853 nset = asize(setnm);
1855 /* This is nasty code... To extract Pres tensor, Volume and Temperature */
1856 for (j = 0; j < nset; j++)
1858 for (i = 0; i < nre; i++)
1860 if (std::strstr(enm[i].name, setnm[j]))
1868 if (gmx_strcasecmp(setnm[j], "Volume") == 0)
1870 printf("Enter the box volume (" unit_volume "): ");
1871 if (1 != scanf("%lf", &dbl))
1873 gmx_fatal(FARGS, "Error reading user input");
1879 gmx_fatal(FARGS, "Could not find term %s for viscosity calculation", setnm[j]);
1886 set = select_by_name(nre, enm, &nset);
1888 /* Print all the different units once */
1889 sprintf(buf, "(%s)", enm[set[0]].unit);
1890 for (i = 1; i < nset; i++)
1892 for (j = 0; j < i; j++)
1894 if (std::strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
1901 std::strcat(buf, ", (");
1902 std::strcat(buf, enm[set[i]].unit);
1903 std::strcat(buf, ")");
1906 out = xvgropen(opt2fn("-o", NFILE, fnm), "GROMACS Energies", "Time (ps)", buf, oenv);
1908 snew(leg, nset + 1);
1909 for (i = 0; (i < nset); i++)
1911 leg[i] = enm[set[i]].name;
1915 leg[nset] = gmx_strdup("Sum");
1916 xvgr_legend(out, nset + 1, leg, oenv);
1920 xvgr_legend(out, nset, leg, oenv);
1923 snew(bIsEner, nset);
1924 for (i = 0; (i < nset); i++)
1927 for (j = 0; (j <= F_ETOT); j++)
1929 bIsEner[i] = bIsEner[i]
1930 || (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
1933 if (bPrAll && nset > 1)
1935 gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
1940 get_dhdl_parms(ftp2fn(efTPR, NFILE, fnm), ir);
1943 /* Initiate energies and set them to zero */
1947 edat.step = nullptr;
1948 edat.steps = nullptr;
1949 edat.points = nullptr;
1950 edat.bHaveSums = TRUE;
1953 /* Initiate counters */
1954 bFoundStart = FALSE;
1959 /* This loop searches for the first frame (when -b option is given),
1960 * or when this has been found it reads just one energy frame
1964 bCont = do_enx(fp, &(frame[NEXT]));
1967 timecheck = check_times(frame[NEXT].t);
1969 } while (bCont && (timecheck < 0));
1971 if ((timecheck == 0) && bCont)
1973 /* We read a valid frame, so we can use it */
1974 fr = &(frame[NEXT]);
1978 /* The frame contains energies, so update cur */
1981 if (edat.nframes % 1000 == 0)
1983 srenew(edat.step, edat.nframes + 1000);
1984 std::memset(&(edat.step[edat.nframes]), 0, 1000 * sizeof(edat.step[0]));
1985 srenew(edat.steps, edat.nframes + 1000);
1986 std::memset(&(edat.steps[edat.nframes]), 0, 1000 * sizeof(edat.steps[0]));
1987 srenew(edat.points, edat.nframes + 1000);
1988 std::memset(&(edat.points[edat.nframes]), 0, 1000 * sizeof(edat.points[0]));
1990 for (i = 0; i < nset; i++)
1992 srenew(edat.s[i].ener, edat.nframes + 1000);
1993 std::memset(&(edat.s[i].ener[edat.nframes]), 0, 1000 * sizeof(edat.s[i].ener[0]));
1994 srenew(edat.s[i].es, edat.nframes + 1000);
1995 std::memset(&(edat.s[i].es[edat.nframes]), 0, 1000 * sizeof(edat.s[i].es[0]));
2000 edat.step[nfr] = fr->step;
2005 /* Initiate the previous step data */
2006 start_step = fr->step;
2008 /* Initiate the energy sums */
2009 edat.steps[nfr] = 1;
2010 edat.points[nfr] = 1;
2011 for (i = 0; i < nset; i++)
2014 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2015 edat.s[i].es[nfr].sum2 = 0;
2022 edat.steps[nfr] = fr->nsteps;
2026 /* mdrun only calculated the energy at energy output
2027 * steps. We don't need to check step intervals.
2029 edat.points[nfr] = 1;
2030 for (i = 0; i < nset; i++)
2033 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2034 edat.s[i].es[nfr].sum2 = 0;
2037 edat.bHaveSums = FALSE;
2039 else if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2041 /* We have statistics to the previous frame */
2042 edat.points[nfr] = fr->nsum;
2043 for (i = 0; i < nset; i++)
2046 edat.s[i].es[nfr].sum = fr->ener[sss].esum;
2047 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2049 edat.npoints += fr->nsum;
2053 /* The interval does not match fr->nsteps:
2054 * can not do exact averages.
2056 edat.bHaveSums = FALSE;
2059 edat.nsteps = fr->step - start_step + 1;
2061 for (i = 0; i < nset; i++)
2063 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2067 * Store energies for analysis afterwards...
2069 if (!bDHDL && (fr->nre > 0))
2071 if (edat.nframes % 1000 == 0)
2073 srenew(time, edat.nframes + 1000);
2075 time[edat.nframes] = fr->t;
2080 do_dhdl(fr, ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2083 /*******************************************
2085 *******************************************/
2092 /* We skip frames with single points (usually only the first frame),
2093 * since they would result in an average plot with outliers.
2097 print_time(out, fr->t);
2098 print1(out, bDp, fr->ener[set[0]].e);
2099 print1(out, bDp, fr->ener[set[0]].esum / fr->nsum);
2100 print1(out, bDp, std::sqrt(fr->ener[set[0]].eav / fr->nsum));
2106 print_time(out, fr->t);
2110 for (i = 0; i < nset; i++)
2112 sum += fr->ener[set[i]].e;
2114 print1(out, bDp, sum / nmol - ezero);
2118 for (i = 0; (i < nset); i++)
2122 print1(out, bDp, (fr->ener[set[i]].e) / nmol - ezero);
2126 print1(out, bDp, fr->ener[set[i]].e);
2135 } while (bCont && (timecheck == 0));
2137 fprintf(stderr, "\n");
2148 gmx_fio_fclose(fp_dhdl);
2149 printf("\n\nWrote %d lambda values with %d samples as ", dh_lambdas, dh_samples);
2152 printf("%d dH histograms ", dh_hists);
2156 printf("%d dH data blocks ", dh_blocks);
2158 printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2162 gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2167 double dt = (frame[cur].t - start_t) / (edat.nframes - 1);
2168 analyse_ener(opt2bSet("-corr", NFILE, fnm),
2169 opt2fn("-corr", NFILE, fnm),
2170 opt2fn("-evisco", NFILE, fnm),
2171 opt2fn("-eviscoi", NFILE, fnm),
2176 opt2fn("-vis", NFILE, fnm),
2196 calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat, nbmin, nbmax);
2199 if (opt2bSet("-f2", NFILE, fnm))
2201 fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm), reftemp, nset, set, leg, &edat, time, oenv);
2204 done_enerdata_t(nset, &edat);
2206 free_enxframe(&frame[0]);
2207 free_enxframe(&frame[1]);
2209 free_enxnms(nre, enm);
2215 const char* nxy = "-nxy";
2217 do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2218 do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2219 do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);
2221 output_env_done(oenv);