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"
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("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("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 nframes, real **sum,
573 real V, real T, int nsteps, double time[],
574 const output_env_t oenv)
577 double av[4], avold[4];
586 dt = (time[1]-time[0]);
589 for (i = 0; i <= nsets; i++)
593 fp0 = xvgropen(fni, "Shear viscosity integral",
594 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
595 fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
596 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
597 for (i = 1; i < nf4; i++)
599 fac = dt*nframes/nsteps;
600 for (m = 0; m <= nsets; m++)
604 for (j = 0; j < nframes-i; j++)
606 for (m = 0; m < nsets; m++)
608 di = sqr(fac*(sum[m][j+i]-sum[m][j]));
611 av[nsets] += di/nsets;
614 /* Convert to SI for the viscosity */
615 fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i);
616 fprintf(fp0, "%10g", time[i]-time[0]);
617 for (m = 0; (m <= nsets); m++)
620 fprintf(fp0, " %10g", av[m]);
623 fprintf(fp1, "%10g", 0.5*(time[i]+time[i-1])-time[0]);
624 for (m = 0; (m <= nsets); m++)
626 fprintf(fp1, " %10g", (av[m]-avold[m])/dt);
649 static void clear_ee_sum(ee_sum_t *ees)
657 static void add_ee_sum(ee_sum_t *ees, double sum, int np)
663 static void add_ee_av(ee_sum_t *ees)
667 av = ees->sum/ees->np;
674 static double calc_ee2(int nb, ee_sum_t *ees)
676 return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
679 static void set_ee_av(ener_ee_t *eee)
683 char buf[STEPSTRSIZE];
684 fprintf(debug, "Storing average for err.est.: %s steps\n",
685 gmx_step_str(eee->nst, buf));
687 add_ee_av(&eee->sum);
689 if (eee->b == 1 || eee->nst < eee->nst_min)
691 eee->nst_min = eee->nst;
696 static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
699 double sum, sum2, sump, see2;
700 gmx_int64_t steps, np, p, bound_nb;
704 double x, sx, sy, sxx, sxy;
707 /* Check if we have exact statistics over all points */
708 for (i = 0; i < nset; i++)
711 ed->bExactStat = FALSE;
712 if (edat->npoints > 0)
714 /* All energy file sum entries 0 signals no exact sums.
715 * But if all energy values are 0, we still have exact sums.
718 for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
720 if (ed->ener[i] != 0)
724 ed->bExactStat = (ed->es[f].sum != 0);
728 ed->bExactStat = TRUE;
734 for (i = 0; i < nset; i++)
745 for (nb = nbmin; nb <= nbmax; nb++)
748 clear_ee_sum(&eee[nb].sum);
752 for (f = 0; f < edat->nframes; f++)
758 /* Add the sum and the sum of variances to the totals. */
764 sum2 += dsqr(sum/np - (sum + es->sum)/(np + p))
770 /* Add a single value to the sum and sum of squares. */
776 /* sum has to be increased after sum2 */
780 /* For the linear regression use variance 1/p.
781 * Note that sump is the sum, not the average, so we don't need p*.
783 x = edat->step[f] - 0.5*(edat->steps[f] - 1);
789 for (nb = nbmin; nb <= nbmax; nb++)
791 /* Check if the current end step is closer to the desired
792 * block boundary than the next end step.
794 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
795 if (eee[nb].nst > 0 &&
796 bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
806 eee[nb].nst += edat->step[f] - edat->step[f-1];
810 add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
814 add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
816 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
817 if (edat->step[f]*nb >= bound_nb)
824 edat->s[i].av = sum/np;
827 edat->s[i].rmsd = sqrt(sum2/np);
831 edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av));
834 if (edat->nframes > 1)
836 edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
840 edat->s[i].slope = 0;
845 for (nb = nbmin; nb <= nbmax; nb++)
847 /* Check if we actually got nb blocks and if the smallest
848 * block is not shorter than 80% of the average.
852 char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
853 fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
855 gmx_step_str(eee[nb].nst_min, buf1),
856 gmx_step_str(edat->nsteps, buf2));
858 if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
860 see2 += calc_ee2(nb, &eee[nb].sum);
866 edat->s[i].ee = sqrt(see2/nee);
876 static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax)
887 snew(s->ener, esum->nframes);
888 snew(s->es, esum->nframes);
890 s->bExactStat = TRUE;
892 for (i = 0; i < nset; i++)
894 if (!edat->s[i].bExactStat)
896 s->bExactStat = FALSE;
898 s->slope += edat->s[i].slope;
901 for (f = 0; f < edat->nframes; f++)
904 for (i = 0; i < nset; i++)
906 sum += edat->s[i].ener[f];
910 for (i = 0; i < nset; i++)
912 sum += edat->s[i].es[f].sum;
918 calc_averages(1, esum, nbmin, nbmax);
923 static char *ee_pr(double ee, char *buf)
930 sprintf(buf, "%s", "--");
934 /* Round to two decimals by printing. */
935 sprintf(tmp, "%.1e", ee);
936 sscanf(tmp, "%lf", &rnd);
937 sprintf(buf, "%g", rnd);
943 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat)
945 /* Remove the drift by performing a fit to y = ax+b.
946 Uses 5 iterations. */
948 double delta, d, sd, sd2;
950 edat->npoints = edat->nframes;
951 edat->nsteps = edat->nframes;
953 for (k = 0; (k < 5); k++)
955 for (i = 0; (i < nset); i++)
957 delta = edat->s[i].slope*dt;
961 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
964 for (j = 0; (j < edat->nframes); j++)
966 edat->s[i].ener[j] -= j*delta;
967 edat->s[i].es[j].sum = 0;
968 edat->s[i].es[j].sum2 = 0;
971 calc_averages(nset, edat, nbmin, nbmax);
975 static void calc_fluctuation_props(FILE *fp,
976 gmx_bool bDriftCorr, real dt,
978 char **leg, enerdata_t *edat,
979 int nbmin, int nbmax)
982 double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
985 eVol, eEnth, eTemp, eEtot, eNR
987 const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
990 NANO3 = NANO*NANO*NANO;
993 fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n"
994 "for spurious drift in the graphs. Note that this is not\n"
995 "a substitute for proper equilibration and sampling!\n");
999 remove_drift(nset, nbmin, nbmax, dt, edat);
1001 for (i = 0; (i < eNR); i++)
1003 for (ii[i] = 0; (ii[i] < nset &&
1004 (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++)
1008 /* if (ii[i] < nset)
1009 fprintf(fp,"Found %s data.\n",my_ener[i]);
1011 /* Compute it all! */
1012 alpha = kappa = cp = dcp = cv = NOTSET;
1016 if (ii[eTemp] < nset)
1018 tt = edat->s[ii[eTemp]].av;
1022 if ((ii[eVol] < nset) && (ii[eTemp] < nset))
1024 vv = edat->s[ii[eVol]].av*NANO3;
1025 varv = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
1026 kappa = (varv/vv)/(BOLTZMANN*tt);
1030 if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
1032 hh = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
1033 varh = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
1034 cp = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
1037 et = varet = NOTSET;
1038 if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
1040 /* Only compute cv in constant volume runs, which we can test
1041 by checking whether the enthalpy was computed.
1043 et = edat->s[ii[eEtot]].av;
1044 varet = sqr(edat->s[ii[eEtot]].rmsd);
1045 cv = KILO*((varet/nmol)/(BOLTZ*tt*tt));
1048 if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
1050 double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
1051 vh_sum = v_sum = h_sum = 0;
1052 for (j = 0; (j < edat->nframes); j++)
1054 v = edat->s[ii[eVol]].ener[j]*NANO3;
1055 h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
1060 vh_aver = vh_sum / edat->nframes;
1061 v_aver = v_sum / edat->nframes;
1062 h_aver = h_sum / edat->nframes;
1063 alpha = (vh_aver-v_aver*h_aver)/(v_aver*BOLTZMANN*tt*tt);
1064 dcp = (v_aver*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
1071 fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n",
1074 fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
1075 fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
1076 fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
1077 fprintf(fp, "please use the g_dos program.\n\n");
1078 fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n"
1079 "a block-averaging error analysis (not implemented in g_energy yet)\n");
1085 fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol);
1090 fprintf(fp, "Volume = %10g m^3/mol\n",
1095 fprintf(fp, "Enthalpy = %10g kJ/mol\n",
1096 hh*AVOGADRO/(KILO*nmol));
1098 if (alpha != NOTSET)
1100 fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
1103 if (kappa != NOTSET)
1105 fprintf(fp, "Isothermal Compressibility Kappa = %10g (J/m^3)\n",
1107 fprintf(fp, "Adiabatic bulk modulus = %10g (m^3/J)\n",
1112 fprintf(fp, "Heat capacity at constant pressure Cp = %10g J/mol K\n",
1117 fprintf(fp, "Heat capacity at constant volume Cv = %10g J/mol K\n",
1122 fprintf(fp, "Cp-Cv = %10g J/mol K\n",
1125 please_cite(fp, "Allen1987a");
1129 fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
1133 static void analyse_ener(gmx_bool bCorr, const char *corrfn,
1134 gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct,
1135 gmx_bool bVisco, const char *visfn, int nmol,
1136 gmx_int64_t start_step, double start_t,
1137 gmx_int64_t step, double t,
1138 double time[], real reftemp,
1140 int nset, int set[], gmx_bool *bIsEner,
1141 char **leg, gmx_enxnm_t *enm,
1142 real Vaver, real ezero,
1143 int nbmin, int nbmax,
1144 const output_env_t oenv)
1147 /* Check out the printed manual for equations! */
1148 double Dt, aver, stddev, errest, delta_t, totaldrift;
1149 enerdata_t *esum = NULL;
1150 real xxx, integral, intBulk, Temp = 0, Pres = 0;
1151 real sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest;
1152 double beta = 0, expE, expEtot, *fee = NULL;
1154 int nexact, nnotexact;
1158 char buf[256], eebuf[100];
1160 nsteps = step - start_step + 1;
1163 fprintf(stdout, "Not enough steps (%s) for statistics\n",
1164 gmx_step_str(nsteps, buf));
1168 /* Calculate the time difference */
1169 delta_t = t - start_t;
1171 fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
1172 gmx_step_str(nsteps, buf), start_t, t, nset);
1174 calc_averages(nset, edat, nbmin, nbmax);
1178 esum = calc_sum(nset, edat, nbmin, nbmax);
1181 if (edat->npoints == 0)
1190 for (i = 0; (i < nset); i++)
1192 if (edat->s[i].bExactStat)
1205 fprintf(stdout, "All statistics are over %s points\n",
1206 gmx_step_str(edat->npoints, buf));
1208 else if (nexact == 0 || edat->npoints == edat->nframes)
1210 fprintf(stdout, "All statistics are over %d points (frames)\n",
1215 fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
1216 for (i = 0; (i < nset); i++)
1218 if (!edat->s[i].bExactStat)
1220 fprintf(stdout, " '%s'", leg[i]);
1223 fprintf(stdout, " %s has statistics over %d points (frames)\n",
1224 nnotexact == 1 ? "is" : "are", edat->nframes);
1225 fprintf(stdout, "All other statistics are over %s points\n",
1226 gmx_step_str(edat->npoints, buf));
1228 fprintf(stdout, "\n");
1230 fprintf(stdout, "%-24s %10s %10s %10s %10s",
1231 "Energy", "Average", "Err.Est.", "RMSD", "Tot-Drift");
1234 fprintf(stdout, " %10s\n", "-kT ln<e^(E/kT)>");
1238 fprintf(stdout, "\n");
1240 fprintf(stdout, "-------------------------------------------------------------------------------\n");
1242 /* Initiate locals, only used with -sum */
1246 beta = 1.0/(BOLTZ*reftemp);
1249 for (i = 0; (i < nset); i++)
1251 aver = edat->s[i].av;
1252 stddev = edat->s[i].rmsd;
1253 errest = edat->s[i].ee;
1258 for (j = 0; (j < edat->nframes); j++)
1260 expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
1264 expEtot += expE/edat->nframes;
1267 fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
1269 if (strstr(leg[i], "empera") != NULL)
1273 else if (strstr(leg[i], "olum") != NULL)
1277 else if (strstr(leg[i], "essure") != NULL)
1283 pr_aver = aver/nmol-ezero;
1284 pr_stddev = stddev/nmol;
1285 pr_errest = errest/nmol;
1294 /* Multiply the slope in steps with the number of steps taken */
1295 totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
1301 fprintf(stdout, "%-24s %10g %10s %10g %10g",
1302 leg[i], pr_aver, ee_pr(pr_errest, eebuf), pr_stddev, totaldrift);
1305 fprintf(stdout, " %10g", fee[i]);
1308 fprintf(stdout, " (%s)\n", enm[set[i]].unit);
1312 for (j = 0; (j < edat->nframes); j++)
1314 edat->s[i].ener[j] -= aver;
1320 totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1321 fprintf(stdout, "%-24s %10g %10s %10s %10g (%s)",
1322 "Total", esum->s[0].av/nmol, ee_pr(esum->s[0].ee/nmol, eebuf),
1323 "--", totaldrift/nmol, enm[set[0]].unit);
1324 /* pr_aver,pr_stddev,a,totaldrift */
1327 fprintf(stdout, " %10g %10g\n",
1328 log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta);
1332 fprintf(stdout, "\n");
1336 /* Do correlation function */
1337 if (edat->nframes > 1)
1339 Dt = delta_t/(edat->nframes - 1);
1347 const char* leg[] = { "Shear", "Bulk" };
1352 /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1354 /* Symmetrise tensor! (and store in first three elements)
1355 * And subtract average pressure!
1358 for (i = 0; i < 12; i++)
1360 snew(eneset[i], edat->nframes);
1363 for (i = 0; i < 3; i++)
1365 snew(enesum[i], edat->nframes);
1367 for (i = 0; (i < edat->nframes); i++)
1369 eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1370 eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1371 eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1372 for (j = 3; j <= 11; j++)
1374 eneset[j][i] = edat->s[j].ener[i];
1376 eneset[11][i] -= Pres;
1377 enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum);
1378 enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
1379 enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
1382 einstein_visco("evisco.xvg", "eviscoi.xvg",
1383 3, edat->nframes, enesum, Vaver, Temp, nsteps, time, oenv);
1385 /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1386 /* Do it for shear viscosity */
1387 strcpy(buf, "Shear Viscosity");
1388 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3,
1389 (edat->nframes+1)/2, eneset, Dt,
1390 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1392 /* Now for bulk viscosity */
1393 strcpy(buf, "Bulk Viscosity");
1394 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1,
1395 (edat->nframes+1)/2, &(eneset[11]), Dt,
1396 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1398 factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1399 fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1400 xvgr_legend(fp, asize(leg), leg, oenv);
1402 /* Use trapezium rule for integration */
1405 nout = get_acfnout();
1406 if ((nout < 2) || (nout >= edat->nframes/2))
1408 nout = edat->nframes/2;
1410 for (i = 1; (i < nout); i++)
1412 integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor;
1413 intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1414 fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk);
1422 strcpy(buf, "Autocorrelation of Energy Fluctuations");
1426 strcpy(buf, "Energy Autocorrelation");
1429 do_autocorr(corrfn, oenv, buf, edat->nframes,
1431 bSum ? &edat->s[nset-1].ener : eneset,
1432 (delta_t/edat->nframes), eacNormal, FALSE);
1438 static void print_time(FILE *fp, double t)
1440 fprintf(fp, "%12.6f", t);
1443 static void print1(FILE *fp, gmx_bool bDp, real e)
1447 fprintf(fp, " %16.12f", e);
1451 fprintf(fp, " %10.6f", e);
1455 static void fec(const char *ene2fn, const char *runavgfn,
1456 real reftemp, int nset, int set[], char *leg[],
1457 enerdata_t *edat, double time[],
1458 const output_env_t oenv)
1460 const char * ravgleg[] = {
1461 "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1462 "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N"
1466 int nre, timecheck, step, nenergy, nenergy2, maxenergy;
1472 gmx_enxnm_t *enm = NULL;
1476 /* read second energy file */
1479 enx = open_enx(ene2fn, "r");
1480 do_enxnms(enx, &(fr->nre), &enm);
1482 snew(eneset2, nset+1);
1488 /* This loop searches for the first frame (when -b option is given),
1489 * or when this has been found it reads just one energy frame
1493 bCont = do_enx(enx, fr);
1497 timecheck = check_times(fr->t);
1501 while (bCont && (timecheck < 0));
1503 /* Store energies for analysis afterwards... */
1504 if ((timecheck == 0) && bCont)
1508 if (nenergy2 >= maxenergy)
1511 for (i = 0; i <= nset; i++)
1513 srenew(eneset2[i], maxenergy);
1516 if (fr->t != time[nenergy2])
1518 fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n",
1519 fr->t, time[nenergy2], gmx_step_str(fr->step, buf));
1521 for (i = 0; i < nset; i++)
1523 eneset2[i][nenergy2] = fr->ener[set[i]].e;
1529 while (bCont && (timecheck == 0));
1532 if (edat->nframes != nenergy2)
1534 fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n",
1535 edat->nframes, nenergy2);
1537 nenergy = min(edat->nframes, nenergy2);
1539 /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1543 fp = xvgropen(runavgfn, "Running average free energy difference",
1544 "Time (" unit_time ")", "\\8D\\4E (" unit_energy ")", oenv);
1545 xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
1547 fprintf(stdout, "\n%-24s %10s\n",
1548 "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
1550 beta = 1.0/(BOLTZ*reftemp);
1551 for (i = 0; i < nset; i++)
1553 if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
1555 fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n",
1556 leg[i], enm[set[i]].name);
1558 for (j = 0; j < nenergy; j++)
1560 dE = eneset2[i][j] - edat->s[i].ener[j];
1561 sum += exp(-dE*beta);
1564 fprintf(fp, "%10g %10g %10g\n",
1565 time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
1568 aver = -BOLTZ*reftemp*log(sum/nenergy);
1569 fprintf(stdout, "%-24s %10g\n", leg[i], aver);
1579 static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl,
1580 const char *filename, gmx_bool bDp,
1581 int *blocks, int *hists, int *samples, int *nlambdas,
1582 const output_env_t oenv)
1584 const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1585 char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
1587 gmx_bool first = FALSE;
1588 int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
1591 double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0;
1592 static int setnr = 0;
1593 double *native_lambda_vec = NULL;
1594 const char **lambda_components = NULL;
1595 int n_lambda_vec = 0;
1596 gmx_bool changing_lambda = FALSE;
1597 int lambda_fep_state;
1599 /* now count the blocks & handle the global dh data */
1600 for (i = 0; i < fr->nblock; i++)
1602 if (fr->block[i].id == enxDHHIST)
1606 else if (fr->block[i].id == enxDH)
1610 else if (fr->block[i].id == enxDHCOLL)
1613 if ( (fr->block[i].nsub < 1) ||
1614 (fr->block[i].sub[0].type != xdr_datatype_double) ||
1615 (fr->block[i].sub[0].nr < 5))
1617 gmx_fatal(FARGS, "Unexpected block data");
1620 /* read the data from the DHCOLL block */
1621 temp = fr->block[i].sub[0].dval[0];
1622 start_time = fr->block[i].sub[0].dval[1];
1623 delta_time = fr->block[i].sub[0].dval[2];
1624 start_lambda = fr->block[i].sub[0].dval[3];
1625 delta_lambda = fr->block[i].sub[0].dval[4];
1626 changing_lambda = (delta_lambda != 0);
1627 if (fr->block[i].nsub > 1)
1629 lambda_fep_state = fr->block[i].sub[1].ival[0];
1630 if (n_lambda_vec == 0)
1632 n_lambda_vec = fr->block[i].sub[1].ival[1];
1636 if (n_lambda_vec != fr->block[i].sub[1].ival[1])
1639 "Unexpected change of basis set in lambda");
1642 if (lambda_components == NULL)
1644 snew(lambda_components, n_lambda_vec);
1646 if (native_lambda_vec == NULL)
1648 snew(native_lambda_vec, n_lambda_vec);
1650 for (j = 0; j < n_lambda_vec; j++)
1652 native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j];
1653 lambda_components[j] =
1654 efpt_singular_names[fr->block[i].sub[1].ival[2+j]];
1660 if (nblock_hist == 0 && nblock_dh == 0)
1662 /* don't do anything */
1665 if (nblock_hist > 0 && nblock_dh > 0)
1667 gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1673 /* we have standard, non-histogram data --
1674 call open_dhdl to open the file */
1675 /* TODO this is an ugly hack that needs to be fixed: this will only
1676 work if the order of data is always the same and if we're
1677 only using the g_energy compiled with the mdrun that produced
1679 *fp_dhdl = open_dhdl(filename, ir, oenv);
1683 sprintf(title, "N(%s)", deltag);
1684 sprintf(label_x, "%s (%s)", deltag, unit_energy);
1685 sprintf(label_y, "Samples");
1686 *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1687 sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
1688 xvgr_subtitle(*fp_dhdl, buf, oenv);
1692 (*hists) += nblock_hist;
1693 (*blocks) += nblock_dh;
1694 (*nlambdas) = nblock_hist+nblock_dh;
1696 /* write the data */
1697 if (nblock_hist > 0)
1699 gmx_int64_t sum = 0;
1701 for (i = 0; i < fr->nblock; i++)
1703 t_enxblock *blk = &(fr->block[i]);
1704 if (blk->id == enxDHHIST)
1706 double foreign_lambda, dx;
1708 int nhist, derivative;
1710 /* check the block types etc. */
1711 if ( (blk->nsub < 2) ||
1712 (blk->sub[0].type != xdr_datatype_double) ||
1713 (blk->sub[1].type != xdr_datatype_int64) ||
1714 (blk->sub[0].nr < 2) ||
1715 (blk->sub[1].nr < 2) )
1717 gmx_fatal(FARGS, "Unexpected block data in file");
1719 foreign_lambda = blk->sub[0].dval[0];
1720 dx = blk->sub[0].dval[1];
1721 nhist = blk->sub[1].lval[0];
1722 derivative = blk->sub[1].lval[1];
1723 for (j = 0; j < nhist; j++)
1726 x0 = blk->sub[1].lval[2+j];
1730 sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1731 deltag, lambda, foreign_lambda,
1732 lambda, start_lambda);
1736 sprintf(legend, "N(%s | %s=%g)",
1737 dhdl, lambda, start_lambda);
1741 xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1743 for (k = 0; k < blk->sub[j+2].nr; k++)
1748 hist = blk->sub[j+2].ival[k];
1751 fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist,
1755 /* multiple histogram data blocks in one histogram
1756 mean that the second one is the reverse of the first one:
1757 for dhdl derivatives, it's important to know both the
1758 maximum and minimum values */
1763 (*samples) += (int)(sum/nblock_hist);
1769 char **setnames = NULL;
1770 int nnames = nblock_dh;
1772 for (i = 0; i < fr->nblock; i++)
1774 t_enxblock *blk = &(fr->block[i]);
1775 if (blk->id == enxDH)
1779 len = blk->sub[2].nr;
1783 if (len != blk->sub[2].nr)
1785 gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1792 for (i = 0; i < len; i++)
1794 double time = start_time + delta_time*i;
1796 fprintf(*fp_dhdl, "%.4f ", time);
1798 for (j = 0; j < fr->nblock; j++)
1800 t_enxblock *blk = &(fr->block[j]);
1801 if (blk->id == enxDH)
1804 if (blk->sub[2].type == xdr_datatype_float)
1806 value = blk->sub[2].fval[i];
1810 value = blk->sub[2].dval[i];
1812 /* we need to decide which data type it is based on the count*/
1814 if (j == 1 && ir->bExpanded)
1816 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! */
1822 fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */
1826 fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */
1831 fprintf(*fp_dhdl, "\n");
1837 int gmx_energy(int argc, char *argv[])
1839 const char *desc[] = {
1840 "[THISMODULE] extracts energy components or distance restraint",
1841 "data from an energy file. The user is prompted to interactively",
1842 "select the desired energy terms.[PAR]",
1844 "Average, RMSD, and drift are calculated with full precision from the",
1845 "simulation (see printed manual). Drift is calculated by performing",
1846 "a least-squares fit of the data to a straight line. The reported total drift",
1847 "is the difference of the fit at the first and last point.",
1848 "An error estimate of the average is given based on a block averages",
1849 "over 5 blocks using the full-precision averages. The error estimate",
1850 "can be performed over multiple block lengths with the options",
1851 "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1852 "[BB]Note[bb] that in most cases the energy files contains averages over all",
1853 "MD steps, or over many more points than the number of frames in",
1854 "energy file. This makes the [THISMODULE] statistics output more accurate",
1855 "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1856 "file, the statistics mentioned above are simply over the single, per-frame",
1857 "energy values.[PAR]",
1859 "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1861 "Some fluctuation-dependent properties can be calculated provided",
1862 "the correct energy terms are selected, and that the command line option",
1863 "[TT]-fluct_props[tt] is given. The following properties",
1864 "will be computed:[BR]",
1865 "Property Energy terms needed[BR]",
1866 "---------------------------------------------------[BR]",
1867 "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp [BR]",
1868 "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp [BR]",
1869 "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1870 "Isothermal compressibility: Vol, Temp [BR]",
1871 "Adiabatic bulk modulus: Vol, Temp [BR]",
1872 "---------------------------------------------------[BR]",
1873 "You always need to set the number of molecules [TT]-nmol[tt].",
1874 "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1875 "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]"
1876 "When the [TT]-viol[tt] option is set, the time averaged",
1877 "violations are plotted and the running time-averaged and",
1878 "instantaneous sum of violations are recalculated. Additionally",
1879 "running time-averaged and instantaneous distances between",
1880 "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1882 "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1883 "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1884 "The first two options plot the orientation, the last three the",
1885 "deviations of the orientations from the experimental values.",
1886 "The options that end on an 'a' plot the average over time",
1887 "as a function of restraint. The options that end on a 't'",
1888 "prompt the user for restraint label numbers and plot the data",
1889 "as a function of time. Option [TT]-odr[tt] plots the RMS",
1890 "deviation as a function of restraint.",
1891 "When the run used time or ensemble averaged orientation restraints,",
1892 "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1893 "not ensemble-averaged orientations and deviations instead of",
1894 "the time and ensemble averages.[PAR]",
1896 "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1897 "tensor for each orientation restraint experiment. With option",
1898 "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1900 "Option [TT]-odh[tt] extracts and plots the free energy data",
1901 "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1902 "from the [TT]ener.edr[tt] file.[PAR]",
1904 "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1905 "difference with an ideal gas state: [BR]",
1906 " [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]",
1907 " [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]",
1908 "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1909 "the average is over the ensemble (or time in a trajectory).",
1910 "Note that this is in principle",
1911 "only correct when averaging over the whole (Boltzmann) ensemble",
1912 "and using the potential energy. This also allows for an entropy",
1913 "estimate using:[BR]",
1914 " [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]",
1915 " [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",
1918 "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1919 "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1920 "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1921 "files, and the average is over the ensemble A. The running average",
1922 "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1923 "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1926 static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1927 static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1928 static int skip = 0, nmol = 1, nbmin = 5, nbmax = 5;
1929 static real reftemp = 300.0, ezero = 0;
1931 { "-fee", FALSE, etBOOL, {&bFee},
1932 "Do a free energy estimate" },
1933 { "-fetemp", FALSE, etREAL, {&reftemp},
1934 "Reference temperature for free energy calculation" },
1935 { "-zero", FALSE, etREAL, {&ezero},
1936 "Subtract a zero-point energy" },
1937 { "-sum", FALSE, etBOOL, {&bSum},
1938 "Sum the energy terms selected rather than display them all" },
1939 { "-dp", FALSE, etBOOL, {&bDp},
1940 "Print energies in high precision" },
1941 { "-nbmin", FALSE, etINT, {&nbmin},
1942 "Minimum number of blocks for error estimate" },
1943 { "-nbmax", FALSE, etINT, {&nbmax},
1944 "Maximum number of blocks for error estimate" },
1945 { "-mutot", FALSE, etBOOL, {&bMutot},
1946 "Compute the total dipole moment from the components" },
1947 { "-skip", FALSE, etINT, {&skip},
1948 "Skip number of frames between data points" },
1949 { "-aver", FALSE, etBOOL, {&bPrAll},
1950 "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1951 { "-nmol", FALSE, etINT, {&nmol},
1952 "Number of molecules in your sample: the energies are divided by this number" },
1953 { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1954 "Compute properties based on energy fluctuations, like heat capacity" },
1955 { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1956 "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1957 { "-fluc", FALSE, etBOOL, {&bFluct},
1958 "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1959 { "-orinst", FALSE, etBOOL, {&bOrinst},
1960 "Analyse instantaneous orientation data" },
1961 { "-ovec", FALSE, etBOOL, {&bOvec},
1962 "Also plot the eigenvectors with [TT]-oten[tt]" }
1964 const char * drleg[] = {
1968 static const char *setnm[] = {
1969 "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1970 "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1971 "Volume", "Pressure"
1974 FILE *out = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL;
1975 FILE *fp_dhdl = NULL;
1980 gmx_localtop_t *top = NULL;
1984 gmx_enxnm_t *enm = NULL;
1985 t_enxframe *frame, *fr = NULL;
1987 #define NEXT (1-cur)
1988 int nre, teller, teller_disre, nfr;
1989 gmx_int64_t start_step;
1990 int nor = 0, nex = 0, norfr = 0, enx_i = 0;
1992 real *bounds = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL;
1993 int *index = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL;
1994 int nbounds = 0, npairs;
1995 gmx_bool bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL;
1996 gmx_bool bFoundStart, bCont, bEDR, bVisco;
1997 double sum, sumaver, sumt, ener, dbl;
1998 double *time = NULL;
2000 int *set = NULL, i, j, k, nset, sss;
2001 gmx_bool *bIsEner = NULL;
2002 char **pairleg, **odtleg, **otenleg;
2005 char *anm_j, *anm_k, *resnm_j, *resnm_k;
2006 int resnr_j, resnr_k;
2007 const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
2010 t_enxblock *blk = NULL;
2011 t_enxblock *blk_disre = NULL;
2013 int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
2016 { efEDR, "-f", NULL, ffREAD },
2017 { efEDR, "-f2", NULL, ffOPTRD },
2018 { efTPX, "-s", NULL, ffOPTRD },
2019 { efXVG, "-o", "energy", ffWRITE },
2020 { efXVG, "-viol", "violaver", ffOPTWR },
2021 { efXVG, "-pairs", "pairs", ffOPTWR },
2022 { efXVG, "-ora", "orienta", ffOPTWR },
2023 { efXVG, "-ort", "orientt", ffOPTWR },
2024 { efXVG, "-oda", "orideva", ffOPTWR },
2025 { efXVG, "-odr", "oridevr", ffOPTWR },
2026 { efXVG, "-odt", "oridevt", ffOPTWR },
2027 { efXVG, "-oten", "oriten", ffOPTWR },
2028 { efXVG, "-corr", "enecorr", ffOPTWR },
2029 { efXVG, "-vis", "visco", ffOPTWR },
2030 { efXVG, "-ravg", "runavgdf", ffOPTWR },
2031 { efXVG, "-odh", "dhdl", ffOPTWR }
2033 #define NFILE asize(fnm)
2038 ppa = add_acf_pargs(&npargs, pa);
2039 if (!parse_common_args(&argc, argv,
2040 PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
2041 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
2046 bDRAll = opt2bSet("-pairs", NFILE, fnm);
2047 bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
2048 bORA = opt2bSet("-ora", NFILE, fnm);
2049 bORT = opt2bSet("-ort", NFILE, fnm);
2050 bODA = opt2bSet("-oda", NFILE, fnm);
2051 bODR = opt2bSet("-odr", NFILE, fnm);
2052 bODT = opt2bSet("-odt", NFILE, fnm);
2053 bORIRE = bORA || bORT || bODA || bODR || bODT;
2054 bOTEN = opt2bSet("-oten", NFILE, fnm);
2055 bDHDL = opt2bSet("-odh", NFILE, fnm);
2060 fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
2061 do_enxnms(fp, &nre, &enm);
2065 bVisco = opt2bSet("-vis", NFILE, fnm);
2067 if ((!bDisRe) && (!bDHDL))
2071 nset = asize(setnm);
2073 /* This is nasty code... To extract Pres tensor, Volume and Temperature */
2074 for (j = 0; j < nset; j++)
2076 for (i = 0; i < nre; i++)
2078 if (strstr(enm[i].name, setnm[j]))
2086 if (gmx_strcasecmp(setnm[j], "Volume") == 0)
2088 printf("Enter the box volume (" unit_volume "): ");
2089 if (1 != scanf("%lf", &dbl))
2091 gmx_fatal(FARGS, "Error reading user input");
2097 gmx_fatal(FARGS, "Could not find term %s for viscosity calculation",
2105 set = select_by_name(nre, enm, &nset);
2107 /* Print all the different units once */
2108 sprintf(buf, "(%s)", enm[set[0]].unit);
2109 for (i = 1; i < nset; i++)
2111 for (j = 0; j < i; j++)
2113 if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
2121 strcat(buf, enm[set[i]].unit);
2125 out = xvgropen(opt2fn("-o", NFILE, fnm), "Gromacs Energies", "Time (ps)", buf,
2129 for (i = 0; (i < nset); i++)
2131 leg[i] = enm[set[i]].name;
2135 leg[nset] = strdup("Sum");
2136 xvgr_legend(out, nset+1, (const char**)leg, oenv);
2140 xvgr_legend(out, nset, (const char**)leg, oenv);
2143 snew(bIsEner, nset);
2144 for (i = 0; (i < nset); i++)
2147 for (j = 0; (j <= F_ETOT); j++)
2149 bIsEner[i] = bIsEner[i] ||
2150 (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
2154 if (bPrAll && nset > 1)
2156 gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
2161 if (bORIRE || bOTEN)
2163 get_orires_parms(ftp2fn(efTPX, NFILE, fnm), &nor, &nex, &or_label, &oobs);
2187 fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n");
2188 fprintf(stderr, "End your selection with 0\n");
2195 if (1 != scanf("%d", &(orsel[j])))
2197 gmx_fatal(FARGS, "Error reading user input");
2200 while (orsel[j] > 0);
2203 fprintf(stderr, "Selecting all %d orientation restraints\n", nor);
2206 for (i = 0; i < nor; i++)
2213 /* Build the selection */
2215 for (i = 0; i < j; i++)
2217 for (k = 0; k < nor; k++)
2219 if (or_label[k] == orsel[i])
2228 fprintf(stderr, "Orientation restraint label %d not found\n",
2233 snew(odtleg, norsel);
2234 for (i = 0; i < norsel; i++)
2236 snew(odtleg[i], 256);
2237 sprintf(odtleg[i], "%d", or_label[orsel[i]]);
2241 fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
2242 "Time (ps)", "", oenv);
2245 fprintf(fort, "%s", orinst_sub);
2247 xvgr_legend(fort, norsel, (const char**)odtleg, oenv);
2251 fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
2252 "Orientation restraint deviation",
2253 "Time (ps)", "", oenv);
2256 fprintf(fodt, "%s", orinst_sub);
2258 xvgr_legend(fodt, norsel, (const char**)odtleg, oenv);
2264 foten = xvgropen(opt2fn("-oten", NFILE, fnm),
2265 "Order tensor", "Time (ps)", "", oenv);
2266 snew(otenleg, bOvec ? nex*12 : nex*3);
2267 for (i = 0; i < nex; i++)
2269 for (j = 0; j < 3; j++)
2271 sprintf(buf, "eig%d", j+1);
2272 otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
2276 for (j = 0; j < 9; j++)
2278 sprintf(buf, "vec%d%s", j/3+1, j%3 == 0 ? "x" : (j%3 == 1 ? "y" : "z"));
2279 otenleg[12*i+3+j] = strdup(buf);
2283 xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv);
2288 nbounds = get_bounds(ftp2fn(efTPX, NFILE, fnm), &bounds, &index, &pair, &npairs,
2290 snew(violaver, npairs);
2291 out = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations",
2292 "Time (ps)", "nm", oenv);
2293 xvgr_legend(out, 2, drleg, oenv);
2296 fp_pairs = xvgropen(opt2fn("-pairs", NFILE, fnm), "Pair Distances",
2297 "Time (ps)", "Distance (nm)", oenv);
2298 if (output_env_get_print_xvgr_codes(oenv))
2300 fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
2307 get_dhdl_parms(ftp2fn(efTPX, NFILE, fnm), &ir);
2310 /* Initiate energies and set them to zero */
2319 /* Initiate counters */
2322 bFoundStart = FALSE;
2327 /* This loop searches for the first frame (when -b option is given),
2328 * or when this has been found it reads just one energy frame
2332 bCont = do_enx(fp, &(frame[NEXT]));
2335 timecheck = check_times(frame[NEXT].t);
2338 while (bCont && (timecheck < 0));
2340 if ((timecheck == 0) && bCont)
2342 /* We read a valid frame, so we can use it */
2343 fr = &(frame[NEXT]);
2347 /* The frame contains energies, so update cur */
2350 if (edat.nframes % 1000 == 0)
2352 srenew(edat.step, edat.nframes+1000);
2353 memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
2354 srenew(edat.steps, edat.nframes+1000);
2355 memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
2356 srenew(edat.points, edat.nframes+1000);
2357 memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
2359 for (i = 0; i < nset; i++)
2361 srenew(edat.s[i].ener, edat.nframes+1000);
2362 memset(&(edat.s[i].ener[edat.nframes]), 0,
2363 1000*sizeof(edat.s[i].ener[0]));
2364 srenew(edat.s[i].es, edat.nframes+1000);
2365 memset(&(edat.s[i].es[edat.nframes]), 0,
2366 1000*sizeof(edat.s[i].es[0]));
2371 edat.step[nfr] = fr->step;
2376 /* Initiate the previous step data */
2377 start_step = fr->step;
2379 /* Initiate the energy sums */
2380 edat.steps[nfr] = 1;
2381 edat.points[nfr] = 1;
2382 for (i = 0; i < nset; i++)
2385 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2386 edat.s[i].es[nfr].sum2 = 0;
2393 edat.steps[nfr] = fr->nsteps;
2395 if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2399 edat.points[nfr] = 1;
2400 for (i = 0; i < nset; i++)
2403 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2404 edat.s[i].es[nfr].sum2 = 0;
2410 edat.points[nfr] = fr->nsum;
2411 for (i = 0; i < nset; i++)
2414 edat.s[i].es[nfr].sum = fr->ener[sss].esum;
2415 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2417 edat.npoints += fr->nsum;
2422 /* The interval does not match fr->nsteps:
2423 * can not do exact averages.
2427 edat.nsteps = fr->step - start_step + 1;
2430 for (i = 0; i < nset; i++)
2432 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2436 * Define distance restraint legends. Can only be done after
2437 * the first frame has been read... (Then we know how many there are)
2439 blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL);
2440 if (bDisRe && bDRAll && !leg && blk_disre)
2445 fa = top->idef.il[F_DISRES].iatoms;
2446 ip = top->idef.iparams;
2447 if (blk_disre->nsub != 2 ||
2448 (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2450 gmx_incons("Number of disre sub-blocks not equal to 2");
2453 ndisre = blk_disre->sub[0].nr;
2454 if (ndisre != top->idef.il[F_DISRES].nr/3)
2456 gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2457 ndisre, top->idef.il[F_DISRES].nr/3);
2459 snew(pairleg, ndisre);
2460 for (i = 0; i < ndisre; i++)
2462 snew(pairleg[i], 30);
2465 gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j);
2466 gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k);
2467 sprintf(pairleg[i], "%d %s %d %s (%d)",
2468 resnr_j, anm_j, resnr_k, anm_k,
2469 ip[fa[3*i]].disres.label);
2471 set = select_it(ndisre, pairleg, &nset);
2473 for (i = 0; (i < nset); i++)
2476 sprintf(leg[2*i], "a %s", pairleg[set[i]]);
2477 snew(leg[2*i+1], 32);
2478 sprintf(leg[2*i+1], "i %s", pairleg[set[i]]);
2480 xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv);
2484 * Store energies for analysis afterwards...
2486 if (!bDisRe && !bDHDL && (fr->nre > 0))
2488 if (edat.nframes % 1000 == 0)
2490 srenew(time, edat.nframes+1000);
2492 time[edat.nframes] = fr->t;
2496 * Printing time, only when we do not want to skip frames
2498 if (!skip || teller % skip == 0)
2502 /*******************************************
2503 * D I S T A N C E R E S T R A I N T S
2504 *******************************************/
2508 float *disre_rt = blk_disre->sub[0].fval;
2509 float *disre_rm3tav = blk_disre->sub[1].fval;
2511 double *disre_rt = blk_disre->sub[0].dval;
2512 double *disre_rm3tav = blk_disre->sub[1].dval;
2515 print_time(out, fr->t);
2516 if (violaver == NULL)
2518 snew(violaver, ndisre);
2521 /* Subtract bounds from distances, to calculate violations */
2522 calc_violations(disre_rt, disre_rm3tav,
2523 nbounds, pair, bounds, violaver, &sumt, &sumaver);
2525 fprintf(out, " %8.4f %8.4f\n", sumaver, sumt);
2528 print_time(fp_pairs, fr->t);
2529 for (i = 0; (i < nset); i++)
2532 fprintf(fp_pairs, " %8.4f", mypow(disre_rm3tav[sss], minthird));
2533 fprintf(fp_pairs, " %8.4f", disre_rt[sss]);
2535 fprintf(fp_pairs, "\n");
2542 do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2545 /*******************************************
2547 *******************************************/
2554 /* We skip frames with single points (usually only the first frame),
2555 * since they would result in an average plot with outliers.
2559 print_time(out, fr->t);
2560 print1(out, bDp, fr->ener[set[0]].e);
2561 print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
2562 print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum));
2568 print_time(out, fr->t);
2572 for (i = 0; i < nset; i++)
2574 sum += fr->ener[set[i]].e;
2576 print1(out, bDp, sum/nmol-ezero);
2580 for (i = 0; (i < nset); i++)
2584 print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero);
2588 print1(out, bDp, fr->ener[set[i]].e);
2595 blk = find_block_id_enxframe(fr, enx_i, NULL);
2599 xdr_datatype dt = xdr_datatype_float;
2601 xdr_datatype dt = xdr_datatype_double;
2605 if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2607 gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2610 vals = blk->sub[0].fval;
2612 vals = blk->sub[0].dval;
2615 if (blk->sub[0].nr != (size_t)nor)
2617 gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2621 for (i = 0; i < nor; i++)
2623 orient[i] += vals[i];
2628 for (i = 0; i < nor; i++)
2630 odrms[i] += sqr(vals[i]-oobs[i]);
2635 fprintf(fort, " %10f", fr->t);
2636 for (i = 0; i < norsel; i++)
2638 fprintf(fort, " %g", vals[orsel[i]]);
2640 fprintf(fort, "\n");
2644 fprintf(fodt, " %10f", fr->t);
2645 for (i = 0; i < norsel; i++)
2647 fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]);
2649 fprintf(fodt, "\n");
2653 blk = find_block_id_enxframe(fr, enxORT, NULL);
2657 xdr_datatype dt = xdr_datatype_float;
2659 xdr_datatype dt = xdr_datatype_double;
2663 if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2665 gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2668 vals = blk->sub[0].fval;
2670 vals = blk->sub[0].dval;
2673 if (blk->sub[0].nr != (size_t)(nex*12))
2675 gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2676 blk->sub[0].nr/12, nex);
2678 fprintf(foten, " %10f", fr->t);
2679 for (i = 0; i < nex; i++)
2681 for (j = 0; j < (bOvec ? 12 : 3); j++)
2683 fprintf(foten, " %g", vals[i*12+j]);
2686 fprintf(foten, "\n");
2693 while (bCont && (timecheck == 0));
2695 fprintf(stderr, "\n");
2704 gmx_ffclose(fp_pairs);
2717 out = xvgropen(opt2fn("-ora", NFILE, fnm),
2718 "Average calculated orientations",
2719 "Restraint label", "", oenv);
2722 fprintf(out, "%s", orinst_sub);
2724 for (i = 0; i < nor; i++)
2726 fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr);
2732 out = xvgropen(opt2fn("-oda", NFILE, fnm),
2733 "Average restraint deviation",
2734 "Restraint label", "", oenv);
2737 fprintf(out, "%s", orinst_sub);
2739 for (i = 0; i < nor; i++)
2741 fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr-oobs[i]);
2747 out = xvgropen(opt2fn("-odr", NFILE, fnm),
2748 "RMS orientation restraint deviations",
2749 "Restraint label", "", oenv);
2752 fprintf(out, "%s", orinst_sub);
2754 for (i = 0; i < nor; i++)
2756 fprintf(out, "%5d %g\n", or_label[i], sqrt(odrms[i]/norfr));
2767 analyse_disre(opt2fn("-viol", NFILE, fnm),
2768 teller_disre, violaver, bounds, index, pair, nbounds, oenv);
2774 gmx_ffclose(fp_dhdl);
2775 printf("\n\nWrote %d lambda values with %d samples as ",
2776 dh_lambdas, dh_samples);
2779 printf("%d dH histograms ", dh_hists);
2783 printf("%d dH data blocks ", dh_blocks);
2785 printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2790 gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2796 double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2797 analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
2798 bFee, bSum, opt2parg_bSet("-nmol", npargs, ppa),
2799 bVisco, opt2fn("-vis", NFILE, fnm),
2801 start_step, start_t, frame[cur].step, frame[cur].t,
2802 time, reftemp, &edat,
2803 nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
2807 calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat,
2811 if (opt2bSet("-f2", NFILE, fnm))
2813 fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm),
2814 reftemp, nset, set, leg, &edat, time, oenv);
2818 const char *nxy = "-nxy";
2820 do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2821 do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2822 do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
2823 do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
2824 do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
2825 do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
2826 do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
2827 do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
2828 do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);