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.
46 #include "gromacs/utility/fatalerror.h"
47 #include "gromacs/math/vec.h"
48 #include "gromacs/utility/cstringutil.h"
49 #include "gromacs/utility/smalloc.h"
50 #include "gromacs/fileio/enxio.h"
51 #include "gromacs/commandline/pargs.h"
55 #include "gromacs/fileio/xvgr.h"
57 #include "gromacs/math/units.h"
58 #include "gromacs/fileio/tpxio.h"
59 #include "gromacs/fileio/trxio.h"
61 #include "gromacs/topology/mtop_util.h"
65 static real minthird = -1.0/3.0, minsixth = -1.0/6.0;
92 static double mypow(double x, double y)
104 static int *select_it(int nre, char *nm[], int *nset)
109 gmx_bool bVerbose = TRUE;
111 if ((getenv("GMX_ENER_VERBOSE")) != NULL)
116 fprintf(stderr, "Select the terms you want from the following list\n");
117 fprintf(stderr, "End your selection with 0\n");
121 for (k = 0; (k < nre); )
123 for (j = 0; (j < 4) && (k < nre); j++, k++)
125 fprintf(stderr, " %3d=%14s", k+1, nm[k]);
127 fprintf(stderr, "\n");
134 if (1 != scanf("%d", &n))
136 gmx_fatal(FARGS, "Error reading user input");
138 if ((n > 0) && (n <= nre))
146 for (i = (*nset) = 0; (i < nre); i++)
159 static void chomp(char *buf)
161 int len = strlen(buf);
163 while ((len > 0) && (buf[len-1] == '\n'))
170 static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
173 int n, k, kk, j, i, nmatch, nind, nss;
175 gmx_bool bEOF, bVerbose = TRUE, bLong = FALSE;
176 char *ptr, buf[STRLEN];
177 const char *fm4 = "%3d %-14s";
178 const char *fm2 = "%3d %-34s";
181 if ((getenv("GMX_ENER_VERBOSE")) != NULL)
186 fprintf(stderr, "\n");
187 fprintf(stderr, "Select the terms you want from the following list by\n");
188 fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
189 fprintf(stderr, "End your selection with an empty line or a zero.\n");
190 fprintf(stderr, "-------------------------------------------------------------------\n");
194 for (k = 0; k < nre; k++)
196 newnm[k] = strdup(nm[k].name);
197 /* Insert dashes in all the names */
198 while ((ptr = strchr(newnm[k], ' ')) != NULL)
208 fprintf(stderr, "\n");
211 for (kk = k; kk < k+4; kk++)
213 if (kk < nre && strlen(nm[kk].name) > 14)
221 fprintf(stderr, " ");
225 fprintf(stderr, fm4, k+1, newnm[k]);
234 fprintf(stderr, fm2, k+1, newnm[k]);
245 fprintf(stderr, "\n\n");
251 while (!bEOF && (fgets2(buf, STRLEN-1, stdin)))
253 /* Remove newlines */
259 /* Empty line means end of input */
260 bEOF = (strlen(buf) == 0);
268 /* First try to read an integer */
269 nss = sscanf(ptr, "%d", &nind);
272 /* Zero means end of input */
277 else if ((1 <= nind) && (nind <= nre))
283 fprintf(stderr, "number %d is out of range\n", nind);
288 /* Now try to read a string */
291 for (nind = 0; nind < nre; nind++)
293 if (gmx_strcasecmp(newnm[nind], ptr) == 0)
303 for (nind = 0; nind < nre; nind++)
305 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
313 fprintf(stderr, "String '%s' does not match anything\n", ptr);
318 /* Look for the first space, and remove spaces from there */
319 if ((ptr = strchr(ptr, ' ')) != NULL)
324 while (!bEOF && (ptr && (strlen(ptr) > 0)));
329 for (i = (*nset) = 0; (i < nre); i++)
341 gmx_fatal(FARGS, "No energy terms selected");
344 for (i = 0; (i < nre); i++)
353 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
360 /* all we need is the ir to be able to write the label */
361 read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, &mtop);
364 static void get_orires_parms(const char *topnm,
365 int *nor, int *nex, int **label, real **obs)
376 read_tpx(topnm, &ir, box, &natoms, NULL, NULL, NULL, &mtop);
377 top = gmx_mtop_generate_local_top(&mtop, &ir);
379 ip = top->idef.iparams;
380 iatom = top->idef.il[F_ORIRES].iatoms;
382 /* Count how many distance restraint there are... */
383 nb = top->idef.il[F_ORIRES].nr;
386 gmx_fatal(FARGS, "No orientation restraints in topology!\n");
393 for (i = 0; i < nb; i += 3)
395 (*label)[i/3] = ip[iatom[i]].orires.label;
396 (*obs)[i/3] = ip[iatom[i]].orires.obs;
397 if (ip[iatom[i]].orires.ex >= *nex)
399 *nex = ip[iatom[i]].orires.ex+1;
402 fprintf(stderr, "Found %d orientation restraints with %d experiments",
406 static int get_bounds(const char *topnm,
407 real **bounds, int **index, int **dr_pair, int *npairs,
408 gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir)
411 t_functype *functype;
413 int natoms, i, j, k, type, ftype, natom;
421 read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, mtop);
423 top = gmx_mtop_generate_local_top(mtop, ir);
426 functype = top->idef.functype;
427 ip = top->idef.iparams;
429 /* Count how many distance restraint there are... */
430 nb = top->idef.il[F_DISRES].nr;
433 gmx_fatal(FARGS, "No distance restraints in topology!\n");
436 /* Allocate memory */
441 /* Fill the bound array */
443 for (i = 0; (i < top->idef.ntypes); i++)
446 if (ftype == F_DISRES)
449 label1 = ip[i].disres.label;
450 b[nb] = ip[i].disres.up1;
457 /* Fill the index array */
459 disres = &(top->idef.il[F_DISRES]);
460 iatom = disres->iatoms;
461 for (i = j = k = 0; (i < disres->nr); )
464 ftype = top->idef.functype[type];
465 natom = interaction_function[ftype].nratoms+1;
466 if (label1 != top->idef.iparams[type].disres.label)
469 label1 = top->idef.iparams[type].disres.label;
479 gmx_incons("get_bounds for distance restraints");
488 static void calc_violations(real rt[], real rav3[], int nb, int index[],
489 real bounds[], real *viol, double *st, double *sa)
491 const real sixth = 1.0/6.0;
493 double rsum, rav, sumaver, sumt;
497 for (i = 0; (i < nb); i++)
501 for (j = index[i]; (j < index[i+1]); j++)
505 viol[j] += mypow(rt[j], -3.0);
508 rsum += mypow(rt[j], -6);
510 rsum = max(0.0, mypow(rsum, -sixth)-bounds[i]);
511 rav = max(0.0, mypow(rav, -sixth)-bounds[i]);
520 static void analyse_disre(const char *voutfn, int nframes,
521 real violaver[], real bounds[], int index[],
522 int pair[], int nbounds,
523 const output_env_t oenv)
526 double sum, sumt, sumaver;
529 /* Subtract bounds from distances, to calculate violations */
530 calc_violations(violaver, violaver,
531 nbounds, pair, bounds, NULL, &sumt, &sumaver);
534 fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
536 fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n",
539 vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm",
543 for (i = 0; (i < nbounds); i++)
545 /* Do ensemble averaging */
547 for (j = pair[i]; (j < pair[i+1]); j++)
549 sumaver += sqr(violaver[j]/nframes);
551 sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]);
554 sum = max(sum, sumaver);
555 fprintf(vout, "%10d %10.5e\n", index[i], sumaver);
558 for (j = 0; (j < dr.ndr); j++)
560 fprintf(vout, "%10d %10.5e\n", j, mypow(violaver[j]/nframes, minthird));
565 fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
567 fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
569 do_view(oenv, voutfn, "-graphtype bar");
572 static void einstein_visco(const char *fn, const char *fni, int nsets,
573 int nint, real **eneint,
574 real V, real T, double dt,
575 const output_env_t oenv)
578 double av[4], avold[4];
584 for (i = 0; i <= nsets; i++)
588 fp0 = xvgropen(fni, "Shear viscosity integral",
589 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
590 fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
591 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
592 for (i = 0; i < nf4; i++)
594 for (m = 0; m <= nsets; m++)
598 for (j = 0; j < nint - i; j++)
600 for (m = 0; m < nsets; m++)
602 di = sqr(eneint[m][j+i] - eneint[m][j]);
605 av[nsets] += di/nsets;
608 /* Convert to SI for the viscosity */
609 fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nint - i);
610 fprintf(fp0, "%10g", i*dt);
611 for (m = 0; (m <= nsets); m++)
614 fprintf(fp0, " %10g", av[m]);
617 fprintf(fp1, "%10g", (i + 0.5)*dt);
618 for (m = 0; (m <= nsets); m++)
620 fprintf(fp1, " %10g", (av[m]-avold[m])/dt);
643 static void clear_ee_sum(ee_sum_t *ees)
651 static void add_ee_sum(ee_sum_t *ees, double sum, int np)
657 static void add_ee_av(ee_sum_t *ees)
661 av = ees->sum/ees->np;
668 static double calc_ee2(int nb, ee_sum_t *ees)
670 return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
673 static void set_ee_av(ener_ee_t *eee)
677 char buf[STEPSTRSIZE];
678 fprintf(debug, "Storing average for err.est.: %s steps\n",
679 gmx_step_str(eee->nst, buf));
681 add_ee_av(&eee->sum);
683 if (eee->b == 1 || eee->nst < eee->nst_min)
685 eee->nst_min = eee->nst;
690 static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
693 double sum, sum2, sump, see2;
694 gmx_int64_t steps, np, p, bound_nb;
698 double x, sx, sy, sxx, sxy;
701 /* Check if we have exact statistics over all points */
702 for (i = 0; i < nset; i++)
705 ed->bExactStat = FALSE;
706 if (edat->npoints > 0)
708 /* All energy file sum entries 0 signals no exact sums.
709 * But if all energy values are 0, we still have exact sums.
712 for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
714 if (ed->ener[i] != 0)
718 ed->bExactStat = (ed->es[f].sum != 0);
722 ed->bExactStat = TRUE;
728 for (i = 0; i < nset; i++)
739 for (nb = nbmin; nb <= nbmax; nb++)
742 clear_ee_sum(&eee[nb].sum);
746 for (f = 0; f < edat->nframes; f++)
752 /* Add the sum and the sum of variances to the totals. */
758 sum2 += dsqr(sum/np - (sum + es->sum)/(np + p))
764 /* Add a single value to the sum and sum of squares. */
770 /* sum has to be increased after sum2 */
774 /* For the linear regression use variance 1/p.
775 * Note that sump is the sum, not the average, so we don't need p*.
777 x = edat->step[f] - 0.5*(edat->steps[f] - 1);
783 for (nb = nbmin; nb <= nbmax; nb++)
785 /* Check if the current end step is closer to the desired
786 * block boundary than the next end step.
788 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
789 if (eee[nb].nst > 0 &&
790 bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
800 eee[nb].nst += edat->step[f] - edat->step[f-1];
804 add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
808 add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
810 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
811 if (edat->step[f]*nb >= bound_nb)
818 edat->s[i].av = sum/np;
821 edat->s[i].rmsd = sqrt(sum2/np);
825 edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av));
828 if (edat->nframes > 1)
830 edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
834 edat->s[i].slope = 0;
839 for (nb = nbmin; nb <= nbmax; nb++)
841 /* Check if we actually got nb blocks and if the smallest
842 * block is not shorter than 80% of the average.
846 char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
847 fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
849 gmx_step_str(eee[nb].nst_min, buf1),
850 gmx_step_str(edat->nsteps, buf2));
852 if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
854 see2 += calc_ee2(nb, &eee[nb].sum);
860 edat->s[i].ee = sqrt(see2/nee);
870 static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax)
881 snew(s->ener, esum->nframes);
882 snew(s->es, esum->nframes);
884 s->bExactStat = TRUE;
886 for (i = 0; i < nset; i++)
888 if (!edat->s[i].bExactStat)
890 s->bExactStat = FALSE;
892 s->slope += edat->s[i].slope;
895 for (f = 0; f < edat->nframes; f++)
898 for (i = 0; i < nset; i++)
900 sum += edat->s[i].ener[f];
904 for (i = 0; i < nset; i++)
906 sum += edat->s[i].es[f].sum;
912 calc_averages(1, esum, nbmin, nbmax);
917 static char *ee_pr(double ee, char *buf)
924 sprintf(buf, "%s", "--");
928 /* Round to two decimals by printing. */
929 sprintf(tmp, "%.1e", ee);
930 sscanf(tmp, "%lf", &rnd);
931 sprintf(buf, "%g", rnd);
937 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat)
939 /* Remove the drift by performing a fit to y = ax+b.
940 Uses 5 iterations. */
942 double delta, d, sd, sd2;
944 edat->npoints = edat->nframes;
945 edat->nsteps = edat->nframes;
947 for (k = 0; (k < 5); k++)
949 for (i = 0; (i < nset); i++)
951 delta = edat->s[i].slope*dt;
955 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
958 for (j = 0; (j < edat->nframes); j++)
960 edat->s[i].ener[j] -= j*delta;
961 edat->s[i].es[j].sum = 0;
962 edat->s[i].es[j].sum2 = 0;
965 calc_averages(nset, edat, nbmin, nbmax);
969 static void calc_fluctuation_props(FILE *fp,
970 gmx_bool bDriftCorr, real dt,
972 char **leg, enerdata_t *edat,
973 int nbmin, int nbmax)
976 double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
979 eVol, eEnth, eTemp, eEtot, eNR
981 const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
984 NANO3 = NANO*NANO*NANO;
987 fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n"
988 "for spurious drift in the graphs. Note that this is not\n"
989 "a substitute for proper equilibration and sampling!\n");
993 remove_drift(nset, nbmin, nbmax, dt, edat);
995 for (i = 0; (i < eNR); i++)
997 for (ii[i] = 0; (ii[i] < nset &&
998 (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++)
1002 /* if (ii[i] < nset)
1003 fprintf(fp,"Found %s data.\n",my_ener[i]);
1005 /* Compute it all! */
1006 alpha = kappa = cp = dcp = cv = NOTSET;
1010 if (ii[eTemp] < nset)
1012 tt = edat->s[ii[eTemp]].av;
1016 if ((ii[eVol] < nset) && (ii[eTemp] < nset))
1018 vv = edat->s[ii[eVol]].av*NANO3;
1019 varv = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
1020 kappa = (varv/vv)/(BOLTZMANN*tt);
1024 if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
1026 hh = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
1027 varh = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
1028 cp = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
1031 et = varet = NOTSET;
1032 if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
1034 /* Only compute cv in constant volume runs, which we can test
1035 by checking whether the enthalpy was computed.
1037 et = edat->s[ii[eEtot]].av;
1038 varet = sqr(edat->s[ii[eEtot]].rmsd);
1039 cv = KILO*((varet/nmol)/(BOLTZ*tt*tt));
1042 if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
1044 double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
1045 vh_sum = v_sum = h_sum = 0;
1046 for (j = 0; (j < edat->nframes); j++)
1048 v = edat->s[ii[eVol]].ener[j]*NANO3;
1049 h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
1054 vh_aver = vh_sum / edat->nframes;
1055 v_aver = v_sum / edat->nframes;
1056 h_aver = h_sum / edat->nframes;
1057 alpha = (vh_aver-v_aver*h_aver)/(v_aver*BOLTZMANN*tt*tt);
1058 dcp = (v_aver*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
1065 fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n",
1068 fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
1069 fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
1070 fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
1071 fprintf(fp, "please use the g_dos program.\n\n");
1072 fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n"
1073 "a block-averaging error analysis (not implemented in g_energy yet)\n");
1079 fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol);
1084 fprintf(fp, "Volume = %10g m^3/mol\n",
1089 fprintf(fp, "Enthalpy = %10g kJ/mol\n",
1090 hh*AVOGADRO/(KILO*nmol));
1092 if (alpha != NOTSET)
1094 fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
1097 if (kappa != NOTSET)
1099 fprintf(fp, "Isothermal Compressibility Kappa = %10g (J/m^3)\n",
1101 fprintf(fp, "Adiabatic bulk modulus = %10g (m^3/J)\n",
1106 fprintf(fp, "Heat capacity at constant pressure Cp = %10g J/mol K\n",
1111 fprintf(fp, "Heat capacity at constant volume Cv = %10g J/mol K\n",
1116 fprintf(fp, "Cp-Cv = %10g J/mol K\n",
1119 please_cite(fp, "Allen1987a");
1123 fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
1127 static void analyse_ener(gmx_bool bCorr, const char *corrfn,
1128 gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct,
1129 gmx_bool bVisco, const char *visfn, int nmol,
1130 gmx_int64_t start_step, double start_t,
1131 gmx_int64_t step, double t,
1134 int nset, int set[], gmx_bool *bIsEner,
1135 char **leg, gmx_enxnm_t *enm,
1136 real Vaver, real ezero,
1137 int nbmin, int nbmax,
1138 const output_env_t oenv)
1141 /* Check out the printed manual for equations! */
1142 double Dt, aver, stddev, errest, delta_t, totaldrift;
1143 enerdata_t *esum = NULL;
1144 real xxx, integral, intBulk, Temp = 0, Pres = 0;
1145 real sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest;
1146 double beta = 0, expE, expEtot, *fee = NULL;
1148 int nexact, nnotexact;
1152 char buf[256], eebuf[100];
1154 nsteps = step - start_step + 1;
1157 fprintf(stdout, "Not enough steps (%s) for statistics\n",
1158 gmx_step_str(nsteps, buf));
1162 /* Calculate the time difference */
1163 delta_t = t - start_t;
1165 fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
1166 gmx_step_str(nsteps, buf), start_t, t, nset);
1168 calc_averages(nset, edat, nbmin, nbmax);
1172 esum = calc_sum(nset, edat, nbmin, nbmax);
1175 if (edat->npoints == 0)
1184 for (i = 0; (i < nset); i++)
1186 if (edat->s[i].bExactStat)
1199 fprintf(stdout, "All statistics are over %s points\n",
1200 gmx_step_str(edat->npoints, buf));
1202 else if (nexact == 0 || edat->npoints == edat->nframes)
1204 fprintf(stdout, "All statistics are over %d points (frames)\n",
1209 fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
1210 for (i = 0; (i < nset); i++)
1212 if (!edat->s[i].bExactStat)
1214 fprintf(stdout, " '%s'", leg[i]);
1217 fprintf(stdout, " %s has statistics over %d points (frames)\n",
1218 nnotexact == 1 ? "is" : "are", edat->nframes);
1219 fprintf(stdout, "All other statistics are over %s points\n",
1220 gmx_step_str(edat->npoints, buf));
1222 fprintf(stdout, "\n");
1224 fprintf(stdout, "%-24s %10s %10s %10s %10s",
1225 "Energy", "Average", "Err.Est.", "RMSD", "Tot-Drift");
1228 fprintf(stdout, " %10s\n", "-kT ln<e^(E/kT)>");
1232 fprintf(stdout, "\n");
1234 fprintf(stdout, "-------------------------------------------------------------------------------\n");
1236 /* Initiate locals, only used with -sum */
1240 beta = 1.0/(BOLTZ*reftemp);
1243 for (i = 0; (i < nset); i++)
1245 aver = edat->s[i].av;
1246 stddev = edat->s[i].rmsd;
1247 errest = edat->s[i].ee;
1252 for (j = 0; (j < edat->nframes); j++)
1254 expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
1258 expEtot += expE/edat->nframes;
1261 fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
1263 if (strstr(leg[i], "empera") != NULL)
1267 else if (strstr(leg[i], "olum") != NULL)
1271 else if (strstr(leg[i], "essure") != NULL)
1277 pr_aver = aver/nmol-ezero;
1278 pr_stddev = stddev/nmol;
1279 pr_errest = errest/nmol;
1288 /* Multiply the slope in steps with the number of steps taken */
1289 totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
1295 fprintf(stdout, "%-24s %10g %10s %10g %10g",
1296 leg[i], pr_aver, ee_pr(pr_errest, eebuf), pr_stddev, totaldrift);
1299 fprintf(stdout, " %10g", fee[i]);
1302 fprintf(stdout, " (%s)\n", enm[set[i]].unit);
1306 for (j = 0; (j < edat->nframes); j++)
1308 edat->s[i].ener[j] -= aver;
1314 totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1315 fprintf(stdout, "%-24s %10g %10s %10s %10g (%s)",
1316 "Total", esum->s[0].av/nmol, ee_pr(esum->s[0].ee/nmol, eebuf),
1317 "--", totaldrift/nmol, enm[set[0]].unit);
1318 /* pr_aver,pr_stddev,a,totaldrift */
1321 fprintf(stdout, " %10g %10g\n",
1322 log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta);
1326 fprintf(stdout, "\n");
1330 /* Do correlation function */
1331 if (edat->nframes > 1)
1333 Dt = delta_t/(edat->nframes - 1);
1341 const char* leg[] = { "Shear", "Bulk" };
1346 /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1348 /* Symmetrise tensor! (and store in first three elements)
1349 * And subtract average pressure!
1352 for (i = 0; i < 12; i++)
1354 snew(eneset[i], edat->nframes);
1356 for (i = 0; (i < edat->nframes); i++)
1358 eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1359 eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1360 eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1361 for (j = 3; j <= 11; j++)
1363 eneset[j][i] = edat->s[j].ener[i];
1365 eneset[11][i] -= Pres;
1368 /* Determine integrals of the off-diagonal pressure elements */
1370 for (i = 0; i < 3; i++)
1372 snew(eneint[i], edat->nframes + 1);
1377 for (i = 0; i < edat->nframes; i++)
1379 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];
1380 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];
1381 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];
1384 einstein_visco("evisco.xvg", "eviscoi.xvg",
1385 3, edat->nframes+1, eneint, Vaver, Temp, Dt, oenv);
1387 for (i = 0; i < 3; i++)
1393 /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1394 /* Do it for shear viscosity */
1395 strcpy(buf, "Shear Viscosity");
1396 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3,
1397 (edat->nframes+1)/2, eneset, Dt,
1398 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1400 /* Now for bulk viscosity */
1401 strcpy(buf, "Bulk Viscosity");
1402 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1,
1403 (edat->nframes+1)/2, &(eneset[11]), Dt,
1404 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1406 factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1407 fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1408 xvgr_legend(fp, asize(leg), leg, oenv);
1410 /* Use trapezium rule for integration */
1413 nout = get_acfnout();
1414 if ((nout < 2) || (nout >= edat->nframes/2))
1416 nout = edat->nframes/2;
1418 for (i = 1; (i < nout); i++)
1420 integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor;
1421 intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1422 fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk);
1426 for (i = 0; i < 12; i++)
1436 strcpy(buf, "Autocorrelation of Energy Fluctuations");
1440 strcpy(buf, "Energy Autocorrelation");
1443 do_autocorr(corrfn, oenv, buf, edat->nframes,
1445 bSum ? &edat->s[nset-1].ener : eneset,
1446 (delta_t/edat->nframes), eacNormal, FALSE);
1452 static void print_time(FILE *fp, double t)
1454 fprintf(fp, "%12.6f", t);
1457 static void print1(FILE *fp, gmx_bool bDp, real e)
1461 fprintf(fp, " %16.12f", e);
1465 fprintf(fp, " %10.6f", e);
1469 static void fec(const char *ene2fn, const char *runavgfn,
1470 real reftemp, int nset, int set[], char *leg[],
1471 enerdata_t *edat, double time[],
1472 const output_env_t oenv)
1474 const char * ravgleg[] = {
1475 "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1476 "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N"
1480 int nre, timecheck, step, nenergy, nenergy2, maxenergy;
1486 gmx_enxnm_t *enm = NULL;
1490 /* read second energy file */
1493 enx = open_enx(ene2fn, "r");
1494 do_enxnms(enx, &(fr->nre), &enm);
1496 snew(eneset2, nset+1);
1502 /* This loop searches for the first frame (when -b option is given),
1503 * or when this has been found it reads just one energy frame
1507 bCont = do_enx(enx, fr);
1511 timecheck = check_times(fr->t);
1515 while (bCont && (timecheck < 0));
1517 /* Store energies for analysis afterwards... */
1518 if ((timecheck == 0) && bCont)
1522 if (nenergy2 >= maxenergy)
1525 for (i = 0; i <= nset; i++)
1527 srenew(eneset2[i], maxenergy);
1530 if (fr->t != time[nenergy2])
1532 fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n",
1533 fr->t, time[nenergy2], gmx_step_str(fr->step, buf));
1535 for (i = 0; i < nset; i++)
1537 eneset2[i][nenergy2] = fr->ener[set[i]].e;
1543 while (bCont && (timecheck == 0));
1546 if (edat->nframes != nenergy2)
1548 fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n",
1549 edat->nframes, nenergy2);
1551 nenergy = min(edat->nframes, nenergy2);
1553 /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1557 fp = xvgropen(runavgfn, "Running average free energy difference",
1558 "Time (" unit_time ")", "\\8D\\4E (" unit_energy ")", oenv);
1559 xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
1561 fprintf(stdout, "\n%-24s %10s\n",
1562 "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
1564 beta = 1.0/(BOLTZ*reftemp);
1565 for (i = 0; i < nset; i++)
1567 if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
1569 fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n",
1570 leg[i], enm[set[i]].name);
1572 for (j = 0; j < nenergy; j++)
1574 dE = eneset2[i][j] - edat->s[i].ener[j];
1575 sum += exp(-dE*beta);
1578 fprintf(fp, "%10g %10g %10g\n",
1579 time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
1582 aver = -BOLTZ*reftemp*log(sum/nenergy);
1583 fprintf(stdout, "%-24s %10g\n", leg[i], aver);
1593 static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl,
1594 const char *filename, gmx_bool bDp,
1595 int *blocks, int *hists, int *samples, int *nlambdas,
1596 const output_env_t oenv)
1598 const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1599 char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
1601 gmx_bool first = FALSE;
1602 int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
1605 double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0;
1606 static int setnr = 0;
1607 double *native_lambda_vec = NULL;
1608 const char **lambda_components = NULL;
1609 int n_lambda_vec = 0;
1610 gmx_bool changing_lambda = FALSE;
1611 int lambda_fep_state;
1613 /* now count the blocks & handle the global dh data */
1614 for (i = 0; i < fr->nblock; i++)
1616 if (fr->block[i].id == enxDHHIST)
1620 else if (fr->block[i].id == enxDH)
1624 else if (fr->block[i].id == enxDHCOLL)
1627 if ( (fr->block[i].nsub < 1) ||
1628 (fr->block[i].sub[0].type != xdr_datatype_double) ||
1629 (fr->block[i].sub[0].nr < 5))
1631 gmx_fatal(FARGS, "Unexpected block data");
1634 /* read the data from the DHCOLL block */
1635 temp = fr->block[i].sub[0].dval[0];
1636 start_time = fr->block[i].sub[0].dval[1];
1637 delta_time = fr->block[i].sub[0].dval[2];
1638 start_lambda = fr->block[i].sub[0].dval[3];
1639 delta_lambda = fr->block[i].sub[0].dval[4];
1640 changing_lambda = (delta_lambda != 0);
1641 if (fr->block[i].nsub > 1)
1643 lambda_fep_state = fr->block[i].sub[1].ival[0];
1644 if (n_lambda_vec == 0)
1646 n_lambda_vec = fr->block[i].sub[1].ival[1];
1650 if (n_lambda_vec != fr->block[i].sub[1].ival[1])
1653 "Unexpected change of basis set in lambda");
1656 if (lambda_components == NULL)
1658 snew(lambda_components, n_lambda_vec);
1660 if (native_lambda_vec == NULL)
1662 snew(native_lambda_vec, n_lambda_vec);
1664 for (j = 0; j < n_lambda_vec; j++)
1666 native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j];
1667 lambda_components[j] =
1668 efpt_singular_names[fr->block[i].sub[1].ival[2+j]];
1674 if (nblock_hist == 0 && nblock_dh == 0)
1676 /* don't do anything */
1679 if (nblock_hist > 0 && nblock_dh > 0)
1681 gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1687 /* we have standard, non-histogram data --
1688 call open_dhdl to open the file */
1689 /* TODO this is an ugly hack that needs to be fixed: this will only
1690 work if the order of data is always the same and if we're
1691 only using the g_energy compiled with the mdrun that produced
1693 *fp_dhdl = open_dhdl(filename, ir, oenv);
1697 sprintf(title, "N(%s)", deltag);
1698 sprintf(label_x, "%s (%s)", deltag, unit_energy);
1699 sprintf(label_y, "Samples");
1700 *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1701 sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
1702 xvgr_subtitle(*fp_dhdl, buf, oenv);
1706 (*hists) += nblock_hist;
1707 (*blocks) += nblock_dh;
1708 (*nlambdas) = nblock_hist+nblock_dh;
1710 /* write the data */
1711 if (nblock_hist > 0)
1713 gmx_int64_t sum = 0;
1715 for (i = 0; i < fr->nblock; i++)
1717 t_enxblock *blk = &(fr->block[i]);
1718 if (blk->id == enxDHHIST)
1720 double foreign_lambda, dx;
1722 int nhist, derivative;
1724 /* check the block types etc. */
1725 if ( (blk->nsub < 2) ||
1726 (blk->sub[0].type != xdr_datatype_double) ||
1727 (blk->sub[1].type != xdr_datatype_int64) ||
1728 (blk->sub[0].nr < 2) ||
1729 (blk->sub[1].nr < 2) )
1731 gmx_fatal(FARGS, "Unexpected block data in file");
1733 foreign_lambda = blk->sub[0].dval[0];
1734 dx = blk->sub[0].dval[1];
1735 nhist = blk->sub[1].lval[0];
1736 derivative = blk->sub[1].lval[1];
1737 for (j = 0; j < nhist; j++)
1740 x0 = blk->sub[1].lval[2+j];
1744 sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1745 deltag, lambda, foreign_lambda,
1746 lambda, start_lambda);
1750 sprintf(legend, "N(%s | %s=%g)",
1751 dhdl, lambda, start_lambda);
1755 xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1757 for (k = 0; k < blk->sub[j+2].nr; k++)
1762 hist = blk->sub[j+2].ival[k];
1765 fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist,
1769 /* multiple histogram data blocks in one histogram
1770 mean that the second one is the reverse of the first one:
1771 for dhdl derivatives, it's important to know both the
1772 maximum and minimum values */
1777 (*samples) += (int)(sum/nblock_hist);
1783 char **setnames = NULL;
1784 int nnames = nblock_dh;
1786 for (i = 0; i < fr->nblock; i++)
1788 t_enxblock *blk = &(fr->block[i]);
1789 if (blk->id == enxDH)
1793 len = blk->sub[2].nr;
1797 if (len != blk->sub[2].nr)
1799 gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1806 for (i = 0; i < len; i++)
1808 double time = start_time + delta_time*i;
1810 fprintf(*fp_dhdl, "%.4f ", time);
1812 for (j = 0; j < fr->nblock; j++)
1814 t_enxblock *blk = &(fr->block[j]);
1815 if (blk->id == enxDH)
1818 if (blk->sub[2].type == xdr_datatype_float)
1820 value = blk->sub[2].fval[i];
1824 value = blk->sub[2].dval[i];
1826 /* we need to decide which data type it is based on the count*/
1828 if (j == 1 && ir->bExpanded)
1830 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! */
1836 fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */
1840 fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */
1845 fprintf(*fp_dhdl, "\n");
1851 int gmx_energy(int argc, char *argv[])
1853 const char *desc[] = {
1854 "[THISMODULE] extracts energy components or distance restraint",
1855 "data from an energy file. The user is prompted to interactively",
1856 "select the desired energy terms.[PAR]",
1858 "Average, RMSD, and drift are calculated with full precision from the",
1859 "simulation (see printed manual). Drift is calculated by performing",
1860 "a least-squares fit of the data to a straight line. The reported total drift",
1861 "is the difference of the fit at the first and last point.",
1862 "An error estimate of the average is given based on a block averages",
1863 "over 5 blocks using the full-precision averages. The error estimate",
1864 "can be performed over multiple block lengths with the options",
1865 "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1866 "[BB]Note[bb] that in most cases the energy files contains averages over all",
1867 "MD steps, or over many more points than the number of frames in",
1868 "energy file. This makes the [THISMODULE] statistics output more accurate",
1869 "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1870 "file, the statistics mentioned above are simply over the single, per-frame",
1871 "energy values.[PAR]",
1873 "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1875 "Some fluctuation-dependent properties can be calculated provided",
1876 "the correct energy terms are selected, and that the command line option",
1877 "[TT]-fluct_props[tt] is given. The following properties",
1878 "will be computed:[BR]",
1879 "Property Energy terms needed[BR]",
1880 "---------------------------------------------------[BR]",
1881 "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp [BR]",
1882 "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp [BR]",
1883 "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1884 "Isothermal compressibility: Vol, Temp [BR]",
1885 "Adiabatic bulk modulus: Vol, Temp [BR]",
1886 "---------------------------------------------------[BR]",
1887 "You always need to set the number of molecules [TT]-nmol[tt].",
1888 "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1889 "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]"
1890 "When the [TT]-viol[tt] option is set, the time averaged",
1891 "violations are plotted and the running time-averaged and",
1892 "instantaneous sum of violations are recalculated. Additionally",
1893 "running time-averaged and instantaneous distances between",
1894 "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1896 "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1897 "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1898 "The first two options plot the orientation, the last three the",
1899 "deviations of the orientations from the experimental values.",
1900 "The options that end on an 'a' plot the average over time",
1901 "as a function of restraint. The options that end on a 't'",
1902 "prompt the user for restraint label numbers and plot the data",
1903 "as a function of time. Option [TT]-odr[tt] plots the RMS",
1904 "deviation as a function of restraint.",
1905 "When the run used time or ensemble averaged orientation restraints,",
1906 "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1907 "not ensemble-averaged orientations and deviations instead of",
1908 "the time and ensemble averages.[PAR]",
1910 "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1911 "tensor for each orientation restraint experiment. With option",
1912 "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1914 "Option [TT]-odh[tt] extracts and plots the free energy data",
1915 "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1916 "from the [TT]ener.edr[tt] file.[PAR]",
1918 "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1919 "difference with an ideal gas state: [BR]",
1920 " [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]",
1921 " [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]",
1922 "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1923 "the average is over the ensemble (or time in a trajectory).",
1924 "Note that this is in principle",
1925 "only correct when averaging over the whole (Boltzmann) ensemble",
1926 "and using the potential energy. This also allows for an entropy",
1927 "estimate using:[BR]",
1928 " [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]",
1929 " [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",
1932 "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1933 "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1934 "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1935 "files, and the average is over the ensemble A. The running average",
1936 "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1937 "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1940 static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1941 static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1942 static int skip = 0, nmol = 1, nbmin = 5, nbmax = 5;
1943 static real reftemp = 300.0, ezero = 0;
1945 { "-fee", FALSE, etBOOL, {&bFee},
1946 "Do a free energy estimate" },
1947 { "-fetemp", FALSE, etREAL, {&reftemp},
1948 "Reference temperature for free energy calculation" },
1949 { "-zero", FALSE, etREAL, {&ezero},
1950 "Subtract a zero-point energy" },
1951 { "-sum", FALSE, etBOOL, {&bSum},
1952 "Sum the energy terms selected rather than display them all" },
1953 { "-dp", FALSE, etBOOL, {&bDp},
1954 "Print energies in high precision" },
1955 { "-nbmin", FALSE, etINT, {&nbmin},
1956 "Minimum number of blocks for error estimate" },
1957 { "-nbmax", FALSE, etINT, {&nbmax},
1958 "Maximum number of blocks for error estimate" },
1959 { "-mutot", FALSE, etBOOL, {&bMutot},
1960 "Compute the total dipole moment from the components" },
1961 { "-skip", FALSE, etINT, {&skip},
1962 "Skip number of frames between data points" },
1963 { "-aver", FALSE, etBOOL, {&bPrAll},
1964 "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1965 { "-nmol", FALSE, etINT, {&nmol},
1966 "Number of molecules in your sample: the energies are divided by this number" },
1967 { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1968 "Compute properties based on energy fluctuations, like heat capacity" },
1969 { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1970 "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1971 { "-fluc", FALSE, etBOOL, {&bFluct},
1972 "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1973 { "-orinst", FALSE, etBOOL, {&bOrinst},
1974 "Analyse instantaneous orientation data" },
1975 { "-ovec", FALSE, etBOOL, {&bOvec},
1976 "Also plot the eigenvectors with [TT]-oten[tt]" }
1978 const char * drleg[] = {
1982 static const char *setnm[] = {
1983 "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1984 "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1985 "Volume", "Pressure"
1988 FILE *out = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL;
1989 FILE *fp_dhdl = NULL;
1994 gmx_localtop_t *top = NULL;
1998 gmx_enxnm_t *enm = NULL;
1999 t_enxframe *frame, *fr = NULL;
2001 #define NEXT (1-cur)
2002 int nre, teller, teller_disre, nfr;
2003 gmx_int64_t start_step;
2004 int nor = 0, nex = 0, norfr = 0, enx_i = 0;
2006 real *bounds = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL;
2007 int *index = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL;
2008 int nbounds = 0, npairs;
2009 gmx_bool bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL;
2010 gmx_bool bFoundStart, bCont, bEDR, bVisco;
2011 double sum, sumaver, sumt, ener, dbl;
2012 double *time = NULL;
2014 int *set = NULL, i, j, k, nset, sss;
2015 gmx_bool *bIsEner = NULL;
2016 char **pairleg, **odtleg, **otenleg;
2019 char *anm_j, *anm_k, *resnm_j, *resnm_k;
2020 int resnr_j, resnr_k;
2021 const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
2024 t_enxblock *blk = NULL;
2025 t_enxblock *blk_disre = NULL;
2027 int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
2030 { efEDR, "-f", NULL, ffREAD },
2031 { efEDR, "-f2", NULL, ffOPTRD },
2032 { efTPX, "-s", NULL, ffOPTRD },
2033 { efXVG, "-o", "energy", ffWRITE },
2034 { efXVG, "-viol", "violaver", ffOPTWR },
2035 { efXVG, "-pairs", "pairs", ffOPTWR },
2036 { efXVG, "-ora", "orienta", ffOPTWR },
2037 { efXVG, "-ort", "orientt", ffOPTWR },
2038 { efXVG, "-oda", "orideva", ffOPTWR },
2039 { efXVG, "-odr", "oridevr", ffOPTWR },
2040 { efXVG, "-odt", "oridevt", ffOPTWR },
2041 { efXVG, "-oten", "oriten", ffOPTWR },
2042 { efXVG, "-corr", "enecorr", ffOPTWR },
2043 { efXVG, "-vis", "visco", ffOPTWR },
2044 { efXVG, "-ravg", "runavgdf", ffOPTWR },
2045 { efXVG, "-odh", "dhdl", ffOPTWR }
2047 #define NFILE asize(fnm)
2052 ppa = add_acf_pargs(&npargs, pa);
2053 if (!parse_common_args(&argc, argv,
2054 PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
2055 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
2060 bDRAll = opt2bSet("-pairs", NFILE, fnm);
2061 bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
2062 bORA = opt2bSet("-ora", NFILE, fnm);
2063 bORT = opt2bSet("-ort", NFILE, fnm);
2064 bODA = opt2bSet("-oda", NFILE, fnm);
2065 bODR = opt2bSet("-odr", NFILE, fnm);
2066 bODT = opt2bSet("-odt", NFILE, fnm);
2067 bORIRE = bORA || bORT || bODA || bODR || bODT;
2068 bOTEN = opt2bSet("-oten", NFILE, fnm);
2069 bDHDL = opt2bSet("-odh", NFILE, fnm);
2074 fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
2075 do_enxnms(fp, &nre, &enm);
2079 bVisco = opt2bSet("-vis", NFILE, fnm);
2081 if ((!bDisRe) && (!bDHDL))
2085 nset = asize(setnm);
2087 /* This is nasty code... To extract Pres tensor, Volume and Temperature */
2088 for (j = 0; j < nset; j++)
2090 for (i = 0; i < nre; i++)
2092 if (strstr(enm[i].name, setnm[j]))
2100 if (gmx_strcasecmp(setnm[j], "Volume") == 0)
2102 printf("Enter the box volume (" unit_volume "): ");
2103 if (1 != scanf("%lf", &dbl))
2105 gmx_fatal(FARGS, "Error reading user input");
2111 gmx_fatal(FARGS, "Could not find term %s for viscosity calculation",
2119 set = select_by_name(nre, enm, &nset);
2121 /* Print all the different units once */
2122 sprintf(buf, "(%s)", enm[set[0]].unit);
2123 for (i = 1; i < nset; i++)
2125 for (j = 0; j < i; j++)
2127 if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
2135 strcat(buf, enm[set[i]].unit);
2139 out = xvgropen(opt2fn("-o", NFILE, fnm), "Gromacs Energies", "Time (ps)", buf,
2143 for (i = 0; (i < nset); i++)
2145 leg[i] = enm[set[i]].name;
2149 leg[nset] = strdup("Sum");
2150 xvgr_legend(out, nset+1, (const char**)leg, oenv);
2154 xvgr_legend(out, nset, (const char**)leg, oenv);
2157 snew(bIsEner, nset);
2158 for (i = 0; (i < nset); i++)
2161 for (j = 0; (j <= F_ETOT); j++)
2163 bIsEner[i] = bIsEner[i] ||
2164 (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
2168 if (bPrAll && nset > 1)
2170 gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
2175 if (bORIRE || bOTEN)
2177 get_orires_parms(ftp2fn(efTPX, NFILE, fnm), &nor, &nex, &or_label, &oobs);
2201 fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n");
2202 fprintf(stderr, "End your selection with 0\n");
2209 if (1 != scanf("%d", &(orsel[j])))
2211 gmx_fatal(FARGS, "Error reading user input");
2214 while (orsel[j] > 0);
2217 fprintf(stderr, "Selecting all %d orientation restraints\n", nor);
2220 for (i = 0; i < nor; i++)
2227 /* Build the selection */
2229 for (i = 0; i < j; i++)
2231 for (k = 0; k < nor; k++)
2233 if (or_label[k] == orsel[i])
2242 fprintf(stderr, "Orientation restraint label %d not found\n",
2247 snew(odtleg, norsel);
2248 for (i = 0; i < norsel; i++)
2250 snew(odtleg[i], 256);
2251 sprintf(odtleg[i], "%d", or_label[orsel[i]]);
2255 fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
2256 "Time (ps)", "", oenv);
2259 fprintf(fort, "%s", orinst_sub);
2261 xvgr_legend(fort, norsel, (const char**)odtleg, oenv);
2265 fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
2266 "Orientation restraint deviation",
2267 "Time (ps)", "", oenv);
2270 fprintf(fodt, "%s", orinst_sub);
2272 xvgr_legend(fodt, norsel, (const char**)odtleg, oenv);
2278 foten = xvgropen(opt2fn("-oten", NFILE, fnm),
2279 "Order tensor", "Time (ps)", "", oenv);
2280 snew(otenleg, bOvec ? nex*12 : nex*3);
2281 for (i = 0; i < nex; i++)
2283 for (j = 0; j < 3; j++)
2285 sprintf(buf, "eig%d", j+1);
2286 otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
2290 for (j = 0; j < 9; j++)
2292 sprintf(buf, "vec%d%s", j/3+1, j%3 == 0 ? "x" : (j%3 == 1 ? "y" : "z"));
2293 otenleg[12*i+3+j] = strdup(buf);
2297 xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv);
2302 nbounds = get_bounds(ftp2fn(efTPX, NFILE, fnm), &bounds, &index, &pair, &npairs,
2304 snew(violaver, npairs);
2305 out = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations",
2306 "Time (ps)", "nm", oenv);
2307 xvgr_legend(out, 2, drleg, oenv);
2310 fp_pairs = xvgropen(opt2fn("-pairs", NFILE, fnm), "Pair Distances",
2311 "Time (ps)", "Distance (nm)", oenv);
2312 if (output_env_get_print_xvgr_codes(oenv))
2314 fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
2321 get_dhdl_parms(ftp2fn(efTPX, NFILE, fnm), &ir);
2324 /* Initiate energies and set them to zero */
2333 /* Initiate counters */
2336 bFoundStart = FALSE;
2341 /* This loop searches for the first frame (when -b option is given),
2342 * or when this has been found it reads just one energy frame
2346 bCont = do_enx(fp, &(frame[NEXT]));
2349 timecheck = check_times(frame[NEXT].t);
2352 while (bCont && (timecheck < 0));
2354 if ((timecheck == 0) && bCont)
2356 /* We read a valid frame, so we can use it */
2357 fr = &(frame[NEXT]);
2361 /* The frame contains energies, so update cur */
2364 if (edat.nframes % 1000 == 0)
2366 srenew(edat.step, edat.nframes+1000);
2367 memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
2368 srenew(edat.steps, edat.nframes+1000);
2369 memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
2370 srenew(edat.points, edat.nframes+1000);
2371 memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
2373 for (i = 0; i < nset; i++)
2375 srenew(edat.s[i].ener, edat.nframes+1000);
2376 memset(&(edat.s[i].ener[edat.nframes]), 0,
2377 1000*sizeof(edat.s[i].ener[0]));
2378 srenew(edat.s[i].es, edat.nframes+1000);
2379 memset(&(edat.s[i].es[edat.nframes]), 0,
2380 1000*sizeof(edat.s[i].es[0]));
2385 edat.step[nfr] = fr->step;
2390 /* Initiate the previous step data */
2391 start_step = fr->step;
2393 /* Initiate the energy sums */
2394 edat.steps[nfr] = 1;
2395 edat.points[nfr] = 1;
2396 for (i = 0; i < nset; i++)
2399 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2400 edat.s[i].es[nfr].sum2 = 0;
2407 edat.steps[nfr] = fr->nsteps;
2409 if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2413 edat.points[nfr] = 1;
2414 for (i = 0; i < nset; i++)
2417 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2418 edat.s[i].es[nfr].sum2 = 0;
2424 edat.points[nfr] = fr->nsum;
2425 for (i = 0; i < nset; i++)
2428 edat.s[i].es[nfr].sum = fr->ener[sss].esum;
2429 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2431 edat.npoints += fr->nsum;
2436 /* The interval does not match fr->nsteps:
2437 * can not do exact averages.
2441 edat.nsteps = fr->step - start_step + 1;
2444 for (i = 0; i < nset; i++)
2446 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2450 * Define distance restraint legends. Can only be done after
2451 * the first frame has been read... (Then we know how many there are)
2453 blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL);
2454 if (bDisRe && bDRAll && !leg && blk_disre)
2459 fa = top->idef.il[F_DISRES].iatoms;
2460 ip = top->idef.iparams;
2461 if (blk_disre->nsub != 2 ||
2462 (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2464 gmx_incons("Number of disre sub-blocks not equal to 2");
2467 ndisre = blk_disre->sub[0].nr;
2468 if (ndisre != top->idef.il[F_DISRES].nr/3)
2470 gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2471 ndisre, top->idef.il[F_DISRES].nr/3);
2473 snew(pairleg, ndisre);
2474 for (i = 0; i < ndisre; i++)
2476 snew(pairleg[i], 30);
2479 gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j);
2480 gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k);
2481 sprintf(pairleg[i], "%d %s %d %s (%d)",
2482 resnr_j, anm_j, resnr_k, anm_k,
2483 ip[fa[3*i]].disres.label);
2485 set = select_it(ndisre, pairleg, &nset);
2487 for (i = 0; (i < nset); i++)
2490 sprintf(leg[2*i], "a %s", pairleg[set[i]]);
2491 snew(leg[2*i+1], 32);
2492 sprintf(leg[2*i+1], "i %s", pairleg[set[i]]);
2494 xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv);
2498 * Store energies for analysis afterwards...
2500 if (!bDisRe && !bDHDL && (fr->nre > 0))
2502 if (edat.nframes % 1000 == 0)
2504 srenew(time, edat.nframes+1000);
2506 time[edat.nframes] = fr->t;
2510 * Printing time, only when we do not want to skip frames
2512 if (!skip || teller % skip == 0)
2516 /*******************************************
2517 * D I S T A N C E R E S T R A I N T S
2518 *******************************************/
2522 float *disre_rt = blk_disre->sub[0].fval;
2523 float *disre_rm3tav = blk_disre->sub[1].fval;
2525 double *disre_rt = blk_disre->sub[0].dval;
2526 double *disre_rm3tav = blk_disre->sub[1].dval;
2529 print_time(out, fr->t);
2530 if (violaver == NULL)
2532 snew(violaver, ndisre);
2535 /* Subtract bounds from distances, to calculate violations */
2536 calc_violations(disre_rt, disre_rm3tav,
2537 nbounds, pair, bounds, violaver, &sumt, &sumaver);
2539 fprintf(out, " %8.4f %8.4f\n", sumaver, sumt);
2542 print_time(fp_pairs, fr->t);
2543 for (i = 0; (i < nset); i++)
2546 fprintf(fp_pairs, " %8.4f", mypow(disre_rm3tav[sss], minthird));
2547 fprintf(fp_pairs, " %8.4f", disre_rt[sss]);
2549 fprintf(fp_pairs, "\n");
2556 do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2559 /*******************************************
2561 *******************************************/
2568 /* We skip frames with single points (usually only the first frame),
2569 * since they would result in an average plot with outliers.
2573 print_time(out, fr->t);
2574 print1(out, bDp, fr->ener[set[0]].e);
2575 print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
2576 print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum));
2582 print_time(out, fr->t);
2586 for (i = 0; i < nset; i++)
2588 sum += fr->ener[set[i]].e;
2590 print1(out, bDp, sum/nmol-ezero);
2594 for (i = 0; (i < nset); i++)
2598 print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero);
2602 print1(out, bDp, fr->ener[set[i]].e);
2609 blk = find_block_id_enxframe(fr, enx_i, NULL);
2613 xdr_datatype dt = xdr_datatype_float;
2615 xdr_datatype dt = xdr_datatype_double;
2619 if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2621 gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2624 vals = blk->sub[0].fval;
2626 vals = blk->sub[0].dval;
2629 if (blk->sub[0].nr != (size_t)nor)
2631 gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2635 for (i = 0; i < nor; i++)
2637 orient[i] += vals[i];
2642 for (i = 0; i < nor; i++)
2644 odrms[i] += sqr(vals[i]-oobs[i]);
2649 fprintf(fort, " %10f", fr->t);
2650 for (i = 0; i < norsel; i++)
2652 fprintf(fort, " %g", vals[orsel[i]]);
2654 fprintf(fort, "\n");
2658 fprintf(fodt, " %10f", fr->t);
2659 for (i = 0; i < norsel; i++)
2661 fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]);
2663 fprintf(fodt, "\n");
2667 blk = find_block_id_enxframe(fr, enxORT, NULL);
2671 xdr_datatype dt = xdr_datatype_float;
2673 xdr_datatype dt = xdr_datatype_double;
2677 if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2679 gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2682 vals = blk->sub[0].fval;
2684 vals = blk->sub[0].dval;
2687 if (blk->sub[0].nr != (size_t)(nex*12))
2689 gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2690 blk->sub[0].nr/12, nex);
2692 fprintf(foten, " %10f", fr->t);
2693 for (i = 0; i < nex; i++)
2695 for (j = 0; j < (bOvec ? 12 : 3); j++)
2697 fprintf(foten, " %g", vals[i*12+j]);
2700 fprintf(foten, "\n");
2707 while (bCont && (timecheck == 0));
2709 fprintf(stderr, "\n");
2718 gmx_ffclose(fp_pairs);
2731 out = xvgropen(opt2fn("-ora", NFILE, fnm),
2732 "Average calculated orientations",
2733 "Restraint label", "", oenv);
2736 fprintf(out, "%s", orinst_sub);
2738 for (i = 0; i < nor; i++)
2740 fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr);
2746 out = xvgropen(opt2fn("-oda", NFILE, fnm),
2747 "Average restraint deviation",
2748 "Restraint label", "", oenv);
2751 fprintf(out, "%s", orinst_sub);
2753 for (i = 0; i < nor; i++)
2755 fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr-oobs[i]);
2761 out = xvgropen(opt2fn("-odr", NFILE, fnm),
2762 "RMS orientation restraint deviations",
2763 "Restraint label", "", oenv);
2766 fprintf(out, "%s", orinst_sub);
2768 for (i = 0; i < nor; i++)
2770 fprintf(out, "%5d %g\n", or_label[i], sqrt(odrms[i]/norfr));
2781 analyse_disre(opt2fn("-viol", NFILE, fnm),
2782 teller_disre, violaver, bounds, index, pair, nbounds, oenv);
2788 gmx_ffclose(fp_dhdl);
2789 printf("\n\nWrote %d lambda values with %d samples as ",
2790 dh_lambdas, dh_samples);
2793 printf("%d dH histograms ", dh_hists);
2797 printf("%d dH data blocks ", dh_blocks);
2799 printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2804 gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2810 double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2811 analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
2812 bFee, bSum, opt2parg_bSet("-nmol", npargs, ppa),
2813 bVisco, opt2fn("-vis", NFILE, fnm),
2815 start_step, start_t, frame[cur].step, frame[cur].t,
2817 nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
2821 calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat,
2825 if (opt2bSet("-f2", NFILE, fnm))
2827 fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm),
2828 reftemp, nset, set, leg, &edat, time, oenv);
2832 const char *nxy = "-nxy";
2834 do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2835 do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2836 do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
2837 do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
2838 do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
2839 do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
2840 do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
2841 do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
2842 do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);