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.
45 #include "gmx_fatal.h"
47 #include "gromacs/utility/cstringutil.h"
48 #include "gromacs/utility/smalloc.h"
49 #include "gromacs/fileio/enxio.h"
50 #include "gromacs/commandline/pargs.h"
57 #include "gromacs/fileio/tpxio.h"
58 #include "gromacs/fileio/trxio.h"
60 #include "mtop_util.h"
64 static real minthird = -1.0/3.0, minsixth = -1.0/6.0;
91 static double mypow(double x, double y)
103 static int *select_it(int nre, char *nm[], int *nset)
108 gmx_bool bVerbose = TRUE;
110 if ((getenv("GMX_ENER_VERBOSE")) != NULL)
115 fprintf(stderr, "Select the terms you want from the following list\n");
116 fprintf(stderr, "End your selection with 0\n");
120 for (k = 0; (k < nre); )
122 for (j = 0; (j < 4) && (k < nre); j++, k++)
124 fprintf(stderr, " %3d=%14s", k+1, nm[k]);
126 fprintf(stderr, "\n");
133 if (1 != scanf("%d", &n))
135 gmx_fatal(FARGS, "Error reading user input");
137 if ((n > 0) && (n <= nre))
145 for (i = (*nset) = 0; (i < nre); i++)
158 static void chomp(char *buf)
160 int len = strlen(buf);
162 while ((len > 0) && (buf[len-1] == '\n'))
169 static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
172 int n, k, kk, j, i, nmatch, nind, nss;
174 gmx_bool bEOF, bVerbose = TRUE, bLong = FALSE;
175 char *ptr, buf[STRLEN];
176 const char *fm4 = "%3d %-14s";
177 const char *fm2 = "%3d %-34s";
180 if ((getenv("GMX_ENER_VERBOSE")) != NULL)
185 fprintf(stderr, "\n");
186 fprintf(stderr, "Select the terms you want from the following list by\n");
187 fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
188 fprintf(stderr, "End your selection with an empty line or a zero.\n");
189 fprintf(stderr, "-------------------------------------------------------------------\n");
193 for (k = 0; k < nre; k++)
195 newnm[k] = strdup(nm[k].name);
196 /* Insert dashes in all the names */
197 while ((ptr = strchr(newnm[k], ' ')) != NULL)
207 fprintf(stderr, "\n");
210 for (kk = k; kk < k+4; kk++)
212 if (kk < nre && strlen(nm[kk].name) > 14)
220 fprintf(stderr, " ");
224 fprintf(stderr, fm4, k+1, newnm[k]);
233 fprintf(stderr, fm2, k+1, newnm[k]);
244 fprintf(stderr, "\n\n");
250 while (!bEOF && (fgets2(buf, STRLEN-1, stdin)))
252 /* Remove newlines */
258 /* Empty line means end of input */
259 bEOF = (strlen(buf) == 0);
267 /* First try to read an integer */
268 nss = sscanf(ptr, "%d", &nind);
271 /* Zero means end of input */
276 else if ((1 <= nind) && (nind <= nre))
282 fprintf(stderr, "number %d is out of range\n", nind);
287 /* Now try to read a string */
290 for (nind = 0; nind < nre; nind++)
292 if (gmx_strcasecmp(newnm[nind], ptr) == 0)
302 for (nind = 0; nind < nre; nind++)
304 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
312 fprintf(stderr, "String '%s' does not match anything\n", ptr);
317 /* Look for the first space, and remove spaces from there */
318 if ((ptr = strchr(ptr, ' ')) != NULL)
323 while (!bEOF && (ptr && (strlen(ptr) > 0)));
328 for (i = (*nset) = 0; (i < nre); i++)
340 gmx_fatal(FARGS, "No energy terms selected");
343 for (i = 0; (i < nre); i++)
352 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
359 /* all we need is the ir to be able to write the label */
360 read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, &mtop);
363 static void get_orires_parms(const char *topnm,
364 int *nor, int *nex, int **label, real **obs)
375 read_tpx(topnm, &ir, box, &natoms, NULL, NULL, NULL, &mtop);
376 top = gmx_mtop_generate_local_top(&mtop, &ir);
378 ip = top->idef.iparams;
379 iatom = top->idef.il[F_ORIRES].iatoms;
381 /* Count how many distance restraint there are... */
382 nb = top->idef.il[F_ORIRES].nr;
385 gmx_fatal(FARGS, "No orientation restraints in topology!\n");
392 for (i = 0; i < nb; i += 3)
394 (*label)[i/3] = ip[iatom[i]].orires.label;
395 (*obs)[i/3] = ip[iatom[i]].orires.obs;
396 if (ip[iatom[i]].orires.ex >= *nex)
398 *nex = ip[iatom[i]].orires.ex+1;
401 fprintf(stderr, "Found %d orientation restraints with %d experiments",
405 static int get_bounds(const char *topnm,
406 real **bounds, int **index, int **dr_pair, int *npairs,
407 gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir)
410 t_functype *functype;
412 int natoms, i, j, k, type, ftype, natom;
420 read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, mtop);
422 top = gmx_mtop_generate_local_top(mtop, ir);
425 functype = top->idef.functype;
426 ip = top->idef.iparams;
428 /* Count how many distance restraint there are... */
429 nb = top->idef.il[F_DISRES].nr;
432 gmx_fatal(FARGS, "No distance restraints in topology!\n");
435 /* Allocate memory */
440 /* Fill the bound array */
442 for (i = 0; (i < top->idef.ntypes); i++)
445 if (ftype == F_DISRES)
448 label1 = ip[i].disres.label;
449 b[nb] = ip[i].disres.up1;
456 /* Fill the index array */
458 disres = &(top->idef.il[F_DISRES]);
459 iatom = disres->iatoms;
460 for (i = j = k = 0; (i < disres->nr); )
463 ftype = top->idef.functype[type];
464 natom = interaction_function[ftype].nratoms+1;
465 if (label1 != top->idef.iparams[type].disres.label)
468 label1 = top->idef.iparams[type].disres.label;
478 gmx_incons("get_bounds for distance restraints");
487 static void calc_violations(real rt[], real rav3[], int nb, int index[],
488 real bounds[], real *viol, double *st, double *sa)
490 const real sixth = 1.0/6.0;
492 double rsum, rav, sumaver, sumt;
496 for (i = 0; (i < nb); i++)
500 for (j = index[i]; (j < index[i+1]); j++)
504 viol[j] += mypow(rt[j], -3.0);
507 rsum += mypow(rt[j], -6);
509 rsum = max(0.0, mypow(rsum, -sixth)-bounds[i]);
510 rav = max(0.0, mypow(rav, -sixth)-bounds[i]);
519 static void analyse_disre(const char *voutfn, int nframes,
520 real violaver[], real bounds[], int index[],
521 int pair[], int nbounds,
522 const output_env_t oenv)
525 double sum, sumt, sumaver;
528 /* Subtract bounds from distances, to calculate violations */
529 calc_violations(violaver, violaver,
530 nbounds, pair, bounds, NULL, &sumt, &sumaver);
533 fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
535 fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n",
538 vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm",
542 for (i = 0; (i < nbounds); i++)
544 /* Do ensemble averaging */
546 for (j = pair[i]; (j < pair[i+1]); j++)
548 sumaver += sqr(violaver[j]/nframes);
550 sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]);
553 sum = max(sum, sumaver);
554 fprintf(vout, "%10d %10.5e\n", index[i], sumaver);
557 for (j = 0; (j < dr.ndr); j++)
559 fprintf(vout, "%10d %10.5e\n", j, mypow(violaver[j]/nframes, minthird));
564 fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
566 fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
568 do_view(oenv, voutfn, "-graphtype bar");
571 static void einstein_visco(const char *fn, const char *fni, int nsets,
572 int nint, real **eneint,
573 real V, real T, double dt,
574 const output_env_t oenv)
577 double av[4], avold[4];
583 for (i = 0; i <= nsets; i++)
587 fp0 = xvgropen(fni, "Shear viscosity integral",
588 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
589 fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
590 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
591 for (i = 0; i < nf4; i++)
593 for (m = 0; m <= nsets; m++)
597 for (j = 0; j < nint - i; j++)
599 for (m = 0; m < nsets; m++)
601 di = sqr(eneint[m][j+i] - eneint[m][j]);
604 av[nsets] += di/nsets;
607 /* Convert to SI for the viscosity */
608 fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nint - i);
609 fprintf(fp0, "%10g", i*dt);
610 for (m = 0; (m <= nsets); m++)
613 fprintf(fp0, " %10g", av[m]);
616 fprintf(fp1, "%10g", (i + 0.5)*dt);
617 for (m = 0; (m <= nsets); m++)
619 fprintf(fp1, " %10g", (av[m]-avold[m])/dt);
642 static void clear_ee_sum(ee_sum_t *ees)
650 static void add_ee_sum(ee_sum_t *ees, double sum, int np)
656 static void add_ee_av(ee_sum_t *ees)
660 av = ees->sum/ees->np;
667 static double calc_ee2(int nb, ee_sum_t *ees)
669 return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
672 static void set_ee_av(ener_ee_t *eee)
676 char buf[STEPSTRSIZE];
677 fprintf(debug, "Storing average for err.est.: %s steps\n",
678 gmx_step_str(eee->nst, buf));
680 add_ee_av(&eee->sum);
682 if (eee->b == 1 || eee->nst < eee->nst_min)
684 eee->nst_min = eee->nst;
689 static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
692 double sum, sum2, sump, see2;
693 gmx_int64_t steps, np, p, bound_nb;
697 double x, sx, sy, sxx, sxy;
700 /* Check if we have exact statistics over all points */
701 for (i = 0; i < nset; i++)
704 ed->bExactStat = FALSE;
705 if (edat->npoints > 0)
707 /* All energy file sum entries 0 signals no exact sums.
708 * But if all energy values are 0, we still have exact sums.
711 for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
713 if (ed->ener[i] != 0)
717 ed->bExactStat = (ed->es[f].sum != 0);
721 ed->bExactStat = TRUE;
727 for (i = 0; i < nset; i++)
738 for (nb = nbmin; nb <= nbmax; nb++)
741 clear_ee_sum(&eee[nb].sum);
745 for (f = 0; f < edat->nframes; f++)
751 /* Add the sum and the sum of variances to the totals. */
757 sum2 += dsqr(sum/np - (sum + es->sum)/(np + p))
763 /* Add a single value to the sum and sum of squares. */
769 /* sum has to be increased after sum2 */
773 /* For the linear regression use variance 1/p.
774 * Note that sump is the sum, not the average, so we don't need p*.
776 x = edat->step[f] - 0.5*(edat->steps[f] - 1);
782 for (nb = nbmin; nb <= nbmax; nb++)
784 /* Check if the current end step is closer to the desired
785 * block boundary than the next end step.
787 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
788 if (eee[nb].nst > 0 &&
789 bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
799 eee[nb].nst += edat->step[f] - edat->step[f-1];
803 add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
807 add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
809 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
810 if (edat->step[f]*nb >= bound_nb)
817 edat->s[i].av = sum/np;
820 edat->s[i].rmsd = sqrt(sum2/np);
824 edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av));
827 if (edat->nframes > 1)
829 edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
833 edat->s[i].slope = 0;
838 for (nb = nbmin; nb <= nbmax; nb++)
840 /* Check if we actually got nb blocks and if the smallest
841 * block is not shorter than 80% of the average.
845 char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
846 fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
848 gmx_step_str(eee[nb].nst_min, buf1),
849 gmx_step_str(edat->nsteps, buf2));
851 if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
853 see2 += calc_ee2(nb, &eee[nb].sum);
859 edat->s[i].ee = sqrt(see2/nee);
869 static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax)
880 snew(s->ener, esum->nframes);
881 snew(s->es, esum->nframes);
883 s->bExactStat = TRUE;
885 for (i = 0; i < nset; i++)
887 if (!edat->s[i].bExactStat)
889 s->bExactStat = FALSE;
891 s->slope += edat->s[i].slope;
894 for (f = 0; f < edat->nframes; f++)
897 for (i = 0; i < nset; i++)
899 sum += edat->s[i].ener[f];
903 for (i = 0; i < nset; i++)
905 sum += edat->s[i].es[f].sum;
911 calc_averages(1, esum, nbmin, nbmax);
916 static char *ee_pr(double ee, char *buf)
923 sprintf(buf, "%s", "--");
927 /* Round to two decimals by printing. */
928 sprintf(tmp, "%.1e", ee);
929 sscanf(tmp, "%lf", &rnd);
930 sprintf(buf, "%g", rnd);
936 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat)
938 /* Remove the drift by performing a fit to y = ax+b.
939 Uses 5 iterations. */
941 double delta, d, sd, sd2;
943 edat->npoints = edat->nframes;
944 edat->nsteps = edat->nframes;
946 for (k = 0; (k < 5); k++)
948 for (i = 0; (i < nset); i++)
950 delta = edat->s[i].slope*dt;
954 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
957 for (j = 0; (j < edat->nframes); j++)
959 edat->s[i].ener[j] -= j*delta;
960 edat->s[i].es[j].sum = 0;
961 edat->s[i].es[j].sum2 = 0;
964 calc_averages(nset, edat, nbmin, nbmax);
968 static void calc_fluctuation_props(FILE *fp,
969 gmx_bool bDriftCorr, real dt,
971 char **leg, enerdata_t *edat,
972 int nbmin, int nbmax)
975 double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
978 eVol, eEnth, eTemp, eEtot, eNR
980 const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
983 NANO3 = NANO*NANO*NANO;
986 fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n"
987 "for spurious drift in the graphs. Note that this is not\n"
988 "a substitute for proper equilibration and sampling!\n");
992 remove_drift(nset, nbmin, nbmax, dt, edat);
994 for (i = 0; (i < eNR); i++)
996 for (ii[i] = 0; (ii[i] < nset &&
997 (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++)
1001 /* if (ii[i] < nset)
1002 fprintf(fp,"Found %s data.\n",my_ener[i]);
1004 /* Compute it all! */
1005 alpha = kappa = cp = dcp = cv = NOTSET;
1009 if (ii[eTemp] < nset)
1011 tt = edat->s[ii[eTemp]].av;
1015 if ((ii[eVol] < nset) && (ii[eTemp] < nset))
1017 vv = edat->s[ii[eVol]].av*NANO3;
1018 varv = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
1019 kappa = (varv/vv)/(BOLTZMANN*tt);
1023 if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
1025 hh = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
1026 varh = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
1027 cp = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
1030 et = varet = NOTSET;
1031 if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
1033 /* Only compute cv in constant volume runs, which we can test
1034 by checking whether the enthalpy was computed.
1036 et = edat->s[ii[eEtot]].av;
1037 varet = sqr(edat->s[ii[eEtot]].rmsd);
1038 cv = KILO*((varet/nmol)/(BOLTZ*tt*tt));
1041 if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
1043 double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
1044 vh_sum = v_sum = h_sum = 0;
1045 for (j = 0; (j < edat->nframes); j++)
1047 v = edat->s[ii[eVol]].ener[j]*NANO3;
1048 h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
1053 vh_aver = vh_sum / edat->nframes;
1054 v_aver = v_sum / edat->nframes;
1055 h_aver = h_sum / edat->nframes;
1056 alpha = (vh_aver-v_aver*h_aver)/(v_aver*BOLTZMANN*tt*tt);
1057 dcp = (v_aver*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
1064 fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n",
1067 fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
1068 fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
1069 fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
1070 fprintf(fp, "please use the g_dos program.\n\n");
1071 fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n"
1072 "a block-averaging error analysis (not implemented in g_energy yet)\n");
1078 fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol);
1083 fprintf(fp, "Volume = %10g m^3/mol\n",
1088 fprintf(fp, "Enthalpy = %10g kJ/mol\n",
1089 hh*AVOGADRO/(KILO*nmol));
1091 if (alpha != NOTSET)
1093 fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
1096 if (kappa != NOTSET)
1098 fprintf(fp, "Isothermal Compressibility Kappa = %10g (J/m^3)\n",
1100 fprintf(fp, "Adiabatic bulk modulus = %10g (m^3/J)\n",
1105 fprintf(fp, "Heat capacity at constant pressure Cp = %10g J/mol K\n",
1110 fprintf(fp, "Heat capacity at constant volume Cv = %10g J/mol K\n",
1115 fprintf(fp, "Cp-Cv = %10g J/mol K\n",
1118 please_cite(fp, "Allen1987a");
1122 fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
1126 static void analyse_ener(gmx_bool bCorr, const char *corrfn,
1127 gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct,
1128 gmx_bool bVisco, const char *visfn, int nmol,
1129 gmx_int64_t start_step, double start_t,
1130 gmx_int64_t step, double t,
1133 int nset, int set[], gmx_bool *bIsEner,
1134 char **leg, gmx_enxnm_t *enm,
1135 real Vaver, real ezero,
1136 int nbmin, int nbmax,
1137 const output_env_t oenv)
1140 /* Check out the printed manual for equations! */
1141 double Dt, aver, stddev, errest, delta_t, totaldrift;
1142 enerdata_t *esum = NULL;
1143 real xxx, integral, intBulk, Temp = 0, Pres = 0;
1144 real sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest;
1145 double beta = 0, expE, expEtot, *fee = NULL;
1147 int nexact, nnotexact;
1151 char buf[256], eebuf[100];
1153 nsteps = step - start_step + 1;
1156 fprintf(stdout, "Not enough steps (%s) for statistics\n",
1157 gmx_step_str(nsteps, buf));
1161 /* Calculate the time difference */
1162 delta_t = t - start_t;
1164 fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
1165 gmx_step_str(nsteps, buf), start_t, t, nset);
1167 calc_averages(nset, edat, nbmin, nbmax);
1171 esum = calc_sum(nset, edat, nbmin, nbmax);
1174 if (edat->npoints == 0)
1183 for (i = 0; (i < nset); i++)
1185 if (edat->s[i].bExactStat)
1198 fprintf(stdout, "All statistics are over %s points\n",
1199 gmx_step_str(edat->npoints, buf));
1201 else if (nexact == 0 || edat->npoints == edat->nframes)
1203 fprintf(stdout, "All statistics are over %d points (frames)\n",
1208 fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
1209 for (i = 0; (i < nset); i++)
1211 if (!edat->s[i].bExactStat)
1213 fprintf(stdout, " '%s'", leg[i]);
1216 fprintf(stdout, " %s has statistics over %d points (frames)\n",
1217 nnotexact == 1 ? "is" : "are", edat->nframes);
1218 fprintf(stdout, "All other statistics are over %s points\n",
1219 gmx_step_str(edat->npoints, buf));
1221 fprintf(stdout, "\n");
1223 fprintf(stdout, "%-24s %10s %10s %10s %10s",
1224 "Energy", "Average", "Err.Est.", "RMSD", "Tot-Drift");
1227 fprintf(stdout, " %10s\n", "-kT ln<e^(E/kT)>");
1231 fprintf(stdout, "\n");
1233 fprintf(stdout, "-------------------------------------------------------------------------------\n");
1235 /* Initiate locals, only used with -sum */
1239 beta = 1.0/(BOLTZ*reftemp);
1242 for (i = 0; (i < nset); i++)
1244 aver = edat->s[i].av;
1245 stddev = edat->s[i].rmsd;
1246 errest = edat->s[i].ee;
1251 for (j = 0; (j < edat->nframes); j++)
1253 expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
1257 expEtot += expE/edat->nframes;
1260 fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
1262 if (strstr(leg[i], "empera") != NULL)
1266 else if (strstr(leg[i], "olum") != NULL)
1270 else if (strstr(leg[i], "essure") != NULL)
1276 pr_aver = aver/nmol-ezero;
1277 pr_stddev = stddev/nmol;
1278 pr_errest = errest/nmol;
1287 /* Multiply the slope in steps with the number of steps taken */
1288 totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
1294 fprintf(stdout, "%-24s %10g %10s %10g %10g",
1295 leg[i], pr_aver, ee_pr(pr_errest, eebuf), pr_stddev, totaldrift);
1298 fprintf(stdout, " %10g", fee[i]);
1301 fprintf(stdout, " (%s)\n", enm[set[i]].unit);
1305 for (j = 0; (j < edat->nframes); j++)
1307 edat->s[i].ener[j] -= aver;
1313 totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1314 fprintf(stdout, "%-24s %10g %10s %10s %10g (%s)",
1315 "Total", esum->s[0].av/nmol, ee_pr(esum->s[0].ee/nmol, eebuf),
1316 "--", totaldrift/nmol, enm[set[0]].unit);
1317 /* pr_aver,pr_stddev,a,totaldrift */
1320 fprintf(stdout, " %10g %10g\n",
1321 log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta);
1325 fprintf(stdout, "\n");
1329 /* Do correlation function */
1330 if (edat->nframes > 1)
1332 Dt = delta_t/(edat->nframes - 1);
1340 const char* leg[] = { "Shear", "Bulk" };
1345 /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1347 /* Symmetrise tensor! (and store in first three elements)
1348 * And subtract average pressure!
1351 for (i = 0; i < 12; i++)
1353 snew(eneset[i], edat->nframes);
1355 for (i = 0; (i < edat->nframes); i++)
1357 eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1358 eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1359 eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1360 for (j = 3; j <= 11; j++)
1362 eneset[j][i] = edat->s[j].ener[i];
1364 eneset[11][i] -= Pres;
1367 /* Determine integrals of the off-diagonal pressure elements */
1369 for (i = 0; i < 3; i++)
1371 snew(eneint[i], edat->nframes + 1);
1376 for (i = 0; i < edat->nframes; i++)
1378 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];
1379 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];
1380 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];
1383 einstein_visco("evisco.xvg", "eviscoi.xvg",
1384 3, edat->nframes+1, eneint, Vaver, Temp, Dt, oenv);
1386 for (i = 0; i < 3; i++)
1392 /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1393 /* Do it for shear viscosity */
1394 strcpy(buf, "Shear Viscosity");
1395 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3,
1396 (edat->nframes+1)/2, eneset, Dt,
1397 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1399 /* Now for bulk viscosity */
1400 strcpy(buf, "Bulk Viscosity");
1401 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1,
1402 (edat->nframes+1)/2, &(eneset[11]), Dt,
1403 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1405 factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1406 fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1407 xvgr_legend(fp, asize(leg), leg, oenv);
1409 /* Use trapezium rule for integration */
1412 nout = get_acfnout();
1413 if ((nout < 2) || (nout >= edat->nframes/2))
1415 nout = edat->nframes/2;
1417 for (i = 1; (i < nout); i++)
1419 integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor;
1420 intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1421 fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk);
1425 for (i = 0; i < 12; i++)
1435 strcpy(buf, "Autocorrelation of Energy Fluctuations");
1439 strcpy(buf, "Energy Autocorrelation");
1442 do_autocorr(corrfn, oenv, buf, edat->nframes,
1444 bSum ? &edat->s[nset-1].ener : eneset,
1445 (delta_t/edat->nframes), eacNormal, FALSE);
1451 static void print_time(FILE *fp, double t)
1453 fprintf(fp, "%12.6f", t);
1456 static void print1(FILE *fp, gmx_bool bDp, real e)
1460 fprintf(fp, " %16.12f", e);
1464 fprintf(fp, " %10.6f", e);
1468 static void fec(const char *ene2fn, const char *runavgfn,
1469 real reftemp, int nset, int set[], char *leg[],
1470 enerdata_t *edat, double time[],
1471 const output_env_t oenv)
1473 const char * ravgleg[] = {
1474 "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1475 "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N"
1479 int nre, timecheck, step, nenergy, nenergy2, maxenergy;
1485 gmx_enxnm_t *enm = NULL;
1489 /* read second energy file */
1492 enx = open_enx(ene2fn, "r");
1493 do_enxnms(enx, &(fr->nre), &enm);
1495 snew(eneset2, nset+1);
1501 /* This loop searches for the first frame (when -b option is given),
1502 * or when this has been found it reads just one energy frame
1506 bCont = do_enx(enx, fr);
1510 timecheck = check_times(fr->t);
1514 while (bCont && (timecheck < 0));
1516 /* Store energies for analysis afterwards... */
1517 if ((timecheck == 0) && bCont)
1521 if (nenergy2 >= maxenergy)
1524 for (i = 0; i <= nset; i++)
1526 srenew(eneset2[i], maxenergy);
1529 if (fr->t != time[nenergy2])
1531 fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n",
1532 fr->t, time[nenergy2], gmx_step_str(fr->step, buf));
1534 for (i = 0; i < nset; i++)
1536 eneset2[i][nenergy2] = fr->ener[set[i]].e;
1542 while (bCont && (timecheck == 0));
1545 if (edat->nframes != nenergy2)
1547 fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n",
1548 edat->nframes, nenergy2);
1550 nenergy = min(edat->nframes, nenergy2);
1552 /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1556 fp = xvgropen(runavgfn, "Running average free energy difference",
1557 "Time (" unit_time ")", "\\8D\\4E (" unit_energy ")", oenv);
1558 xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
1560 fprintf(stdout, "\n%-24s %10s\n",
1561 "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
1563 beta = 1.0/(BOLTZ*reftemp);
1564 for (i = 0; i < nset; i++)
1566 if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
1568 fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n",
1569 leg[i], enm[set[i]].name);
1571 for (j = 0; j < nenergy; j++)
1573 dE = eneset2[i][j] - edat->s[i].ener[j];
1574 sum += exp(-dE*beta);
1577 fprintf(fp, "%10g %10g %10g\n",
1578 time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
1581 aver = -BOLTZ*reftemp*log(sum/nenergy);
1582 fprintf(stdout, "%-24s %10g\n", leg[i], aver);
1592 static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl,
1593 const char *filename, gmx_bool bDp,
1594 int *blocks, int *hists, int *samples, int *nlambdas,
1595 const output_env_t oenv)
1597 const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1598 char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
1600 gmx_bool first = FALSE;
1601 int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
1604 double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0;
1605 static int setnr = 0;
1606 double *native_lambda_vec = NULL;
1607 const char **lambda_components = NULL;
1608 int n_lambda_vec = 0;
1609 gmx_bool changing_lambda = FALSE;
1610 int lambda_fep_state;
1612 /* now count the blocks & handle the global dh data */
1613 for (i = 0; i < fr->nblock; i++)
1615 if (fr->block[i].id == enxDHHIST)
1619 else if (fr->block[i].id == enxDH)
1623 else if (fr->block[i].id == enxDHCOLL)
1626 if ( (fr->block[i].nsub < 1) ||
1627 (fr->block[i].sub[0].type != xdr_datatype_double) ||
1628 (fr->block[i].sub[0].nr < 5))
1630 gmx_fatal(FARGS, "Unexpected block data");
1633 /* read the data from the DHCOLL block */
1634 temp = fr->block[i].sub[0].dval[0];
1635 start_time = fr->block[i].sub[0].dval[1];
1636 delta_time = fr->block[i].sub[0].dval[2];
1637 start_lambda = fr->block[i].sub[0].dval[3];
1638 delta_lambda = fr->block[i].sub[0].dval[4];
1639 changing_lambda = (delta_lambda != 0);
1640 if (fr->block[i].nsub > 1)
1642 lambda_fep_state = fr->block[i].sub[1].ival[0];
1643 if (n_lambda_vec == 0)
1645 n_lambda_vec = fr->block[i].sub[1].ival[1];
1649 if (n_lambda_vec != fr->block[i].sub[1].ival[1])
1652 "Unexpected change of basis set in lambda");
1655 if (lambda_components == NULL)
1657 snew(lambda_components, n_lambda_vec);
1659 if (native_lambda_vec == NULL)
1661 snew(native_lambda_vec, n_lambda_vec);
1663 for (j = 0; j < n_lambda_vec; j++)
1665 native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j];
1666 lambda_components[j] =
1667 efpt_singular_names[fr->block[i].sub[1].ival[2+j]];
1673 if (nblock_hist == 0 && nblock_dh == 0)
1675 /* don't do anything */
1678 if (nblock_hist > 0 && nblock_dh > 0)
1680 gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1686 /* we have standard, non-histogram data --
1687 call open_dhdl to open the file */
1688 /* TODO this is an ugly hack that needs to be fixed: this will only
1689 work if the order of data is always the same and if we're
1690 only using the g_energy compiled with the mdrun that produced
1692 *fp_dhdl = open_dhdl(filename, ir, oenv);
1696 sprintf(title, "N(%s)", deltag);
1697 sprintf(label_x, "%s (%s)", deltag, unit_energy);
1698 sprintf(label_y, "Samples");
1699 *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1700 sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
1701 xvgr_subtitle(*fp_dhdl, buf, oenv);
1705 (*hists) += nblock_hist;
1706 (*blocks) += nblock_dh;
1707 (*nlambdas) = nblock_hist+nblock_dh;
1709 /* write the data */
1710 if (nblock_hist > 0)
1712 gmx_int64_t sum = 0;
1714 for (i = 0; i < fr->nblock; i++)
1716 t_enxblock *blk = &(fr->block[i]);
1717 if (blk->id == enxDHHIST)
1719 double foreign_lambda, dx;
1721 int nhist, derivative;
1723 /* check the block types etc. */
1724 if ( (blk->nsub < 2) ||
1725 (blk->sub[0].type != xdr_datatype_double) ||
1726 (blk->sub[1].type != xdr_datatype_int64) ||
1727 (blk->sub[0].nr < 2) ||
1728 (blk->sub[1].nr < 2) )
1730 gmx_fatal(FARGS, "Unexpected block data in file");
1732 foreign_lambda = blk->sub[0].dval[0];
1733 dx = blk->sub[0].dval[1];
1734 nhist = blk->sub[1].lval[0];
1735 derivative = blk->sub[1].lval[1];
1736 for (j = 0; j < nhist; j++)
1739 x0 = blk->sub[1].lval[2+j];
1743 sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1744 deltag, lambda, foreign_lambda,
1745 lambda, start_lambda);
1749 sprintf(legend, "N(%s | %s=%g)",
1750 dhdl, lambda, start_lambda);
1754 xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1756 for (k = 0; k < blk->sub[j+2].nr; k++)
1761 hist = blk->sub[j+2].ival[k];
1764 fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist,
1768 /* multiple histogram data blocks in one histogram
1769 mean that the second one is the reverse of the first one:
1770 for dhdl derivatives, it's important to know both the
1771 maximum and minimum values */
1776 (*samples) += (int)(sum/nblock_hist);
1782 char **setnames = NULL;
1783 int nnames = nblock_dh;
1785 for (i = 0; i < fr->nblock; i++)
1787 t_enxblock *blk = &(fr->block[i]);
1788 if (blk->id == enxDH)
1792 len = blk->sub[2].nr;
1796 if (len != blk->sub[2].nr)
1798 gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1805 for (i = 0; i < len; i++)
1807 double time = start_time + delta_time*i;
1809 fprintf(*fp_dhdl, "%.4f ", time);
1811 for (j = 0; j < fr->nblock; j++)
1813 t_enxblock *blk = &(fr->block[j]);
1814 if (blk->id == enxDH)
1817 if (blk->sub[2].type == xdr_datatype_float)
1819 value = blk->sub[2].fval[i];
1823 value = blk->sub[2].dval[i];
1825 /* we need to decide which data type it is based on the count*/
1827 if (j == 1 && ir->bExpanded)
1829 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! */
1835 fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */
1839 fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */
1844 fprintf(*fp_dhdl, "\n");
1850 int gmx_energy(int argc, char *argv[])
1852 const char *desc[] = {
1853 "[THISMODULE] extracts energy components or distance restraint",
1854 "data from an energy file. The user is prompted to interactively",
1855 "select the desired energy terms.[PAR]",
1857 "Average, RMSD, and drift are calculated with full precision from the",
1858 "simulation (see printed manual). Drift is calculated by performing",
1859 "a least-squares fit of the data to a straight line. The reported total drift",
1860 "is the difference of the fit at the first and last point.",
1861 "An error estimate of the average is given based on a block averages",
1862 "over 5 blocks using the full-precision averages. The error estimate",
1863 "can be performed over multiple block lengths with the options",
1864 "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1865 "[BB]Note[bb] that in most cases the energy files contains averages over all",
1866 "MD steps, or over many more points than the number of frames in",
1867 "energy file. This makes the [THISMODULE] statistics output more accurate",
1868 "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1869 "file, the statistics mentioned above are simply over the single, per-frame",
1870 "energy values.[PAR]",
1872 "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1874 "Some fluctuation-dependent properties can be calculated provided",
1875 "the correct energy terms are selected, and that the command line option",
1876 "[TT]-fluct_props[tt] is given. The following properties",
1877 "will be computed:[BR]",
1878 "Property Energy terms needed[BR]",
1879 "---------------------------------------------------[BR]",
1880 "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp [BR]",
1881 "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp [BR]",
1882 "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1883 "Isothermal compressibility: Vol, Temp [BR]",
1884 "Adiabatic bulk modulus: Vol, Temp [BR]",
1885 "---------------------------------------------------[BR]",
1886 "You always need to set the number of molecules [TT]-nmol[tt].",
1887 "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1888 "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]"
1889 "When the [TT]-viol[tt] option is set, the time averaged",
1890 "violations are plotted and the running time-averaged and",
1891 "instantaneous sum of violations are recalculated. Additionally",
1892 "running time-averaged and instantaneous distances between",
1893 "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1895 "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1896 "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1897 "The first two options plot the orientation, the last three the",
1898 "deviations of the orientations from the experimental values.",
1899 "The options that end on an 'a' plot the average over time",
1900 "as a function of restraint. The options that end on a 't'",
1901 "prompt the user for restraint label numbers and plot the data",
1902 "as a function of time. Option [TT]-odr[tt] plots the RMS",
1903 "deviation as a function of restraint.",
1904 "When the run used time or ensemble averaged orientation restraints,",
1905 "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1906 "not ensemble-averaged orientations and deviations instead of",
1907 "the time and ensemble averages.[PAR]",
1909 "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1910 "tensor for each orientation restraint experiment. With option",
1911 "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1913 "Option [TT]-odh[tt] extracts and plots the free energy data",
1914 "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1915 "from the [TT]ener.edr[tt] file.[PAR]",
1917 "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1918 "difference with an ideal gas state: [BR]",
1919 " [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]",
1920 " [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]",
1921 "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1922 "the average is over the ensemble (or time in a trajectory).",
1923 "Note that this is in principle",
1924 "only correct when averaging over the whole (Boltzmann) ensemble",
1925 "and using the potential energy. This also allows for an entropy",
1926 "estimate using:[BR]",
1927 " [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]",
1928 " [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",
1931 "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1932 "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1933 "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1934 "files, and the average is over the ensemble A. The running average",
1935 "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1936 "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1939 static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1940 static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1941 static int skip = 0, nmol = 1, nbmin = 5, nbmax = 5;
1942 static real reftemp = 300.0, ezero = 0;
1944 { "-fee", FALSE, etBOOL, {&bFee},
1945 "Do a free energy estimate" },
1946 { "-fetemp", FALSE, etREAL, {&reftemp},
1947 "Reference temperature for free energy calculation" },
1948 { "-zero", FALSE, etREAL, {&ezero},
1949 "Subtract a zero-point energy" },
1950 { "-sum", FALSE, etBOOL, {&bSum},
1951 "Sum the energy terms selected rather than display them all" },
1952 { "-dp", FALSE, etBOOL, {&bDp},
1953 "Print energies in high precision" },
1954 { "-nbmin", FALSE, etINT, {&nbmin},
1955 "Minimum number of blocks for error estimate" },
1956 { "-nbmax", FALSE, etINT, {&nbmax},
1957 "Maximum number of blocks for error estimate" },
1958 { "-mutot", FALSE, etBOOL, {&bMutot},
1959 "Compute the total dipole moment from the components" },
1960 { "-skip", FALSE, etINT, {&skip},
1961 "Skip number of frames between data points" },
1962 { "-aver", FALSE, etBOOL, {&bPrAll},
1963 "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1964 { "-nmol", FALSE, etINT, {&nmol},
1965 "Number of molecules in your sample: the energies are divided by this number" },
1966 { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1967 "Compute properties based on energy fluctuations, like heat capacity" },
1968 { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1969 "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1970 { "-fluc", FALSE, etBOOL, {&bFluct},
1971 "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1972 { "-orinst", FALSE, etBOOL, {&bOrinst},
1973 "Analyse instantaneous orientation data" },
1974 { "-ovec", FALSE, etBOOL, {&bOvec},
1975 "Also plot the eigenvectors with [TT]-oten[tt]" }
1977 const char * drleg[] = {
1981 static const char *setnm[] = {
1982 "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1983 "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1984 "Volume", "Pressure"
1987 FILE *out = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL;
1988 FILE *fp_dhdl = NULL;
1993 gmx_localtop_t *top = NULL;
1997 gmx_enxnm_t *enm = NULL;
1998 t_enxframe *frame, *fr = NULL;
2000 #define NEXT (1-cur)
2001 int nre, teller, teller_disre, nfr;
2002 gmx_int64_t start_step;
2003 int nor = 0, nex = 0, norfr = 0, enx_i = 0;
2005 real *bounds = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL;
2006 int *index = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL;
2007 int nbounds = 0, npairs;
2008 gmx_bool bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL;
2009 gmx_bool bFoundStart, bCont, bEDR, bVisco;
2010 double sum, sumaver, sumt, ener, dbl;
2011 double *time = NULL;
2013 int *set = NULL, i, j, k, nset, sss;
2014 gmx_bool *bIsEner = NULL;
2015 char **pairleg, **odtleg, **otenleg;
2018 char *anm_j, *anm_k, *resnm_j, *resnm_k;
2019 int resnr_j, resnr_k;
2020 const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
2023 t_enxblock *blk = NULL;
2024 t_enxblock *blk_disre = NULL;
2026 int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
2029 { efEDR, "-f", NULL, ffREAD },
2030 { efEDR, "-f2", NULL, ffOPTRD },
2031 { efTPX, "-s", NULL, ffOPTRD },
2032 { efXVG, "-o", "energy", ffWRITE },
2033 { efXVG, "-viol", "violaver", ffOPTWR },
2034 { efXVG, "-pairs", "pairs", ffOPTWR },
2035 { efXVG, "-ora", "orienta", ffOPTWR },
2036 { efXVG, "-ort", "orientt", ffOPTWR },
2037 { efXVG, "-oda", "orideva", ffOPTWR },
2038 { efXVG, "-odr", "oridevr", ffOPTWR },
2039 { efXVG, "-odt", "oridevt", ffOPTWR },
2040 { efXVG, "-oten", "oriten", ffOPTWR },
2041 { efXVG, "-corr", "enecorr", ffOPTWR },
2042 { efXVG, "-vis", "visco", ffOPTWR },
2043 { efXVG, "-ravg", "runavgdf", ffOPTWR },
2044 { efXVG, "-odh", "dhdl", ffOPTWR }
2046 #define NFILE asize(fnm)
2051 ppa = add_acf_pargs(&npargs, pa);
2052 if (!parse_common_args(&argc, argv,
2053 PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
2054 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
2059 bDRAll = opt2bSet("-pairs", NFILE, fnm);
2060 bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
2061 bORA = opt2bSet("-ora", NFILE, fnm);
2062 bORT = opt2bSet("-ort", NFILE, fnm);
2063 bODA = opt2bSet("-oda", NFILE, fnm);
2064 bODR = opt2bSet("-odr", NFILE, fnm);
2065 bODT = opt2bSet("-odt", NFILE, fnm);
2066 bORIRE = bORA || bORT || bODA || bODR || bODT;
2067 bOTEN = opt2bSet("-oten", NFILE, fnm);
2068 bDHDL = opt2bSet("-odh", NFILE, fnm);
2073 fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
2074 do_enxnms(fp, &nre, &enm);
2078 bVisco = opt2bSet("-vis", NFILE, fnm);
2080 if ((!bDisRe) && (!bDHDL))
2084 nset = asize(setnm);
2086 /* This is nasty code... To extract Pres tensor, Volume and Temperature */
2087 for (j = 0; j < nset; j++)
2089 for (i = 0; i < nre; i++)
2091 if (strstr(enm[i].name, setnm[j]))
2099 if (gmx_strcasecmp(setnm[j], "Volume") == 0)
2101 printf("Enter the box volume (" unit_volume "): ");
2102 if (1 != scanf("%lf", &dbl))
2104 gmx_fatal(FARGS, "Error reading user input");
2110 gmx_fatal(FARGS, "Could not find term %s for viscosity calculation",
2118 set = select_by_name(nre, enm, &nset);
2120 /* Print all the different units once */
2121 sprintf(buf, "(%s)", enm[set[0]].unit);
2122 for (i = 1; i < nset; i++)
2124 for (j = 0; j < i; j++)
2126 if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
2134 strcat(buf, enm[set[i]].unit);
2138 out = xvgropen(opt2fn("-o", NFILE, fnm), "Gromacs Energies", "Time (ps)", buf,
2142 for (i = 0; (i < nset); i++)
2144 leg[i] = enm[set[i]].name;
2148 leg[nset] = strdup("Sum");
2149 xvgr_legend(out, nset+1, (const char**)leg, oenv);
2153 xvgr_legend(out, nset, (const char**)leg, oenv);
2156 snew(bIsEner, nset);
2157 for (i = 0; (i < nset); i++)
2160 for (j = 0; (j <= F_ETOT); j++)
2162 bIsEner[i] = bIsEner[i] ||
2163 (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
2167 if (bPrAll && nset > 1)
2169 gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
2174 if (bORIRE || bOTEN)
2176 get_orires_parms(ftp2fn(efTPX, NFILE, fnm), &nor, &nex, &or_label, &oobs);
2200 fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n");
2201 fprintf(stderr, "End your selection with 0\n");
2208 if (1 != scanf("%d", &(orsel[j])))
2210 gmx_fatal(FARGS, "Error reading user input");
2213 while (orsel[j] > 0);
2216 fprintf(stderr, "Selecting all %d orientation restraints\n", nor);
2219 for (i = 0; i < nor; i++)
2226 /* Build the selection */
2228 for (i = 0; i < j; i++)
2230 for (k = 0; k < nor; k++)
2232 if (or_label[k] == orsel[i])
2241 fprintf(stderr, "Orientation restraint label %d not found\n",
2246 snew(odtleg, norsel);
2247 for (i = 0; i < norsel; i++)
2249 snew(odtleg[i], 256);
2250 sprintf(odtleg[i], "%d", or_label[orsel[i]]);
2254 fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
2255 "Time (ps)", "", oenv);
2256 if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2258 fprintf(fort, "%s", orinst_sub);
2260 xvgr_legend(fort, norsel, (const char**)odtleg, oenv);
2264 fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
2265 "Orientation restraint deviation",
2266 "Time (ps)", "", oenv);
2267 if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2269 fprintf(fodt, "%s", orinst_sub);
2271 xvgr_legend(fodt, norsel, (const char**)odtleg, oenv);
2277 foten = xvgropen(opt2fn("-oten", NFILE, fnm),
2278 "Order tensor", "Time (ps)", "", oenv);
2279 snew(otenleg, bOvec ? nex*12 : nex*3);
2280 for (i = 0; i < nex; i++)
2282 for (j = 0; j < 3; j++)
2284 sprintf(buf, "eig%d", j+1);
2285 otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
2289 for (j = 0; j < 9; j++)
2291 sprintf(buf, "vec%d%s", j/3+1, j%3 == 0 ? "x" : (j%3 == 1 ? "y" : "z"));
2292 otenleg[12*i+3+j] = strdup(buf);
2296 xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv);
2301 nbounds = get_bounds(ftp2fn(efTPX, NFILE, fnm), &bounds, &index, &pair, &npairs,
2303 snew(violaver, npairs);
2304 out = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations",
2305 "Time (ps)", "nm", oenv);
2306 xvgr_legend(out, 2, drleg, oenv);
2309 fp_pairs = xvgropen(opt2fn("-pairs", NFILE, fnm), "Pair Distances",
2310 "Time (ps)", "Distance (nm)", oenv);
2311 if (output_env_get_print_xvgr_codes(oenv))
2313 fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
2320 get_dhdl_parms(ftp2fn(efTPX, NFILE, fnm), &ir);
2323 /* Initiate energies and set them to zero */
2332 /* Initiate counters */
2335 bFoundStart = FALSE;
2340 /* This loop searches for the first frame (when -b option is given),
2341 * or when this has been found it reads just one energy frame
2345 bCont = do_enx(fp, &(frame[NEXT]));
2348 timecheck = check_times(frame[NEXT].t);
2351 while (bCont && (timecheck < 0));
2353 if ((timecheck == 0) && bCont)
2355 /* We read a valid frame, so we can use it */
2356 fr = &(frame[NEXT]);
2360 /* The frame contains energies, so update cur */
2363 if (edat.nframes % 1000 == 0)
2365 srenew(edat.step, edat.nframes+1000);
2366 memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
2367 srenew(edat.steps, edat.nframes+1000);
2368 memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
2369 srenew(edat.points, edat.nframes+1000);
2370 memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
2372 for (i = 0; i < nset; i++)
2374 srenew(edat.s[i].ener, edat.nframes+1000);
2375 memset(&(edat.s[i].ener[edat.nframes]), 0,
2376 1000*sizeof(edat.s[i].ener[0]));
2377 srenew(edat.s[i].es, edat.nframes+1000);
2378 memset(&(edat.s[i].es[edat.nframes]), 0,
2379 1000*sizeof(edat.s[i].es[0]));
2384 edat.step[nfr] = fr->step;
2389 /* Initiate the previous step data */
2390 start_step = fr->step;
2392 /* Initiate the energy sums */
2393 edat.steps[nfr] = 1;
2394 edat.points[nfr] = 1;
2395 for (i = 0; i < nset; i++)
2398 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2399 edat.s[i].es[nfr].sum2 = 0;
2406 edat.steps[nfr] = fr->nsteps;
2408 if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2412 edat.points[nfr] = 1;
2413 for (i = 0; i < nset; i++)
2416 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2417 edat.s[i].es[nfr].sum2 = 0;
2423 edat.points[nfr] = fr->nsum;
2424 for (i = 0; i < nset; i++)
2427 edat.s[i].es[nfr].sum = fr->ener[sss].esum;
2428 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2430 edat.npoints += fr->nsum;
2435 /* The interval does not match fr->nsteps:
2436 * can not do exact averages.
2440 edat.nsteps = fr->step - start_step + 1;
2443 for (i = 0; i < nset; i++)
2445 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2449 * Define distance restraint legends. Can only be done after
2450 * the first frame has been read... (Then we know how many there are)
2452 blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL);
2453 if (bDisRe && bDRAll && !leg && blk_disre)
2458 fa = top->idef.il[F_DISRES].iatoms;
2459 ip = top->idef.iparams;
2460 if (blk_disre->nsub != 2 ||
2461 (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2463 gmx_incons("Number of disre sub-blocks not equal to 2");
2466 ndisre = blk_disre->sub[0].nr;
2467 if (ndisre != top->idef.il[F_DISRES].nr/3)
2469 gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2470 ndisre, top->idef.il[F_DISRES].nr/3);
2472 snew(pairleg, ndisre);
2473 for (i = 0; i < ndisre; i++)
2475 snew(pairleg[i], 30);
2478 gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j);
2479 gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k);
2480 sprintf(pairleg[i], "%d %s %d %s (%d)",
2481 resnr_j, anm_j, resnr_k, anm_k,
2482 ip[fa[3*i]].disres.label);
2484 set = select_it(ndisre, pairleg, &nset);
2486 for (i = 0; (i < nset); i++)
2489 sprintf(leg[2*i], "a %s", pairleg[set[i]]);
2490 snew(leg[2*i+1], 32);
2491 sprintf(leg[2*i+1], "i %s", pairleg[set[i]]);
2493 xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv);
2497 * Store energies for analysis afterwards...
2499 if (!bDisRe && !bDHDL && (fr->nre > 0))
2501 if (edat.nframes % 1000 == 0)
2503 srenew(time, edat.nframes+1000);
2505 time[edat.nframes] = fr->t;
2509 * Printing time, only when we do not want to skip frames
2511 if (!skip || teller % skip == 0)
2515 /*******************************************
2516 * D I S T A N C E R E S T R A I N T S
2517 *******************************************/
2521 float *disre_rt = blk_disre->sub[0].fval;
2522 float *disre_rm3tav = blk_disre->sub[1].fval;
2524 double *disre_rt = blk_disre->sub[0].dval;
2525 double *disre_rm3tav = blk_disre->sub[1].dval;
2528 print_time(out, fr->t);
2529 if (violaver == NULL)
2531 snew(violaver, ndisre);
2534 /* Subtract bounds from distances, to calculate violations */
2535 calc_violations(disre_rt, disre_rm3tav,
2536 nbounds, pair, bounds, violaver, &sumt, &sumaver);
2538 fprintf(out, " %8.4f %8.4f\n", sumaver, sumt);
2541 print_time(fp_pairs, fr->t);
2542 for (i = 0; (i < nset); i++)
2545 fprintf(fp_pairs, " %8.4f", mypow(disre_rm3tav[sss], minthird));
2546 fprintf(fp_pairs, " %8.4f", disre_rt[sss]);
2548 fprintf(fp_pairs, "\n");
2555 do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2558 /*******************************************
2560 *******************************************/
2567 /* We skip frames with single points (usually only the first frame),
2568 * since they would result in an average plot with outliers.
2572 print_time(out, fr->t);
2573 print1(out, bDp, fr->ener[set[0]].e);
2574 print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
2575 print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum));
2581 print_time(out, fr->t);
2585 for (i = 0; i < nset; i++)
2587 sum += fr->ener[set[i]].e;
2589 print1(out, bDp, sum/nmol-ezero);
2593 for (i = 0; (i < nset); i++)
2597 print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero);
2601 print1(out, bDp, fr->ener[set[i]].e);
2608 blk = find_block_id_enxframe(fr, enx_i, NULL);
2612 xdr_datatype dt = xdr_datatype_float;
2614 xdr_datatype dt = xdr_datatype_double;
2618 if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2620 gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2623 vals = blk->sub[0].fval;
2625 vals = blk->sub[0].dval;
2628 if (blk->sub[0].nr != (size_t)nor)
2630 gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2634 for (i = 0; i < nor; i++)
2636 orient[i] += vals[i];
2641 for (i = 0; i < nor; i++)
2643 odrms[i] += sqr(vals[i]-oobs[i]);
2648 fprintf(fort, " %10f", fr->t);
2649 for (i = 0; i < norsel; i++)
2651 fprintf(fort, " %g", vals[orsel[i]]);
2653 fprintf(fort, "\n");
2657 fprintf(fodt, " %10f", fr->t);
2658 for (i = 0; i < norsel; i++)
2660 fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]);
2662 fprintf(fodt, "\n");
2666 blk = find_block_id_enxframe(fr, enxORT, NULL);
2670 xdr_datatype dt = xdr_datatype_float;
2672 xdr_datatype dt = xdr_datatype_double;
2676 if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2678 gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2681 vals = blk->sub[0].fval;
2683 vals = blk->sub[0].dval;
2686 if (blk->sub[0].nr != (size_t)(nex*12))
2688 gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2689 blk->sub[0].nr/12, nex);
2691 fprintf(foten, " %10f", fr->t);
2692 for (i = 0; i < nex; i++)
2694 for (j = 0; j < (bOvec ? 12 : 3); j++)
2696 fprintf(foten, " %g", vals[i*12+j]);
2699 fprintf(foten, "\n");
2706 while (bCont && (timecheck == 0));
2708 fprintf(stderr, "\n");
2717 gmx_ffclose(fp_pairs);
2730 out = xvgropen(opt2fn("-ora", NFILE, fnm),
2731 "Average calculated orientations",
2732 "Restraint label", "", oenv);
2733 if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2735 fprintf(out, "%s", orinst_sub);
2737 for (i = 0; i < nor; i++)
2739 fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr);
2745 out = xvgropen(opt2fn("-oda", NFILE, fnm),
2746 "Average restraint deviation",
2747 "Restraint label", "", oenv);
2748 if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2750 fprintf(out, "%s", orinst_sub);
2752 for (i = 0; i < nor; i++)
2754 fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr-oobs[i]);
2760 out = xvgropen(opt2fn("-odr", NFILE, fnm),
2761 "RMS orientation restraint deviations",
2762 "Restraint label", "", oenv);
2763 if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2765 fprintf(out, "%s", orinst_sub);
2767 for (i = 0; i < nor; i++)
2769 fprintf(out, "%5d %g\n", or_label[i], sqrt(odrms[i]/norfr));
2780 analyse_disre(opt2fn("-viol", NFILE, fnm),
2781 teller_disre, violaver, bounds, index, pair, nbounds, oenv);
2787 gmx_ffclose(fp_dhdl);
2788 printf("\n\nWrote %d lambda values with %d samples as ",
2789 dh_lambdas, dh_samples);
2792 printf("%d dH histograms ", dh_hists);
2796 printf("%d dH data blocks ", dh_blocks);
2798 printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2803 gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2809 double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2810 analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
2811 bFee, bSum, opt2parg_bSet("-nmol", npargs, ppa),
2812 bVisco, opt2fn("-vis", NFILE, fnm),
2814 start_step, start_t, frame[cur].step, frame[cur].t,
2816 nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
2820 calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat,
2824 if (opt2bSet("-f2", NFILE, fnm))
2826 fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm),
2827 reftemp, nset, set, leg, &edat, time, oenv);
2831 const char *nxy = "-nxy";
2833 do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2834 do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2835 do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
2836 do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
2837 do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
2838 do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
2839 do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
2840 do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
2841 do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);