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, 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.
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.
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.
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.
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.
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.
44 #include "gromacs/utility/fatalerror.h"
45 #include "gromacs/math/vec.h"
46 #include "gromacs/utility/cstringutil.h"
47 #include "gromacs/utility/smalloc.h"
48 #include "gromacs/fileio/enxio.h"
49 #include "gromacs/commandline/pargs.h"
53 #include "gromacs/fileio/xvgr.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/fileio/tpxio.h"
57 #include "gromacs/fileio/trxio.h"
59 #include "gromacs/topology/mtop_util.h"
63 static real minthird = -1.0/3.0, minsixth = -1.0/6.0;
90 static double mypow(double x, double y)
102 static int *select_it(int nre, char *nm[], int *nset)
107 gmx_bool bVerbose = TRUE;
109 if ((getenv("GMX_ENER_VERBOSE")) != NULL)
114 fprintf(stderr, "Select the terms you want from the following list\n");
115 fprintf(stderr, "End your selection with 0\n");
119 for (k = 0; (k < nre); )
121 for (j = 0; (j < 4) && (k < nre); j++, k++)
123 fprintf(stderr, " %3d=%14s", k+1, nm[k]);
125 fprintf(stderr, "\n");
132 if (1 != scanf("%d", &n))
134 gmx_fatal(FARGS, "Error reading user input");
136 if ((n > 0) && (n <= nre))
144 for (i = (*nset) = 0; (i < nre); i++)
157 static void chomp(char *buf)
159 int len = strlen(buf);
161 while ((len > 0) && (buf[len-1] == '\n'))
168 static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
171 int n, k, kk, j, i, nmatch, nind, nss;
173 gmx_bool bEOF, bVerbose = TRUE, bLong = FALSE;
174 char *ptr, buf[STRLEN];
175 const char *fm4 = "%3d %-14s";
176 const char *fm2 = "%3d %-34s";
179 if ((getenv("GMX_ENER_VERBOSE")) != NULL)
184 fprintf(stderr, "\n");
185 fprintf(stderr, "Select the terms you want from the following list by\n");
186 fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
187 fprintf(stderr, "End your selection with an empty line or a zero.\n");
188 fprintf(stderr, "-------------------------------------------------------------------\n");
192 for (k = 0; k < nre; k++)
194 newnm[k] = strdup(nm[k].name);
195 /* Insert dashes in all the names */
196 while ((ptr = strchr(newnm[k], ' ')) != NULL)
206 fprintf(stderr, "\n");
209 for (kk = k; kk < k+4; kk++)
211 if (kk < nre && strlen(nm[kk].name) > 14)
219 fprintf(stderr, " ");
223 fprintf(stderr, fm4, k+1, newnm[k]);
232 fprintf(stderr, fm2, k+1, newnm[k]);
243 fprintf(stderr, "\n\n");
249 while (!bEOF && (fgets2(buf, STRLEN-1, stdin)))
251 /* Remove newlines */
257 /* Empty line means end of input */
258 bEOF = (strlen(buf) == 0);
266 /* First try to read an integer */
267 nss = sscanf(ptr, "%d", &nind);
270 /* Zero means end of input */
275 else if ((1 <= nind) && (nind <= nre))
281 fprintf(stderr, "number %d is out of range\n", nind);
286 /* Now try to read a string */
289 for (nind = 0; nind < nre; nind++)
291 if (gmx_strcasecmp(newnm[nind], ptr) == 0)
301 for (nind = 0; nind < nre; nind++)
303 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
311 fprintf(stderr, "String '%s' does not match anything\n", ptr);
316 /* Look for the first space, and remove spaces from there */
317 if ((ptr = strchr(ptr, ' ')) != NULL)
322 while (!bEOF && (ptr && (strlen(ptr) > 0)));
327 for (i = (*nset) = 0; (i < nre); i++)
339 gmx_fatal(FARGS, "No energy terms selected");
342 for (i = 0; (i < nre); i++)
351 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
358 /* all we need is the ir to be able to write the label */
359 read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, &mtop);
362 static void get_orires_parms(const char *topnm,
363 int *nor, int *nex, int **label, real **obs)
374 read_tpx(topnm, &ir, box, &natoms, NULL, NULL, NULL, &mtop);
375 top = gmx_mtop_generate_local_top(&mtop, &ir);
377 ip = top->idef.iparams;
378 iatom = top->idef.il[F_ORIRES].iatoms;
380 /* Count how many distance restraint there are... */
381 nb = top->idef.il[F_ORIRES].nr;
384 gmx_fatal(FARGS, "No orientation restraints in topology!\n");
391 for (i = 0; i < nb; i += 3)
393 (*label)[i/3] = ip[iatom[i]].orires.label;
394 (*obs)[i/3] = ip[iatom[i]].orires.obs;
395 if (ip[iatom[i]].orires.ex >= *nex)
397 *nex = ip[iatom[i]].orires.ex+1;
400 fprintf(stderr, "Found %d orientation restraints with %d experiments",
404 static int get_bounds(const char *topnm,
405 real **bounds, int **index, int **dr_pair, int *npairs,
406 gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir)
409 t_functype *functype;
411 int natoms, i, j, k, type, ftype, natom;
419 read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, mtop);
421 top = gmx_mtop_generate_local_top(mtop, ir);
424 functype = top->idef.functype;
425 ip = top->idef.iparams;
427 /* Count how many distance restraint there are... */
428 nb = top->idef.il[F_DISRES].nr;
431 gmx_fatal(FARGS, "No distance restraints in topology!\n");
434 /* Allocate memory */
439 /* Fill the bound array */
441 for (i = 0; (i < top->idef.ntypes); i++)
444 if (ftype == F_DISRES)
447 label1 = ip[i].disres.label;
448 b[nb] = ip[i].disres.up1;
455 /* Fill the index array */
457 disres = &(top->idef.il[F_DISRES]);
458 iatom = disres->iatoms;
459 for (i = j = k = 0; (i < disres->nr); )
462 ftype = top->idef.functype[type];
463 natom = interaction_function[ftype].nratoms+1;
464 if (label1 != top->idef.iparams[type].disres.label)
467 label1 = top->idef.iparams[type].disres.label;
477 gmx_incons("get_bounds for distance restraints");
486 static void calc_violations(real rt[], real rav3[], int nb, int index[],
487 real bounds[], real *viol, double *st, double *sa)
489 const real sixth = 1.0/6.0;
491 double rsum, rav, sumaver, sumt;
495 for (i = 0; (i < nb); i++)
499 for (j = index[i]; (j < index[i+1]); j++)
503 viol[j] += mypow(rt[j], -3.0);
506 rsum += mypow(rt[j], -6);
508 rsum = max(0.0, mypow(rsum, -sixth)-bounds[i]);
509 rav = max(0.0, mypow(rav, -sixth)-bounds[i]);
518 static void analyse_disre(const char *voutfn, int nframes,
519 real violaver[], real bounds[], int index[],
520 int pair[], int nbounds,
521 const output_env_t oenv)
524 double sum, sumt, sumaver;
527 /* Subtract bounds from distances, to calculate violations */
528 calc_violations(violaver, violaver,
529 nbounds, pair, bounds, NULL, &sumt, &sumaver);
532 fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
534 fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n",
537 vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm",
541 for (i = 0; (i < nbounds); i++)
543 /* Do ensemble averaging */
545 for (j = pair[i]; (j < pair[i+1]); j++)
547 sumaver += sqr(violaver[j]/nframes);
549 sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]);
552 sum = max(sum, sumaver);
553 fprintf(vout, "%10d %10.5e\n", index[i], sumaver);
556 for (j = 0; (j < dr.ndr); j++)
558 fprintf(vout, "%10d %10.5e\n", j, mypow(violaver[j]/nframes, minthird));
563 fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
565 fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
567 do_view(oenv, voutfn, "-graphtype bar");
570 static void einstein_visco(const char *fn, const char *fni, int nsets,
571 int nint, real **eneint,
572 real V, real T, double dt,
573 const output_env_t oenv)
576 double av[4], avold[4];
582 for (i = 0; i <= nsets; i++)
586 fp0 = xvgropen(fni, "Shear viscosity integral",
587 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
588 fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
589 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
590 for (i = 0; i < nf4; i++)
592 for (m = 0; m <= nsets; m++)
596 for (j = 0; j < nint - i; j++)
598 for (m = 0; m < nsets; m++)
600 di = sqr(eneint[m][j+i] - eneint[m][j]);
603 av[nsets] += di/nsets;
606 /* Convert to SI for the viscosity */
607 fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nint - i);
608 fprintf(fp0, "%10g", i*dt);
609 for (m = 0; (m <= nsets); m++)
612 fprintf(fp0, " %10g", av[m]);
615 fprintf(fp1, "%10g", (i + 0.5)*dt);
616 for (m = 0; (m <= nsets); m++)
618 fprintf(fp1, " %10g", (av[m]-avold[m])/dt);
641 static void clear_ee_sum(ee_sum_t *ees)
649 static void add_ee_sum(ee_sum_t *ees, double sum, int np)
655 static void add_ee_av(ee_sum_t *ees)
659 av = ees->sum/ees->np;
666 static double calc_ee2(int nb, ee_sum_t *ees)
668 return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
671 static void set_ee_av(ener_ee_t *eee)
675 char buf[STEPSTRSIZE];
676 fprintf(debug, "Storing average for err.est.: %s steps\n",
677 gmx_step_str(eee->nst, buf));
679 add_ee_av(&eee->sum);
681 if (eee->b == 1 || eee->nst < eee->nst_min)
683 eee->nst_min = eee->nst;
688 static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
691 double sum, sum2, sump, see2;
692 gmx_int64_t steps, np, p, bound_nb;
696 double x, sx, sy, sxx, sxy;
699 /* Check if we have exact statistics over all points */
700 for (i = 0; i < nset; i++)
703 ed->bExactStat = FALSE;
704 if (edat->npoints > 0)
706 /* All energy file sum entries 0 signals no exact sums.
707 * But if all energy values are 0, we still have exact sums.
710 for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
712 if (ed->ener[i] != 0)
716 ed->bExactStat = (ed->es[f].sum != 0);
720 ed->bExactStat = TRUE;
726 for (i = 0; i < nset; i++)
737 for (nb = nbmin; nb <= nbmax; nb++)
740 clear_ee_sum(&eee[nb].sum);
744 for (f = 0; f < edat->nframes; f++)
750 /* Add the sum and the sum of variances to the totals. */
756 sum2 += dsqr(sum/np - (sum + es->sum)/(np + p))
762 /* Add a single value to the sum and sum of squares. */
768 /* sum has to be increased after sum2 */
772 /* For the linear regression use variance 1/p.
773 * Note that sump is the sum, not the average, so we don't need p*.
775 x = edat->step[f] - 0.5*(edat->steps[f] - 1);
781 for (nb = nbmin; nb <= nbmax; nb++)
783 /* Check if the current end step is closer to the desired
784 * block boundary than the next end step.
786 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
787 if (eee[nb].nst > 0 &&
788 bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
798 eee[nb].nst += edat->step[f] - edat->step[f-1];
802 add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
806 add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
808 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
809 if (edat->step[f]*nb >= bound_nb)
816 edat->s[i].av = sum/np;
819 edat->s[i].rmsd = sqrt(sum2/np);
823 edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av));
826 if (edat->nframes > 1)
828 edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
832 edat->s[i].slope = 0;
837 for (nb = nbmin; nb <= nbmax; nb++)
839 /* Check if we actually got nb blocks and if the smallest
840 * block is not shorter than 80% of the average.
844 char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
845 fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
847 gmx_step_str(eee[nb].nst_min, buf1),
848 gmx_step_str(edat->nsteps, buf2));
850 if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
852 see2 += calc_ee2(nb, &eee[nb].sum);
858 edat->s[i].ee = sqrt(see2/nee);
868 static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax)
879 snew(s->ener, esum->nframes);
880 snew(s->es, esum->nframes);
882 s->bExactStat = TRUE;
884 for (i = 0; i < nset; i++)
886 if (!edat->s[i].bExactStat)
888 s->bExactStat = FALSE;
890 s->slope += edat->s[i].slope;
893 for (f = 0; f < edat->nframes; f++)
896 for (i = 0; i < nset; i++)
898 sum += edat->s[i].ener[f];
902 for (i = 0; i < nset; i++)
904 sum += edat->s[i].es[f].sum;
910 calc_averages(1, esum, nbmin, nbmax);
915 static char *ee_pr(double ee, char *buf)
922 sprintf(buf, "%s", "--");
926 /* Round to two decimals by printing. */
927 sprintf(tmp, "%.1e", ee);
928 sscanf(tmp, "%lf", &rnd);
929 sprintf(buf, "%g", rnd);
935 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat)
937 /* Remove the drift by performing a fit to y = ax+b.
938 Uses 5 iterations. */
940 double delta, d, sd, sd2;
942 edat->npoints = edat->nframes;
943 edat->nsteps = edat->nframes;
945 for (k = 0; (k < 5); k++)
947 for (i = 0; (i < nset); i++)
949 delta = edat->s[i].slope*dt;
953 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
956 for (j = 0; (j < edat->nframes); j++)
958 edat->s[i].ener[j] -= j*delta;
959 edat->s[i].es[j].sum = 0;
960 edat->s[i].es[j].sum2 = 0;
963 calc_averages(nset, edat, nbmin, nbmax);
967 static void calc_fluctuation_props(FILE *fp,
968 gmx_bool bDriftCorr, real dt,
970 char **leg, enerdata_t *edat,
971 int nbmin, int nbmax)
974 double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
977 eVol, eEnth, eTemp, eEtot, eNR
979 const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
982 NANO3 = NANO*NANO*NANO;
985 fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n"
986 "for spurious drift in the graphs. Note that this is not\n"
987 "a substitute for proper equilibration and sampling!\n");
991 remove_drift(nset, nbmin, nbmax, dt, edat);
993 for (i = 0; (i < eNR); i++)
995 for (ii[i] = 0; (ii[i] < nset &&
996 (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++)
1000 /* if (ii[i] < nset)
1001 fprintf(fp,"Found %s data.\n",my_ener[i]);
1003 /* Compute it all! */
1004 alpha = kappa = cp = dcp = cv = NOTSET;
1008 if (ii[eTemp] < nset)
1010 tt = edat->s[ii[eTemp]].av;
1014 if ((ii[eVol] < nset) && (ii[eTemp] < nset))
1016 vv = edat->s[ii[eVol]].av*NANO3;
1017 varv = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
1018 kappa = (varv/vv)/(BOLTZMANN*tt);
1022 if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
1024 hh = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
1025 varh = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
1026 cp = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
1029 et = varet = NOTSET;
1030 if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
1032 /* Only compute cv in constant volume runs, which we can test
1033 by checking whether the enthalpy was computed.
1035 et = edat->s[ii[eEtot]].av;
1036 varet = sqr(edat->s[ii[eEtot]].rmsd);
1037 cv = KILO*((varet/nmol)/(BOLTZ*tt*tt));
1040 if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
1042 double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
1043 vh_sum = v_sum = h_sum = 0;
1044 for (j = 0; (j < edat->nframes); j++)
1046 v = edat->s[ii[eVol]].ener[j]*NANO3;
1047 h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
1052 vh_aver = vh_sum / edat->nframes;
1053 v_aver = v_sum / edat->nframes;
1054 h_aver = h_sum / edat->nframes;
1055 alpha = (vh_aver-v_aver*h_aver)/(v_aver*BOLTZMANN*tt*tt);
1056 dcp = (v_aver*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
1063 fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n",
1066 fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
1067 fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
1068 fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
1069 fprintf(fp, "please use the g_dos program.\n\n");
1070 fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n"
1071 "a block-averaging error analysis (not implemented in g_energy yet)\n");
1077 fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol);
1082 fprintf(fp, "Volume = %10g m^3/mol\n",
1087 fprintf(fp, "Enthalpy = %10g kJ/mol\n",
1088 hh*AVOGADRO/(KILO*nmol));
1090 if (alpha != NOTSET)
1092 fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
1095 if (kappa != NOTSET)
1097 fprintf(fp, "Isothermal Compressibility Kappa = %10g (J/m^3)\n",
1099 fprintf(fp, "Adiabatic bulk modulus = %10g (m^3/J)\n",
1104 fprintf(fp, "Heat capacity at constant pressure Cp = %10g J/mol K\n",
1109 fprintf(fp, "Heat capacity at constant volume Cv = %10g J/mol K\n",
1114 fprintf(fp, "Cp-Cv = %10g J/mol K\n",
1117 please_cite(fp, "Allen1987a");
1121 fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
1125 static void analyse_ener(gmx_bool bCorr, const char *corrfn,
1126 gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct,
1127 gmx_bool bVisco, const char *visfn, int nmol,
1128 gmx_int64_t start_step, double start_t,
1129 gmx_int64_t step, double t,
1132 int nset, int set[], gmx_bool *bIsEner,
1133 char **leg, gmx_enxnm_t *enm,
1134 real Vaver, real ezero,
1135 int nbmin, int nbmax,
1136 const output_env_t oenv)
1139 /* Check out the printed manual for equations! */
1140 double Dt, aver, stddev, errest, delta_t, totaldrift;
1141 enerdata_t *esum = NULL;
1142 real xxx, integral, intBulk, Temp = 0, Pres = 0;
1143 real sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest;
1144 double beta = 0, expE, expEtot, *fee = NULL;
1146 int nexact, nnotexact;
1150 char buf[256], eebuf[100];
1152 nsteps = step - start_step + 1;
1155 fprintf(stdout, "Not enough steps (%s) for statistics\n",
1156 gmx_step_str(nsteps, buf));
1160 /* Calculate the time difference */
1161 delta_t = t - start_t;
1163 fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
1164 gmx_step_str(nsteps, buf), start_t, t, nset);
1166 calc_averages(nset, edat, nbmin, nbmax);
1170 esum = calc_sum(nset, edat, nbmin, nbmax);
1173 if (edat->npoints == 0)
1182 for (i = 0; (i < nset); i++)
1184 if (edat->s[i].bExactStat)
1197 fprintf(stdout, "All statistics are over %s points\n",
1198 gmx_step_str(edat->npoints, buf));
1200 else if (nexact == 0 || edat->npoints == edat->nframes)
1202 fprintf(stdout, "All statistics are over %d points (frames)\n",
1207 fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
1208 for (i = 0; (i < nset); i++)
1210 if (!edat->s[i].bExactStat)
1212 fprintf(stdout, " '%s'", leg[i]);
1215 fprintf(stdout, " %s has statistics over %d points (frames)\n",
1216 nnotexact == 1 ? "is" : "are", edat->nframes);
1217 fprintf(stdout, "All other statistics are over %s points\n",
1218 gmx_step_str(edat->npoints, buf));
1220 fprintf(stdout, "\n");
1222 fprintf(stdout, "%-24s %10s %10s %10s %10s",
1223 "Energy", "Average", "Err.Est.", "RMSD", "Tot-Drift");
1226 fprintf(stdout, " %10s\n", "-kT ln<e^(E/kT)>");
1230 fprintf(stdout, "\n");
1232 fprintf(stdout, "-------------------------------------------------------------------------------\n");
1234 /* Initiate locals, only used with -sum */
1238 beta = 1.0/(BOLTZ*reftemp);
1241 for (i = 0; (i < nset); i++)
1243 aver = edat->s[i].av;
1244 stddev = edat->s[i].rmsd;
1245 errest = edat->s[i].ee;
1250 for (j = 0; (j < edat->nframes); j++)
1252 expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
1256 expEtot += expE/edat->nframes;
1259 fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
1261 if (strstr(leg[i], "empera") != NULL)
1265 else if (strstr(leg[i], "olum") != NULL)
1269 else if (strstr(leg[i], "essure") != NULL)
1275 pr_aver = aver/nmol-ezero;
1276 pr_stddev = stddev/nmol;
1277 pr_errest = errest/nmol;
1286 /* Multiply the slope in steps with the number of steps taken */
1287 totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
1293 fprintf(stdout, "%-24s %10g %10s %10g %10g",
1294 leg[i], pr_aver, ee_pr(pr_errest, eebuf), pr_stddev, totaldrift);
1297 fprintf(stdout, " %10g", fee[i]);
1300 fprintf(stdout, " (%s)\n", enm[set[i]].unit);
1304 for (j = 0; (j < edat->nframes); j++)
1306 edat->s[i].ener[j] -= aver;
1312 totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1313 fprintf(stdout, "%-24s %10g %10s %10s %10g (%s)",
1314 "Total", esum->s[0].av/nmol, ee_pr(esum->s[0].ee/nmol, eebuf),
1315 "--", totaldrift/nmol, enm[set[0]].unit);
1316 /* pr_aver,pr_stddev,a,totaldrift */
1319 fprintf(stdout, " %10g %10g\n",
1320 log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta);
1324 fprintf(stdout, "\n");
1328 /* Do correlation function */
1329 if (edat->nframes > 1)
1331 Dt = delta_t/(edat->nframes - 1);
1339 const char* leg[] = { "Shear", "Bulk" };
1344 /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1346 /* Symmetrise tensor! (and store in first three elements)
1347 * And subtract average pressure!
1350 for (i = 0; i < 12; i++)
1352 snew(eneset[i], edat->nframes);
1354 for (i = 0; (i < edat->nframes); i++)
1356 eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1357 eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1358 eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1359 for (j = 3; j <= 11; j++)
1361 eneset[j][i] = edat->s[j].ener[i];
1363 eneset[11][i] -= Pres;
1366 /* Determine integrals of the off-diagonal pressure elements */
1368 for (i = 0; i < 3; i++)
1370 snew(eneint[i], edat->nframes + 1);
1375 for (i = 0; i < edat->nframes; i++)
1377 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];
1378 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];
1379 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];
1382 einstein_visco("evisco.xvg", "eviscoi.xvg",
1383 3, edat->nframes+1, eneint, Vaver, Temp, Dt, oenv);
1385 for (i = 0; i < 3; i++)
1391 /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1392 /* Do it for shear viscosity */
1393 strcpy(buf, "Shear Viscosity");
1394 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3,
1395 (edat->nframes+1)/2, eneset, Dt,
1396 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1398 /* Now for bulk viscosity */
1399 strcpy(buf, "Bulk Viscosity");
1400 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1,
1401 (edat->nframes+1)/2, &(eneset[11]), Dt,
1402 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1404 factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1405 fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1406 xvgr_legend(fp, asize(leg), leg, oenv);
1408 /* Use trapezium rule for integration */
1411 nout = get_acfnout();
1412 if ((nout < 2) || (nout >= edat->nframes/2))
1414 nout = edat->nframes/2;
1416 for (i = 1; (i < nout); i++)
1418 integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor;
1419 intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1420 fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk);
1424 for (i = 0; i < 12; i++)
1434 strcpy(buf, "Autocorrelation of Energy Fluctuations");
1438 strcpy(buf, "Energy Autocorrelation");
1441 do_autocorr(corrfn, oenv, buf, edat->nframes,
1443 bSum ? &edat->s[nset-1].ener : eneset,
1444 (delta_t/edat->nframes), eacNormal, FALSE);
1450 static void print_time(FILE *fp, double t)
1452 fprintf(fp, "%12.6f", t);
1455 static void print1(FILE *fp, gmx_bool bDp, real e)
1459 fprintf(fp, " %16.12f", e);
1463 fprintf(fp, " %10.6f", e);
1467 static void fec(const char *ene2fn, const char *runavgfn,
1468 real reftemp, int nset, int set[], char *leg[],
1469 enerdata_t *edat, double time[],
1470 const output_env_t oenv)
1472 const char * ravgleg[] = {
1473 "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1474 "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N"
1478 int nre, timecheck, step, nenergy, nenergy2, maxenergy;
1484 gmx_enxnm_t *enm = NULL;
1488 /* read second energy file */
1491 enx = open_enx(ene2fn, "r");
1492 do_enxnms(enx, &(fr->nre), &enm);
1494 snew(eneset2, nset+1);
1500 /* This loop searches for the first frame (when -b option is given),
1501 * or when this has been found it reads just one energy frame
1505 bCont = do_enx(enx, fr);
1509 timecheck = check_times(fr->t);
1513 while (bCont && (timecheck < 0));
1515 /* Store energies for analysis afterwards... */
1516 if ((timecheck == 0) && bCont)
1520 if (nenergy2 >= maxenergy)
1523 for (i = 0; i <= nset; i++)
1525 srenew(eneset2[i], maxenergy);
1528 if (fr->t != time[nenergy2])
1530 fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n",
1531 fr->t, time[nenergy2], gmx_step_str(fr->step, buf));
1533 for (i = 0; i < nset; i++)
1535 eneset2[i][nenergy2] = fr->ener[set[i]].e;
1541 while (bCont && (timecheck == 0));
1544 if (edat->nframes != nenergy2)
1546 fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n",
1547 edat->nframes, nenergy2);
1549 nenergy = min(edat->nframes, nenergy2);
1551 /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1555 fp = xvgropen(runavgfn, "Running average free energy difference",
1556 "Time (" unit_time ")", "\\8D\\4E (" unit_energy ")", oenv);
1557 xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
1559 fprintf(stdout, "\n%-24s %10s\n",
1560 "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
1562 beta = 1.0/(BOLTZ*reftemp);
1563 for (i = 0; i < nset; i++)
1565 if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
1567 fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n",
1568 leg[i], enm[set[i]].name);
1570 for (j = 0; j < nenergy; j++)
1572 dE = eneset2[i][j] - edat->s[i].ener[j];
1573 sum += exp(-dE*beta);
1576 fprintf(fp, "%10g %10g %10g\n",
1577 time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
1580 aver = -BOLTZ*reftemp*log(sum/nenergy);
1581 fprintf(stdout, "%-24s %10g\n", leg[i], aver);
1591 static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl,
1592 const char *filename, gmx_bool bDp,
1593 int *blocks, int *hists, int *samples, int *nlambdas,
1594 const output_env_t oenv)
1596 const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1597 char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
1599 gmx_bool first = FALSE;
1600 int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
1603 double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0;
1604 static int setnr = 0;
1605 double *native_lambda_vec = NULL;
1606 const char **lambda_components = NULL;
1607 int n_lambda_vec = 0;
1608 gmx_bool changing_lambda = FALSE;
1609 int lambda_fep_state;
1611 /* now count the blocks & handle the global dh data */
1612 for (i = 0; i < fr->nblock; i++)
1614 if (fr->block[i].id == enxDHHIST)
1618 else if (fr->block[i].id == enxDH)
1622 else if (fr->block[i].id == enxDHCOLL)
1625 if ( (fr->block[i].nsub < 1) ||
1626 (fr->block[i].sub[0].type != xdr_datatype_double) ||
1627 (fr->block[i].sub[0].nr < 5))
1629 gmx_fatal(FARGS, "Unexpected block data");
1632 /* read the data from the DHCOLL block */
1633 temp = fr->block[i].sub[0].dval[0];
1634 start_time = fr->block[i].sub[0].dval[1];
1635 delta_time = fr->block[i].sub[0].dval[2];
1636 start_lambda = fr->block[i].sub[0].dval[3];
1637 delta_lambda = fr->block[i].sub[0].dval[4];
1638 changing_lambda = (delta_lambda != 0);
1639 if (fr->block[i].nsub > 1)
1641 lambda_fep_state = fr->block[i].sub[1].ival[0];
1642 if (n_lambda_vec == 0)
1644 n_lambda_vec = fr->block[i].sub[1].ival[1];
1648 if (n_lambda_vec != fr->block[i].sub[1].ival[1])
1651 "Unexpected change of basis set in lambda");
1654 if (lambda_components == NULL)
1656 snew(lambda_components, n_lambda_vec);
1658 if (native_lambda_vec == NULL)
1660 snew(native_lambda_vec, n_lambda_vec);
1662 for (j = 0; j < n_lambda_vec; j++)
1664 native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j];
1665 lambda_components[j] =
1666 efpt_singular_names[fr->block[i].sub[1].ival[2+j]];
1672 if (nblock_hist == 0 && nblock_dh == 0)
1674 /* don't do anything */
1677 if (nblock_hist > 0 && nblock_dh > 0)
1679 gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1685 /* we have standard, non-histogram data --
1686 call open_dhdl to open the file */
1687 /* TODO this is an ugly hack that needs to be fixed: this will only
1688 work if the order of data is always the same and if we're
1689 only using the g_energy compiled with the mdrun that produced
1691 *fp_dhdl = open_dhdl(filename, ir, oenv);
1695 sprintf(title, "N(%s)", deltag);
1696 sprintf(label_x, "%s (%s)", deltag, unit_energy);
1697 sprintf(label_y, "Samples");
1698 *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1699 sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
1700 xvgr_subtitle(*fp_dhdl, buf, oenv);
1704 (*hists) += nblock_hist;
1705 (*blocks) += nblock_dh;
1706 (*nlambdas) = nblock_hist+nblock_dh;
1708 /* write the data */
1709 if (nblock_hist > 0)
1711 gmx_int64_t sum = 0;
1713 for (i = 0; i < fr->nblock; i++)
1715 t_enxblock *blk = &(fr->block[i]);
1716 if (blk->id == enxDHHIST)
1718 double foreign_lambda, dx;
1720 int nhist, derivative;
1722 /* check the block types etc. */
1723 if ( (blk->nsub < 2) ||
1724 (blk->sub[0].type != xdr_datatype_double) ||
1725 (blk->sub[1].type != xdr_datatype_int64) ||
1726 (blk->sub[0].nr < 2) ||
1727 (blk->sub[1].nr < 2) )
1729 gmx_fatal(FARGS, "Unexpected block data in file");
1731 foreign_lambda = blk->sub[0].dval[0];
1732 dx = blk->sub[0].dval[1];
1733 nhist = blk->sub[1].lval[0];
1734 derivative = blk->sub[1].lval[1];
1735 for (j = 0; j < nhist; j++)
1738 x0 = blk->sub[1].lval[2+j];
1742 sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1743 deltag, lambda, foreign_lambda,
1744 lambda, start_lambda);
1748 sprintf(legend, "N(%s | %s=%g)",
1749 dhdl, lambda, start_lambda);
1753 xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1755 for (k = 0; k < blk->sub[j+2].nr; k++)
1760 hist = blk->sub[j+2].ival[k];
1763 fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist,
1767 /* multiple histogram data blocks in one histogram
1768 mean that the second one is the reverse of the first one:
1769 for dhdl derivatives, it's important to know both the
1770 maximum and minimum values */
1775 (*samples) += (int)(sum/nblock_hist);
1781 char **setnames = NULL;
1782 int nnames = nblock_dh;
1784 for (i = 0; i < fr->nblock; i++)
1786 t_enxblock *blk = &(fr->block[i]);
1787 if (blk->id == enxDH)
1791 len = blk->sub[2].nr;
1795 if (len != blk->sub[2].nr)
1797 gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1804 for (i = 0; i < len; i++)
1806 double time = start_time + delta_time*i;
1808 fprintf(*fp_dhdl, "%.4f ", time);
1810 for (j = 0; j < fr->nblock; j++)
1812 t_enxblock *blk = &(fr->block[j]);
1813 if (blk->id == enxDH)
1816 if (blk->sub[2].type == xdr_datatype_float)
1818 value = blk->sub[2].fval[i];
1822 value = blk->sub[2].dval[i];
1824 /* we need to decide which data type it is based on the count*/
1826 if (j == 1 && ir->bExpanded)
1828 fprintf(*fp_dhdl, "%4d", (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! */
1834 fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */
1838 fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */
1843 fprintf(*fp_dhdl, "\n");
1849 int gmx_energy(int argc, char *argv[])
1851 const char *desc[] = {
1852 "[THISMODULE] extracts energy components or distance restraint",
1853 "data from an energy file. The user is prompted to interactively",
1854 "select the desired energy terms.[PAR]",
1856 "Average, RMSD, and drift are calculated with full precision from the",
1857 "simulation (see printed manual). Drift is calculated by performing",
1858 "a least-squares fit of the data to a straight line. The reported total drift",
1859 "is the difference of the fit at the first and last point.",
1860 "An error estimate of the average is given based on a block averages",
1861 "over 5 blocks using the full-precision averages. The error estimate",
1862 "can be performed over multiple block lengths with the options",
1863 "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1864 "[BB]Note[bb] that in most cases the energy files contains averages over all",
1865 "MD steps, or over many more points than the number of frames in",
1866 "energy file. This makes the [THISMODULE] statistics output more accurate",
1867 "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1868 "file, the statistics mentioned above are simply over the single, per-frame",
1869 "energy values.[PAR]",
1871 "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1873 "Some fluctuation-dependent properties can be calculated provided",
1874 "the correct energy terms are selected, and that the command line option",
1875 "[TT]-fluct_props[tt] is given. The following properties",
1876 "will be computed:[BR]",
1877 "Property Energy terms needed[BR]",
1878 "---------------------------------------------------[BR]",
1879 "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp [BR]",
1880 "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp [BR]",
1881 "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1882 "Isothermal compressibility: Vol, Temp [BR]",
1883 "Adiabatic bulk modulus: Vol, Temp [BR]",
1884 "---------------------------------------------------[BR]",
1885 "You always need to set the number of molecules [TT]-nmol[tt].",
1886 "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1887 "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]"
1888 "When the [TT]-viol[tt] option is set, the time averaged",
1889 "violations are plotted and the running time-averaged and",
1890 "instantaneous sum of violations are recalculated. Additionally",
1891 "running time-averaged and instantaneous distances between",
1892 "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1894 "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1895 "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1896 "The first two options plot the orientation, the last three the",
1897 "deviations of the orientations from the experimental values.",
1898 "The options that end on an 'a' plot the average over time",
1899 "as a function of restraint. The options that end on a 't'",
1900 "prompt the user for restraint label numbers and plot the data",
1901 "as a function of time. Option [TT]-odr[tt] plots the RMS",
1902 "deviation as a function of restraint.",
1903 "When the run used time or ensemble averaged orientation restraints,",
1904 "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1905 "not ensemble-averaged orientations and deviations instead of",
1906 "the time and ensemble averages.[PAR]",
1908 "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1909 "tensor for each orientation restraint experiment. With option",
1910 "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1912 "Option [TT]-odh[tt] extracts and plots the free energy data",
1913 "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1914 "from the [TT]ener.edr[tt] file.[PAR]",
1916 "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1917 "difference with an ideal gas state: [BR]",
1918 " [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][BR]",
1919 " [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][BR]",
1920 "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1921 "the average is over the ensemble (or time in a trajectory).",
1922 "Note that this is in principle",
1923 "only correct when averaging over the whole (Boltzmann) ensemble",
1924 "and using the potential energy. This also allows for an entropy",
1925 "estimate using:[BR]",
1926 " [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[BR]",
1927 " [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",
1930 "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1931 "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1932 "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1933 "files, and the average is over the ensemble A. The running average",
1934 "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1935 "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1938 static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1939 static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1940 static int skip = 0, nmol = 1, nbmin = 5, nbmax = 5;
1941 static real reftemp = 300.0, ezero = 0;
1943 { "-fee", FALSE, etBOOL, {&bFee},
1944 "Do a free energy estimate" },
1945 { "-fetemp", FALSE, etREAL, {&reftemp},
1946 "Reference temperature for free energy calculation" },
1947 { "-zero", FALSE, etREAL, {&ezero},
1948 "Subtract a zero-point energy" },
1949 { "-sum", FALSE, etBOOL, {&bSum},
1950 "Sum the energy terms selected rather than display them all" },
1951 { "-dp", FALSE, etBOOL, {&bDp},
1952 "Print energies in high precision" },
1953 { "-nbmin", FALSE, etINT, {&nbmin},
1954 "Minimum number of blocks for error estimate" },
1955 { "-nbmax", FALSE, etINT, {&nbmax},
1956 "Maximum number of blocks for error estimate" },
1957 { "-mutot", FALSE, etBOOL, {&bMutot},
1958 "Compute the total dipole moment from the components" },
1959 { "-skip", FALSE, etINT, {&skip},
1960 "Skip number of frames between data points" },
1961 { "-aver", FALSE, etBOOL, {&bPrAll},
1962 "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1963 { "-nmol", FALSE, etINT, {&nmol},
1964 "Number of molecules in your sample: the energies are divided by this number" },
1965 { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1966 "Compute properties based on energy fluctuations, like heat capacity" },
1967 { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1968 "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1969 { "-fluc", FALSE, etBOOL, {&bFluct},
1970 "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1971 { "-orinst", FALSE, etBOOL, {&bOrinst},
1972 "Analyse instantaneous orientation data" },
1973 { "-ovec", FALSE, etBOOL, {&bOvec},
1974 "Also plot the eigenvectors with [TT]-oten[tt]" }
1976 const char * drleg[] = {
1980 static const char *setnm[] = {
1981 "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1982 "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1983 "Volume", "Pressure"
1986 FILE *out = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL;
1987 FILE *fp_dhdl = NULL;
1992 gmx_localtop_t *top = NULL;
1996 gmx_enxnm_t *enm = NULL;
1997 t_enxframe *frame, *fr = NULL;
1999 #define NEXT (1-cur)
2000 int nre, teller, teller_disre, nfr;
2001 gmx_int64_t start_step;
2002 int nor = 0, nex = 0, norfr = 0, enx_i = 0;
2004 real *bounds = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL;
2005 int *index = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL;
2006 int nbounds = 0, npairs;
2007 gmx_bool bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL;
2008 gmx_bool bFoundStart, bCont, bEDR, bVisco;
2009 double sum, sumaver, sumt, ener, dbl;
2010 double *time = NULL;
2012 int *set = NULL, i, j, k, nset, sss;
2013 gmx_bool *bIsEner = NULL;
2014 char **pairleg, **odtleg, **otenleg;
2017 char *anm_j, *anm_k, *resnm_j, *resnm_k;
2018 int resnr_j, resnr_k;
2019 const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
2022 t_enxblock *blk = NULL;
2023 t_enxblock *blk_disre = NULL;
2025 int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
2028 { efEDR, "-f", NULL, ffREAD },
2029 { efEDR, "-f2", NULL, ffOPTRD },
2030 { efTPX, "-s", NULL, ffOPTRD },
2031 { efXVG, "-o", "energy", ffWRITE },
2032 { efXVG, "-viol", "violaver", ffOPTWR },
2033 { efXVG, "-pairs", "pairs", ffOPTWR },
2034 { efXVG, "-ora", "orienta", ffOPTWR },
2035 { efXVG, "-ort", "orientt", ffOPTWR },
2036 { efXVG, "-oda", "orideva", ffOPTWR },
2037 { efXVG, "-odr", "oridevr", ffOPTWR },
2038 { efXVG, "-odt", "oridevt", ffOPTWR },
2039 { efXVG, "-oten", "oriten", ffOPTWR },
2040 { efXVG, "-corr", "enecorr", ffOPTWR },
2041 { efXVG, "-vis", "visco", ffOPTWR },
2042 { efXVG, "-ravg", "runavgdf", ffOPTWR },
2043 { efXVG, "-odh", "dhdl", ffOPTWR }
2045 #define NFILE asize(fnm)
2050 ppa = add_acf_pargs(&npargs, pa);
2051 if (!parse_common_args(&argc, argv,
2052 PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
2053 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
2058 bDRAll = opt2bSet("-pairs", NFILE, fnm);
2059 bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
2060 bORA = opt2bSet("-ora", NFILE, fnm);
2061 bORT = opt2bSet("-ort", NFILE, fnm);
2062 bODA = opt2bSet("-oda", NFILE, fnm);
2063 bODR = opt2bSet("-odr", NFILE, fnm);
2064 bODT = opt2bSet("-odt", NFILE, fnm);
2065 bORIRE = bORA || bORT || bODA || bODR || bODT;
2066 bOTEN = opt2bSet("-oten", NFILE, fnm);
2067 bDHDL = opt2bSet("-odh", NFILE, fnm);
2072 fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
2073 do_enxnms(fp, &nre, &enm);
2077 bVisco = opt2bSet("-vis", NFILE, fnm);
2079 if ((!bDisRe) && (!bDHDL))
2083 nset = asize(setnm);
2085 /* This is nasty code... To extract Pres tensor, Volume and Temperature */
2086 for (j = 0; j < nset; j++)
2088 for (i = 0; i < nre; i++)
2090 if (strstr(enm[i].name, setnm[j]))
2098 if (gmx_strcasecmp(setnm[j], "Volume") == 0)
2100 printf("Enter the box volume (" unit_volume "): ");
2101 if (1 != scanf("%lf", &dbl))
2103 gmx_fatal(FARGS, "Error reading user input");
2109 gmx_fatal(FARGS, "Could not find term %s for viscosity calculation",
2117 set = select_by_name(nre, enm, &nset);
2119 /* Print all the different units once */
2120 sprintf(buf, "(%s)", enm[set[0]].unit);
2121 for (i = 1; i < nset; i++)
2123 for (j = 0; j < i; j++)
2125 if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
2133 strcat(buf, enm[set[i]].unit);
2137 out = xvgropen(opt2fn("-o", NFILE, fnm), "Gromacs Energies", "Time (ps)", buf,
2141 for (i = 0; (i < nset); i++)
2143 leg[i] = enm[set[i]].name;
2147 leg[nset] = strdup("Sum");
2148 xvgr_legend(out, nset+1, (const char**)leg, oenv);
2152 xvgr_legend(out, nset, (const char**)leg, oenv);
2155 snew(bIsEner, nset);
2156 for (i = 0; (i < nset); i++)
2159 for (j = 0; (j <= F_ETOT); j++)
2161 bIsEner[i] = bIsEner[i] ||
2162 (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
2166 if (bPrAll && nset > 1)
2168 gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
2173 if (bORIRE || bOTEN)
2175 get_orires_parms(ftp2fn(efTPX, NFILE, fnm), &nor, &nex, &or_label, &oobs);
2199 fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n");
2200 fprintf(stderr, "End your selection with 0\n");
2207 if (1 != scanf("%d", &(orsel[j])))
2209 gmx_fatal(FARGS, "Error reading user input");
2212 while (orsel[j] > 0);
2215 fprintf(stderr, "Selecting all %d orientation restraints\n", nor);
2218 for (i = 0; i < nor; i++)
2225 /* Build the selection */
2227 for (i = 0; i < j; i++)
2229 for (k = 0; k < nor; k++)
2231 if (or_label[k] == orsel[i])
2240 fprintf(stderr, "Orientation restraint label %d not found\n",
2245 snew(odtleg, norsel);
2246 for (i = 0; i < norsel; i++)
2248 snew(odtleg[i], 256);
2249 sprintf(odtleg[i], "%d", or_label[orsel[i]]);
2253 fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
2254 "Time (ps)", "", oenv);
2255 if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2257 fprintf(fort, "%s", orinst_sub);
2259 xvgr_legend(fort, norsel, (const char**)odtleg, oenv);
2263 fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
2264 "Orientation restraint deviation",
2265 "Time (ps)", "", oenv);
2266 if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2268 fprintf(fodt, "%s", orinst_sub);
2270 xvgr_legend(fodt, norsel, (const char**)odtleg, oenv);
2276 foten = xvgropen(opt2fn("-oten", NFILE, fnm),
2277 "Order tensor", "Time (ps)", "", oenv);
2278 snew(otenleg, bOvec ? nex*12 : nex*3);
2279 for (i = 0; i < nex; i++)
2281 for (j = 0; j < 3; j++)
2283 sprintf(buf, "eig%d", j+1);
2284 otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
2288 for (j = 0; j < 9; j++)
2290 sprintf(buf, "vec%d%s", j/3+1, j%3 == 0 ? "x" : (j%3 == 1 ? "y" : "z"));
2291 otenleg[12*i+3+j] = strdup(buf);
2295 xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv);
2300 nbounds = get_bounds(ftp2fn(efTPX, NFILE, fnm), &bounds, &index, &pair, &npairs,
2302 snew(violaver, npairs);
2303 out = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations",
2304 "Time (ps)", "nm", oenv);
2305 xvgr_legend(out, 2, drleg, oenv);
2308 fp_pairs = xvgropen(opt2fn("-pairs", NFILE, fnm), "Pair Distances",
2309 "Time (ps)", "Distance (nm)", oenv);
2310 if (output_env_get_print_xvgr_codes(oenv))
2312 fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
2319 get_dhdl_parms(ftp2fn(efTPX, NFILE, fnm), &ir);
2322 /* Initiate energies and set them to zero */
2331 /* Initiate counters */
2334 bFoundStart = FALSE;
2339 /* This loop searches for the first frame (when -b option is given),
2340 * or when this has been found it reads just one energy frame
2344 bCont = do_enx(fp, &(frame[NEXT]));
2347 timecheck = check_times(frame[NEXT].t);
2350 while (bCont && (timecheck < 0));
2352 if ((timecheck == 0) && bCont)
2354 /* We read a valid frame, so we can use it */
2355 fr = &(frame[NEXT]);
2359 /* The frame contains energies, so update cur */
2362 if (edat.nframes % 1000 == 0)
2364 srenew(edat.step, edat.nframes+1000);
2365 memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
2366 srenew(edat.steps, edat.nframes+1000);
2367 memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
2368 srenew(edat.points, edat.nframes+1000);
2369 memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
2371 for (i = 0; i < nset; i++)
2373 srenew(edat.s[i].ener, edat.nframes+1000);
2374 memset(&(edat.s[i].ener[edat.nframes]), 0,
2375 1000*sizeof(edat.s[i].ener[0]));
2376 srenew(edat.s[i].es, edat.nframes+1000);
2377 memset(&(edat.s[i].es[edat.nframes]), 0,
2378 1000*sizeof(edat.s[i].es[0]));
2383 edat.step[nfr] = fr->step;
2388 /* Initiate the previous step data */
2389 start_step = fr->step;
2391 /* Initiate the energy sums */
2392 edat.steps[nfr] = 1;
2393 edat.points[nfr] = 1;
2394 for (i = 0; i < nset; i++)
2397 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2398 edat.s[i].es[nfr].sum2 = 0;
2405 edat.steps[nfr] = fr->nsteps;
2407 if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2411 edat.points[nfr] = 1;
2412 for (i = 0; i < nset; i++)
2415 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2416 edat.s[i].es[nfr].sum2 = 0;
2422 edat.points[nfr] = fr->nsum;
2423 for (i = 0; i < nset; i++)
2426 edat.s[i].es[nfr].sum = fr->ener[sss].esum;
2427 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2429 edat.npoints += fr->nsum;
2434 /* The interval does not match fr->nsteps:
2435 * can not do exact averages.
2439 edat.nsteps = fr->step - start_step + 1;
2442 for (i = 0; i < nset; i++)
2444 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2448 * Define distance restraint legends. Can only be done after
2449 * the first frame has been read... (Then we know how many there are)
2451 blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL);
2452 if (bDisRe && bDRAll && !leg && blk_disre)
2457 fa = top->idef.il[F_DISRES].iatoms;
2458 ip = top->idef.iparams;
2459 if (blk_disre->nsub != 2 ||
2460 (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2462 gmx_incons("Number of disre sub-blocks not equal to 2");
2465 ndisre = blk_disre->sub[0].nr;
2466 if (ndisre != top->idef.il[F_DISRES].nr/3)
2468 gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2469 ndisre, top->idef.il[F_DISRES].nr/3);
2471 snew(pairleg, ndisre);
2472 for (i = 0; i < ndisre; i++)
2474 snew(pairleg[i], 30);
2477 gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j);
2478 gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k);
2479 sprintf(pairleg[i], "%d %s %d %s (%d)",
2480 resnr_j, anm_j, resnr_k, anm_k,
2481 ip[fa[3*i]].disres.label);
2483 set = select_it(ndisre, pairleg, &nset);
2485 for (i = 0; (i < nset); i++)
2488 sprintf(leg[2*i], "a %s", pairleg[set[i]]);
2489 snew(leg[2*i+1], 32);
2490 sprintf(leg[2*i+1], "i %s", pairleg[set[i]]);
2492 xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv);
2496 * Store energies for analysis afterwards...
2498 if (!bDisRe && !bDHDL && (fr->nre > 0))
2500 if (edat.nframes % 1000 == 0)
2502 srenew(time, edat.nframes+1000);
2504 time[edat.nframes] = fr->t;
2508 * Printing time, only when we do not want to skip frames
2510 if (!skip || teller % skip == 0)
2514 /*******************************************
2515 * D I S T A N C E R E S T R A I N T S
2516 *******************************************/
2520 float *disre_rt = blk_disre->sub[0].fval;
2521 float *disre_rm3tav = blk_disre->sub[1].fval;
2523 double *disre_rt = blk_disre->sub[0].dval;
2524 double *disre_rm3tav = blk_disre->sub[1].dval;
2527 print_time(out, fr->t);
2528 if (violaver == NULL)
2530 snew(violaver, ndisre);
2533 /* Subtract bounds from distances, to calculate violations */
2534 calc_violations(disre_rt, disre_rm3tav,
2535 nbounds, pair, bounds, violaver, &sumt, &sumaver);
2537 fprintf(out, " %8.4f %8.4f\n", sumaver, sumt);
2540 print_time(fp_pairs, fr->t);
2541 for (i = 0; (i < nset); i++)
2544 fprintf(fp_pairs, " %8.4f", mypow(disre_rm3tav[sss], minthird));
2545 fprintf(fp_pairs, " %8.4f", disre_rt[sss]);
2547 fprintf(fp_pairs, "\n");
2554 do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2557 /*******************************************
2559 *******************************************/
2566 /* We skip frames with single points (usually only the first frame),
2567 * since they would result in an average plot with outliers.
2571 print_time(out, fr->t);
2572 print1(out, bDp, fr->ener[set[0]].e);
2573 print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
2574 print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum));
2580 print_time(out, fr->t);
2584 for (i = 0; i < nset; i++)
2586 sum += fr->ener[set[i]].e;
2588 print1(out, bDp, sum/nmol-ezero);
2592 for (i = 0; (i < nset); i++)
2596 print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero);
2600 print1(out, bDp, fr->ener[set[i]].e);
2607 blk = find_block_id_enxframe(fr, enx_i, NULL);
2611 xdr_datatype dt = xdr_datatype_float;
2613 xdr_datatype dt = xdr_datatype_double;
2617 if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2619 gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2622 vals = blk->sub[0].fval;
2624 vals = blk->sub[0].dval;
2627 if (blk->sub[0].nr != (size_t)nor)
2629 gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2633 for (i = 0; i < nor; i++)
2635 orient[i] += vals[i];
2640 for (i = 0; i < nor; i++)
2642 odrms[i] += sqr(vals[i]-oobs[i]);
2647 fprintf(fort, " %10f", fr->t);
2648 for (i = 0; i < norsel; i++)
2650 fprintf(fort, " %g", vals[orsel[i]]);
2652 fprintf(fort, "\n");
2656 fprintf(fodt, " %10f", fr->t);
2657 for (i = 0; i < norsel; i++)
2659 fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]);
2661 fprintf(fodt, "\n");
2665 blk = find_block_id_enxframe(fr, enxORT, NULL);
2669 xdr_datatype dt = xdr_datatype_float;
2671 xdr_datatype dt = xdr_datatype_double;
2675 if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2677 gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2680 vals = blk->sub[0].fval;
2682 vals = blk->sub[0].dval;
2685 if (blk->sub[0].nr != (size_t)(nex*12))
2687 gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2688 blk->sub[0].nr/12, nex);
2690 fprintf(foten, " %10f", fr->t);
2691 for (i = 0; i < nex; i++)
2693 for (j = 0; j < (bOvec ? 12 : 3); j++)
2695 fprintf(foten, " %g", vals[i*12+j]);
2698 fprintf(foten, "\n");
2705 while (bCont && (timecheck == 0));
2707 fprintf(stderr, "\n");
2716 gmx_ffclose(fp_pairs);
2729 out = xvgropen(opt2fn("-ora", NFILE, fnm),
2730 "Average calculated orientations",
2731 "Restraint label", "", oenv);
2732 if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2734 fprintf(out, "%s", orinst_sub);
2736 for (i = 0; i < nor; i++)
2738 fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr);
2744 out = xvgropen(opt2fn("-oda", NFILE, fnm),
2745 "Average restraint deviation",
2746 "Restraint label", "", oenv);
2747 if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2749 fprintf(out, "%s", orinst_sub);
2751 for (i = 0; i < nor; i++)
2753 fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr-oobs[i]);
2759 out = xvgropen(opt2fn("-odr", NFILE, fnm),
2760 "RMS orientation restraint deviations",
2761 "Restraint label", "", oenv);
2762 if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2764 fprintf(out, "%s", orinst_sub);
2766 for (i = 0; i < nor; i++)
2768 fprintf(out, "%5d %g\n", or_label[i], sqrt(odrms[i]/norfr));
2779 analyse_disre(opt2fn("-viol", NFILE, fnm),
2780 teller_disre, violaver, bounds, index, pair, nbounds, oenv);
2786 gmx_ffclose(fp_dhdl);
2787 printf("\n\nWrote %d lambda values with %d samples as ",
2788 dh_lambdas, dh_samples);
2791 printf("%d dH histograms ", dh_hists);
2795 printf("%d dH data blocks ", dh_blocks);
2797 printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2802 gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2808 double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2809 analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
2810 bFee, bSum, opt2parg_bSet("-nmol", npargs, ppa),
2811 bVisco, opt2fn("-vis", NFILE, fnm),
2813 start_step, start_t, frame[cur].step, frame[cur].t,
2815 nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
2819 calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat,
2823 if (opt2bSet("-f2", NFILE, fnm))
2825 fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm),
2826 reftemp, nset, set, leg, &edat, time, oenv);
2830 const char *nxy = "-nxy";
2832 do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2833 do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2834 do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
2835 do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
2836 do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
2837 do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
2838 do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
2839 do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
2840 do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);