2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
46 #include "gmx_fatal.h"
50 #include "gromacs/fileio/enxio.h"
51 #include "gromacs/commandline/pargs.h"
58 #include "gromacs/fileio/tpxio.h"
59 #include "gromacs/fileio/trxio.h"
61 #include "mtop_util.h"
65 static real minthird = -1.0/3.0, minsixth = -1.0/6.0;
92 static double mypow(double x, double y)
104 static int *select_it(int nre, char *nm[], int *nset)
109 gmx_bool bVerbose = TRUE;
111 if ((getenv("VERBOSE")) != NULL)
116 fprintf(stderr, "Select the terms you want from the following list\n");
117 fprintf(stderr, "End your selection with 0\n");
121 for (k = 0; (k < nre); )
123 for (j = 0; (j < 4) && (k < nre); j++, k++)
125 fprintf(stderr, " %3d=%14s", k+1, nm[k]);
127 fprintf(stderr, "\n");
134 if (1 != scanf("%d", &n))
136 gmx_fatal(FARGS, "Error reading user input");
138 if ((n > 0) && (n <= nre))
146 for (i = (*nset) = 0; (i < nre); i++)
159 static void chomp(char *buf)
161 int len = strlen(buf);
163 while ((len > 0) && (buf[len-1] == '\n'))
170 static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
173 int n, k, kk, j, i, nmatch, nind, nss;
175 gmx_bool bEOF, bVerbose = TRUE, bLong = FALSE;
176 char *ptr, buf[STRLEN];
177 const char *fm4 = "%3d %-14s";
178 const char *fm2 = "%3d %-34s";
181 if ((getenv("VERBOSE")) != NULL)
186 fprintf(stderr, "\n");
187 fprintf(stderr, "Select the terms you want from the following list by\n");
188 fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
189 fprintf(stderr, "End your selection with an empty line or a zero.\n");
190 fprintf(stderr, "-------------------------------------------------------------------\n");
194 for (k = 0; k < nre; k++)
196 newnm[k] = strdup(nm[k].name);
197 /* Insert dashes in all the names */
198 while ((ptr = strchr(newnm[k], ' ')) != NULL)
208 fprintf(stderr, "\n");
211 for (kk = k; kk < k+4; kk++)
213 if (kk < nre && strlen(nm[kk].name) > 14)
221 fprintf(stderr, " ");
225 fprintf(stderr, fm4, k+1, newnm[k]);
234 fprintf(stderr, fm2, k+1, newnm[k]);
245 fprintf(stderr, "\n\n");
251 while (!bEOF && (fgets2(buf, STRLEN-1, stdin)))
253 /* Remove newlines */
259 /* Empty line means end of input */
260 bEOF = (strlen(buf) == 0);
268 /* First try to read an integer */
269 nss = sscanf(ptr, "%d", &nind);
272 /* Zero means end of input */
277 else if ((1 <= nind) && (nind <= nre))
283 fprintf(stderr, "number %d is out of range\n", nind);
288 /* Now try to read a string */
291 for (nind = 0; nind < nre; nind++)
293 if (gmx_strcasecmp(newnm[nind], ptr) == 0)
303 for (nind = 0; nind < nre; nind++)
305 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
313 fprintf(stderr, "String '%s' does not match anything\n", ptr);
318 /* Look for the first space, and remove spaces from there */
319 if ((ptr = strchr(ptr, ' ')) != NULL)
324 while (!bEOF && (ptr && (strlen(ptr) > 0)));
329 for (i = (*nset) = 0; (i < nre); i++)
341 gmx_fatal(FARGS, "No energy terms selected");
344 for (i = 0; (i < nre); i++)
353 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
360 /* all we need is the ir to be able to write the label */
361 read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, &mtop);
364 static void get_orires_parms(const char *topnm,
365 int *nor, int *nex, int **label, real **obs)
376 read_tpx(topnm, &ir, box, &natoms, NULL, NULL, NULL, &mtop);
377 top = gmx_mtop_generate_local_top(&mtop, &ir);
379 ip = top->idef.iparams;
380 iatom = top->idef.il[F_ORIRES].iatoms;
382 /* Count how many distance restraint there are... */
383 nb = top->idef.il[F_ORIRES].nr;
386 gmx_fatal(FARGS, "No orientation restraints in topology!\n");
393 for (i = 0; i < nb; i += 3)
395 (*label)[i/3] = ip[iatom[i]].orires.label;
396 (*obs)[i/3] = ip[iatom[i]].orires.obs;
397 if (ip[iatom[i]].orires.ex >= *nex)
399 *nex = ip[iatom[i]].orires.ex+1;
402 fprintf(stderr, "Found %d orientation restraints with %d experiments",
406 static int get_bounds(const char *topnm,
407 real **bounds, int **index, int **dr_pair, int *npairs,
408 gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir)
411 t_functype *functype;
413 int natoms, i, j, k, type, ftype, natom;
421 read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, mtop);
423 top = gmx_mtop_generate_local_top(mtop, ir);
426 functype = top->idef.functype;
427 ip = top->idef.iparams;
429 /* Count how many distance restraint there are... */
430 nb = top->idef.il[F_DISRES].nr;
433 gmx_fatal(FARGS, "No distance restraints in topology!\n");
436 /* Allocate memory */
441 /* Fill the bound array */
443 for (i = 0; (i < top->idef.ntypes); i++)
446 if (ftype == F_DISRES)
449 label1 = ip[i].disres.label;
450 b[nb] = ip[i].disres.up1;
457 /* Fill the index array */
459 disres = &(top->idef.il[F_DISRES]);
460 iatom = disres->iatoms;
461 for (i = j = k = 0; (i < disres->nr); )
464 ftype = top->idef.functype[type];
465 natom = interaction_function[ftype].nratoms+1;
466 if (label1 != top->idef.iparams[type].disres.label)
469 label1 = top->idef.iparams[type].disres.label;
479 gmx_incons("get_bounds for distance restraints");
488 static void calc_violations(real rt[], real rav3[], int nb, int index[],
489 real bounds[], real *viol, double *st, double *sa)
491 const real sixth = 1.0/6.0;
493 double rsum, rav, sumaver, sumt;
497 for (i = 0; (i < nb); i++)
501 for (j = index[i]; (j < index[i+1]); j++)
505 viol[j] += mypow(rt[j], -3.0);
508 rsum += mypow(rt[j], -6);
510 rsum = max(0.0, mypow(rsum, -sixth)-bounds[i]);
511 rav = max(0.0, mypow(rav, -sixth)-bounds[i]);
520 static void analyse_disre(const char *voutfn, int nframes,
521 real violaver[], real bounds[], int index[],
522 int pair[], int nbounds,
523 const output_env_t oenv)
526 double sum, sumt, sumaver;
529 /* Subtract bounds from distances, to calculate violations */
530 calc_violations(violaver, violaver,
531 nbounds, pair, bounds, NULL, &sumt, &sumaver);
534 fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
536 fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n",
539 vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm",
543 for (i = 0; (i < nbounds); i++)
545 /* Do ensemble averaging */
547 for (j = pair[i]; (j < pair[i+1]); j++)
549 sumaver += sqr(violaver[j]/nframes);
551 sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]);
554 sum = max(sum, sumaver);
555 fprintf(vout, "%10d %10.5e\n", index[i], sumaver);
558 for (j = 0; (j < dr.ndr); j++)
560 fprintf(vout, "%10d %10.5e\n", j, mypow(violaver[j]/nframes, minthird));
565 fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
567 fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
569 do_view(oenv, voutfn, "-graphtype bar");
572 static void einstein_visco(const char *fn, const char *fni, int nsets,
573 int nframes, real **sum,
574 real V, real T, int nsteps, double time[],
575 const output_env_t oenv)
578 double av[4], avold[4];
587 dt = (time[1]-time[0]);
590 for (i = 0; i <= nsets; i++)
594 fp0 = xvgropen(fni, "Shear viscosity integral",
595 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
596 fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
597 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
598 for (i = 1; i < nf4; i++)
600 fac = dt*nframes/nsteps;
601 for (m = 0; m <= nsets; m++)
605 for (j = 0; j < nframes-i; j++)
607 for (m = 0; m < nsets; m++)
609 di = sqr(fac*(sum[m][j+i]-sum[m][j]));
612 av[nsets] += di/nsets;
615 /* Convert to SI for the viscosity */
616 fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i);
617 fprintf(fp0, "%10g", time[i]-time[0]);
618 for (m = 0; (m <= nsets); m++)
621 fprintf(fp0, " %10g", av[m]);
624 fprintf(fp1, "%10g", 0.5*(time[i]+time[i-1])-time[0]);
625 for (m = 0; (m <= nsets); m++)
627 fprintf(fp1, " %10g", (av[m]-avold[m])/dt);
650 static void clear_ee_sum(ee_sum_t *ees)
658 static void add_ee_sum(ee_sum_t *ees, double sum, int np)
664 static void add_ee_av(ee_sum_t *ees)
668 av = ees->sum/ees->np;
675 static double calc_ee2(int nb, ee_sum_t *ees)
677 return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
680 static void set_ee_av(ener_ee_t *eee)
684 char buf[STEPSTRSIZE];
685 fprintf(debug, "Storing average for err.est.: %s steps\n",
686 gmx_step_str(eee->nst, buf));
688 add_ee_av(&eee->sum);
690 if (eee->b == 1 || eee->nst < eee->nst_min)
692 eee->nst_min = eee->nst;
697 static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
700 double sum, sum2, sump, see2;
701 gmx_int64_t steps, np, p, bound_nb;
705 double x, sx, sy, sxx, sxy;
708 /* Check if we have exact statistics over all points */
709 for (i = 0; i < nset; i++)
712 ed->bExactStat = FALSE;
713 if (edat->npoints > 0)
715 /* All energy file sum entries 0 signals no exact sums.
716 * But if all energy values are 0, we still have exact sums.
719 for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
721 if (ed->ener[i] != 0)
725 ed->bExactStat = (ed->es[f].sum != 0);
729 ed->bExactStat = TRUE;
735 for (i = 0; i < nset; i++)
746 for (nb = nbmin; nb <= nbmax; nb++)
749 clear_ee_sum(&eee[nb].sum);
753 for (f = 0; f < edat->nframes; f++)
759 /* Add the sum and the sum of variances to the totals. */
765 sum2 += dsqr(sum/np - (sum + es->sum)/(np + p))
771 /* Add a single value to the sum and sum of squares. */
777 /* sum has to be increased after sum2 */
781 /* For the linear regression use variance 1/p.
782 * Note that sump is the sum, not the average, so we don't need p*.
784 x = edat->step[f] - 0.5*(edat->steps[f] - 1);
790 for (nb = nbmin; nb <= nbmax; nb++)
792 /* Check if the current end step is closer to the desired
793 * block boundary than the next end step.
795 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
796 if (eee[nb].nst > 0 &&
797 bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
807 eee[nb].nst += edat->step[f] - edat->step[f-1];
811 add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
815 add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
817 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
818 if (edat->step[f]*nb >= bound_nb)
825 edat->s[i].av = sum/np;
828 edat->s[i].rmsd = sqrt(sum2/np);
832 edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av));
835 if (edat->nframes > 1)
837 edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
841 edat->s[i].slope = 0;
846 for (nb = nbmin; nb <= nbmax; nb++)
848 /* Check if we actually got nb blocks and if the smallest
849 * block is not shorter than 80% of the average.
853 char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
854 fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
856 gmx_step_str(eee[nb].nst_min, buf1),
857 gmx_step_str(edat->nsteps, buf2));
859 if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
861 see2 += calc_ee2(nb, &eee[nb].sum);
867 edat->s[i].ee = sqrt(see2/nee);
877 static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax)
888 snew(s->ener, esum->nframes);
889 snew(s->es, esum->nframes);
891 s->bExactStat = TRUE;
893 for (i = 0; i < nset; i++)
895 if (!edat->s[i].bExactStat)
897 s->bExactStat = FALSE;
899 s->slope += edat->s[i].slope;
902 for (f = 0; f < edat->nframes; f++)
905 for (i = 0; i < nset; i++)
907 sum += edat->s[i].ener[f];
911 for (i = 0; i < nset; i++)
913 sum += edat->s[i].es[f].sum;
919 calc_averages(1, esum, nbmin, nbmax);
924 static char *ee_pr(double ee, char *buf)
931 sprintf(buf, "%s", "--");
935 /* Round to two decimals by printing. */
936 sprintf(tmp, "%.1e", ee);
937 sscanf(tmp, "%lf", &rnd);
938 sprintf(buf, "%g", rnd);
944 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat)
946 /* Remove the drift by performing a fit to y = ax+b.
947 Uses 5 iterations. */
949 double delta, d, sd, sd2;
951 edat->npoints = edat->nframes;
952 edat->nsteps = edat->nframes;
954 for (k = 0; (k < 5); k++)
956 for (i = 0; (i < nset); i++)
958 delta = edat->s[i].slope*dt;
962 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
965 for (j = 0; (j < edat->nframes); j++)
967 edat->s[i].ener[j] -= j*delta;
968 edat->s[i].es[j].sum = 0;
969 edat->s[i].es[j].sum2 = 0;
972 calc_averages(nset, edat, nbmin, nbmax);
976 static void calc_fluctuation_props(FILE *fp,
977 gmx_bool bDriftCorr, real dt,
979 char **leg, enerdata_t *edat,
980 int nbmin, int nbmax)
983 double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
986 eVol, eEnth, eTemp, eEtot, eNR
988 const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
991 NANO3 = NANO*NANO*NANO;
994 fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n"
995 "for spurious drift in the graphs. Note that this is not\n"
996 "a substitute for proper equilibration and sampling!\n");
1000 remove_drift(nset, nbmin, nbmax, dt, edat);
1002 for (i = 0; (i < eNR); i++)
1004 for (ii[i] = 0; (ii[i] < nset &&
1005 (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++)
1009 /* if (ii[i] < nset)
1010 fprintf(fp,"Found %s data.\n",my_ener[i]);
1012 /* Compute it all! */
1013 alpha = kappa = cp = dcp = cv = NOTSET;
1017 if (ii[eTemp] < nset)
1019 tt = edat->s[ii[eTemp]].av;
1023 if ((ii[eVol] < nset) && (ii[eTemp] < nset))
1025 vv = edat->s[ii[eVol]].av*NANO3;
1026 varv = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
1027 kappa = (varv/vv)/(BOLTZMANN*tt);
1031 if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
1033 hh = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
1034 varh = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
1035 cp = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
1038 et = varet = NOTSET;
1039 if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
1041 /* Only compute cv in constant volume runs, which we can test
1042 by checking whether the enthalpy was computed.
1044 et = edat->s[ii[eEtot]].av;
1045 varet = sqr(edat->s[ii[eEtot]].rmsd);
1046 cv = KILO*((varet/nmol)/(BOLTZ*tt*tt));
1049 if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
1051 double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
1052 vh_sum = v_sum = h_sum = 0;
1053 for (j = 0; (j < edat->nframes); j++)
1055 v = edat->s[ii[eVol]].ener[j]*NANO3;
1056 h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
1061 vh_aver = vh_sum / edat->nframes;
1062 v_aver = v_sum / edat->nframes;
1063 h_aver = h_sum / edat->nframes;
1064 alpha = (vh_aver-v_aver*h_aver)/(v_aver*BOLTZMANN*tt*tt);
1065 dcp = (v_aver*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
1072 fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n",
1075 fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
1076 fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
1077 fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
1078 fprintf(fp, "please use the g_dos program.\n\n");
1079 fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n"
1080 "a block-averaging error analysis (not implemented in g_energy yet)\n");
1086 fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol);
1091 fprintf(fp, "Volume = %10g m^3/mol\n",
1096 fprintf(fp, "Enthalpy = %10g kJ/mol\n",
1097 hh*AVOGADRO/(KILO*nmol));
1099 if (alpha != NOTSET)
1101 fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
1104 if (kappa != NOTSET)
1106 fprintf(fp, "Isothermal Compressibility Kappa = %10g (J/m^3)\n",
1108 fprintf(fp, "Adiabatic bulk modulus = %10g (m^3/J)\n",
1113 fprintf(fp, "Heat capacity at constant pressure Cp = %10g J/mol K\n",
1118 fprintf(fp, "Heat capacity at constant volume Cv = %10g J/mol K\n",
1123 fprintf(fp, "Cp-Cv = %10g J/mol K\n",
1126 please_cite(fp, "Allen1987a");
1130 fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
1134 static void analyse_ener(gmx_bool bCorr, const char *corrfn,
1135 gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct,
1136 gmx_bool bVisco, const char *visfn, int nmol,
1137 gmx_int64_t start_step, double start_t,
1138 gmx_int64_t step, double t,
1139 double time[], real reftemp,
1141 int nset, int set[], gmx_bool *bIsEner,
1142 char **leg, gmx_enxnm_t *enm,
1143 real Vaver, real ezero,
1144 int nbmin, int nbmax,
1145 const output_env_t oenv)
1148 /* Check out the printed manual for equations! */
1149 double Dt, aver, stddev, errest, delta_t, totaldrift;
1150 enerdata_t *esum = NULL;
1151 real xxx, integral, intBulk, Temp = 0, Pres = 0;
1152 real sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest;
1153 double beta = 0, expE, expEtot, *fee = NULL;
1155 int nexact, nnotexact;
1159 char buf[256], eebuf[100];
1161 nsteps = step - start_step + 1;
1164 fprintf(stdout, "Not enough steps (%s) for statistics\n",
1165 gmx_step_str(nsteps, buf));
1169 /* Calculate the time difference */
1170 delta_t = t - start_t;
1172 fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
1173 gmx_step_str(nsteps, buf), start_t, t, nset);
1175 calc_averages(nset, edat, nbmin, nbmax);
1179 esum = calc_sum(nset, edat, nbmin, nbmax);
1182 if (edat->npoints == 0)
1191 for (i = 0; (i < nset); i++)
1193 if (edat->s[i].bExactStat)
1206 fprintf(stdout, "All statistics are over %s points\n",
1207 gmx_step_str(edat->npoints, buf));
1209 else if (nexact == 0 || edat->npoints == edat->nframes)
1211 fprintf(stdout, "All statistics are over %d points (frames)\n",
1216 fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
1217 for (i = 0; (i < nset); i++)
1219 if (!edat->s[i].bExactStat)
1221 fprintf(stdout, " '%s'", leg[i]);
1224 fprintf(stdout, " %s has statistics over %d points (frames)\n",
1225 nnotexact == 1 ? "is" : "are", edat->nframes);
1226 fprintf(stdout, "All other statistics are over %s points\n",
1227 gmx_step_str(edat->npoints, buf));
1229 fprintf(stdout, "\n");
1231 fprintf(stdout, "%-24s %10s %10s %10s %10s",
1232 "Energy", "Average", "Err.Est.", "RMSD", "Tot-Drift");
1235 fprintf(stdout, " %10s\n", "-kT ln<e^(E/kT)>");
1239 fprintf(stdout, "\n");
1241 fprintf(stdout, "-------------------------------------------------------------------------------\n");
1243 /* Initiate locals, only used with -sum */
1247 beta = 1.0/(BOLTZ*reftemp);
1250 for (i = 0; (i < nset); i++)
1252 aver = edat->s[i].av;
1253 stddev = edat->s[i].rmsd;
1254 errest = edat->s[i].ee;
1259 for (j = 0; (j < edat->nframes); j++)
1261 expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
1265 expEtot += expE/edat->nframes;
1268 fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
1270 if (strstr(leg[i], "empera") != NULL)
1274 else if (strstr(leg[i], "olum") != NULL)
1278 else if (strstr(leg[i], "essure") != NULL)
1284 pr_aver = aver/nmol-ezero;
1285 pr_stddev = stddev/nmol;
1286 pr_errest = errest/nmol;
1295 /* Multiply the slope in steps with the number of steps taken */
1296 totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
1302 fprintf(stdout, "%-24s %10g %10s %10g %10g",
1303 leg[i], pr_aver, ee_pr(pr_errest, eebuf), pr_stddev, totaldrift);
1306 fprintf(stdout, " %10g", fee[i]);
1309 fprintf(stdout, " (%s)\n", enm[set[i]].unit);
1313 for (j = 0; (j < edat->nframes); j++)
1315 edat->s[i].ener[j] -= aver;
1321 totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1322 fprintf(stdout, "%-24s %10g %10s %10s %10g (%s)",
1323 "Total", esum->s[0].av/nmol, ee_pr(esum->s[0].ee/nmol, eebuf),
1324 "--", totaldrift/nmol, enm[set[0]].unit);
1325 /* pr_aver,pr_stddev,a,totaldrift */
1328 fprintf(stdout, " %10g %10g\n",
1329 log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta);
1333 fprintf(stdout, "\n");
1337 /* Do correlation function */
1338 if (edat->nframes > 1)
1340 Dt = delta_t/(edat->nframes - 1);
1348 const char* leg[] = { "Shear", "Bulk" };
1353 /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1355 /* Symmetrise tensor! (and store in first three elements)
1356 * And subtract average pressure!
1359 for (i = 0; i < 12; i++)
1361 snew(eneset[i], edat->nframes);
1364 for (i = 0; i < 3; i++)
1366 snew(enesum[i], edat->nframes);
1368 for (i = 0; (i < edat->nframes); i++)
1370 eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1371 eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1372 eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1373 for (j = 3; j <= 11; j++)
1375 eneset[j][i] = edat->s[j].ener[i];
1377 eneset[11][i] -= Pres;
1378 enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum);
1379 enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
1380 enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
1383 einstein_visco("evisco.xvg", "eviscoi.xvg",
1384 3, edat->nframes, enesum, Vaver, Temp, nsteps, time, oenv);
1386 /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1387 /* Do it for shear viscosity */
1388 strcpy(buf, "Shear Viscosity");
1389 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3,
1390 (edat->nframes+1)/2, eneset, Dt,
1391 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1393 /* Now for bulk viscosity */
1394 strcpy(buf, "Bulk Viscosity");
1395 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1,
1396 (edat->nframes+1)/2, &(eneset[11]), Dt,
1397 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1399 factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1400 fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1401 xvgr_legend(fp, asize(leg), leg, oenv);
1403 /* Use trapezium rule for integration */
1406 nout = get_acfnout();
1407 if ((nout < 2) || (nout >= edat->nframes/2))
1409 nout = edat->nframes/2;
1411 for (i = 1; (i < nout); i++)
1413 integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor;
1414 intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1415 fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk);
1423 strcpy(buf, "Autocorrelation of Energy Fluctuations");
1427 strcpy(buf, "Energy Autocorrelation");
1430 do_autocorr(corrfn, oenv, buf, edat->nframes,
1432 bSum ? &edat->s[nset-1].ener : eneset,
1433 (delta_t/edat->nframes), eacNormal, FALSE);
1439 static void print_time(FILE *fp, double t)
1441 fprintf(fp, "%12.6f", t);
1444 static void print1(FILE *fp, gmx_bool bDp, real e)
1448 fprintf(fp, " %16.12f", e);
1452 fprintf(fp, " %10.6f", e);
1456 static void fec(const char *ene2fn, const char *runavgfn,
1457 real reftemp, int nset, int set[], char *leg[],
1458 enerdata_t *edat, double time[],
1459 const output_env_t oenv)
1461 const char * ravgleg[] = {
1462 "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1463 "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N"
1467 int nre, timecheck, step, nenergy, nenergy2, maxenergy;
1473 gmx_enxnm_t *enm = NULL;
1477 /* read second energy file */
1480 enx = open_enx(ene2fn, "r");
1481 do_enxnms(enx, &(fr->nre), &enm);
1483 snew(eneset2, nset+1);
1489 /* This loop searches for the first frame (when -b option is given),
1490 * or when this has been found it reads just one energy frame
1494 bCont = do_enx(enx, fr);
1498 timecheck = check_times(fr->t);
1502 while (bCont && (timecheck < 0));
1504 /* Store energies for analysis afterwards... */
1505 if ((timecheck == 0) && bCont)
1509 if (nenergy2 >= maxenergy)
1512 for (i = 0; i <= nset; i++)
1514 srenew(eneset2[i], maxenergy);
1517 if (fr->t != time[nenergy2])
1519 fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n",
1520 fr->t, time[nenergy2], gmx_step_str(fr->step, buf));
1522 for (i = 0; i < nset; i++)
1524 eneset2[i][nenergy2] = fr->ener[set[i]].e;
1530 while (bCont && (timecheck == 0));
1533 if (edat->nframes != nenergy2)
1535 fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n",
1536 edat->nframes, nenergy2);
1538 nenergy = min(edat->nframes, nenergy2);
1540 /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1544 fp = xvgropen(runavgfn, "Running average free energy difference",
1545 "Time (" unit_time ")", "\\8D\\4E (" unit_energy ")", oenv);
1546 xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
1548 fprintf(stdout, "\n%-24s %10s\n",
1549 "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
1551 beta = 1.0/(BOLTZ*reftemp);
1552 for (i = 0; i < nset; i++)
1554 if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
1556 fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n",
1557 leg[i], enm[set[i]].name);
1559 for (j = 0; j < nenergy; j++)
1561 dE = eneset2[i][j] - edat->s[i].ener[j];
1562 sum += exp(-dE*beta);
1565 fprintf(fp, "%10g %10g %10g\n",
1566 time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
1569 aver = -BOLTZ*reftemp*log(sum/nenergy);
1570 fprintf(stdout, "%-24s %10g\n", leg[i], aver);
1580 static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl,
1581 const char *filename, gmx_bool bDp,
1582 int *blocks, int *hists, int *samples, int *nlambdas,
1583 const output_env_t oenv)
1585 const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1586 char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
1588 gmx_bool first = FALSE;
1589 int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
1592 double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0;
1593 static int setnr = 0;
1594 double *native_lambda_vec = NULL;
1595 const char **lambda_components = NULL;
1596 int n_lambda_vec = 0;
1597 gmx_bool changing_lambda = FALSE;
1598 int lambda_fep_state;
1600 /* now count the blocks & handle the global dh data */
1601 for (i = 0; i < fr->nblock; i++)
1603 if (fr->block[i].id == enxDHHIST)
1607 else if (fr->block[i].id == enxDH)
1611 else if (fr->block[i].id == enxDHCOLL)
1614 if ( (fr->block[i].nsub < 1) ||
1615 (fr->block[i].sub[0].type != xdr_datatype_double) ||
1616 (fr->block[i].sub[0].nr < 5))
1618 gmx_fatal(FARGS, "Unexpected block data");
1621 /* read the data from the DHCOLL block */
1622 temp = fr->block[i].sub[0].dval[0];
1623 start_time = fr->block[i].sub[0].dval[1];
1624 delta_time = fr->block[i].sub[0].dval[2];
1625 start_lambda = fr->block[i].sub[0].dval[3];
1626 delta_lambda = fr->block[i].sub[0].dval[4];
1627 changing_lambda = (delta_lambda != 0);
1628 if (fr->block[i].nsub > 1)
1630 lambda_fep_state = fr->block[i].sub[1].ival[0];
1631 if (n_lambda_vec == 0)
1633 n_lambda_vec = fr->block[i].sub[1].ival[1];
1637 if (n_lambda_vec != fr->block[i].sub[1].ival[1])
1640 "Unexpected change of basis set in lambda");
1643 if (lambda_components == NULL)
1645 snew(lambda_components, n_lambda_vec);
1647 if (native_lambda_vec == NULL)
1649 snew(native_lambda_vec, n_lambda_vec);
1651 for (j = 0; j < n_lambda_vec; j++)
1653 native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j];
1654 lambda_components[j] =
1655 efpt_singular_names[fr->block[i].sub[1].ival[2+j]];
1661 if (nblock_hist == 0 && nblock_dh == 0)
1663 /* don't do anything */
1666 if (nblock_hist > 0 && nblock_dh > 0)
1668 gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1674 /* we have standard, non-histogram data --
1675 call open_dhdl to open the file */
1676 /* TODO this is an ugly hack that needs to be fixed: this will only
1677 work if the order of data is always the same and if we're
1678 only using the g_energy compiled with the mdrun that produced
1680 *fp_dhdl = open_dhdl(filename, ir, oenv);
1684 sprintf(title, "N(%s)", deltag);
1685 sprintf(label_x, "%s (%s)", deltag, unit_energy);
1686 sprintf(label_y, "Samples");
1687 *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1688 sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
1689 xvgr_subtitle(*fp_dhdl, buf, oenv);
1693 (*hists) += nblock_hist;
1694 (*blocks) += nblock_dh;
1695 (*nlambdas) = nblock_hist+nblock_dh;
1697 /* write the data */
1698 if (nblock_hist > 0)
1700 gmx_int64_t sum = 0;
1702 for (i = 0; i < fr->nblock; i++)
1704 t_enxblock *blk = &(fr->block[i]);
1705 if (blk->id == enxDHHIST)
1707 double foreign_lambda, dx;
1709 int nhist, derivative;
1711 /* check the block types etc. */
1712 if ( (blk->nsub < 2) ||
1713 (blk->sub[0].type != xdr_datatype_double) ||
1714 (blk->sub[1].type != xdr_datatype_int64) ||
1715 (blk->sub[0].nr < 2) ||
1716 (blk->sub[1].nr < 2) )
1718 gmx_fatal(FARGS, "Unexpected block data in file");
1720 foreign_lambda = blk->sub[0].dval[0];
1721 dx = blk->sub[0].dval[1];
1722 nhist = blk->sub[1].lval[0];
1723 derivative = blk->sub[1].lval[1];
1724 for (j = 0; j < nhist; j++)
1727 x0 = blk->sub[1].lval[2+j];
1731 sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1732 deltag, lambda, foreign_lambda,
1733 lambda, start_lambda);
1737 sprintf(legend, "N(%s | %s=%g)",
1738 dhdl, lambda, start_lambda);
1742 xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1744 for (k = 0; k < blk->sub[j+2].nr; k++)
1749 hist = blk->sub[j+2].ival[k];
1752 fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist,
1756 /* multiple histogram data blocks in one histogram
1757 mean that the second one is the reverse of the first one:
1758 for dhdl derivatives, it's important to know both the
1759 maximum and minimum values */
1764 (*samples) += (int)(sum/nblock_hist);
1770 char **setnames = NULL;
1771 int nnames = nblock_dh;
1773 for (i = 0; i < fr->nblock; i++)
1775 t_enxblock *blk = &(fr->block[i]);
1776 if (blk->id == enxDH)
1780 len = blk->sub[2].nr;
1784 if (len != blk->sub[2].nr)
1786 gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1793 for (i = 0; i < len; i++)
1795 double time = start_time + delta_time*i;
1797 fprintf(*fp_dhdl, "%.4f ", time);
1799 for (j = 0; j < fr->nblock; j++)
1801 t_enxblock *blk = &(fr->block[j]);
1802 if (blk->id == enxDH)
1805 if (blk->sub[2].type == xdr_datatype_float)
1807 value = blk->sub[2].fval[i];
1811 value = blk->sub[2].dval[i];
1813 /* we need to decide which data type it is based on the count*/
1815 if (j == 1 && ir->bExpanded)
1817 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! */
1823 fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */
1827 fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */
1832 fprintf(*fp_dhdl, "\n");
1838 int gmx_energy(int argc, char *argv[])
1840 const char *desc[] = {
1841 "[THISMODULE] extracts energy components or distance restraint",
1842 "data from an energy file. The user is prompted to interactively",
1843 "select the desired energy terms.[PAR]",
1845 "Average, RMSD, and drift are calculated with full precision from the",
1846 "simulation (see printed manual). Drift is calculated by performing",
1847 "a least-squares fit of the data to a straight line. The reported total drift",
1848 "is the difference of the fit at the first and last point.",
1849 "An error estimate of the average is given based on a block averages",
1850 "over 5 blocks using the full-precision averages. The error estimate",
1851 "can be performed over multiple block lengths with the options",
1852 "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1853 "[BB]Note[bb] that in most cases the energy files contains averages over all",
1854 "MD steps, or over many more points than the number of frames in",
1855 "energy file. This makes the [THISMODULE] statistics output more accurate",
1856 "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1857 "file, the statistics mentioned above are simply over the single, per-frame",
1858 "energy values.[PAR]",
1860 "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1862 "Some fluctuation-dependent properties can be calculated provided",
1863 "the correct energy terms are selected, and that the command line option",
1864 "[TT]-fluct_props[tt] is given. The following properties",
1865 "will be computed:[BR]",
1866 "Property Energy terms needed[BR]",
1867 "---------------------------------------------------[BR]",
1868 "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp [BR]",
1869 "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp [BR]",
1870 "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1871 "Isothermal compressibility: Vol, Temp [BR]",
1872 "Adiabatic bulk modulus: Vol, Temp [BR]",
1873 "---------------------------------------------------[BR]",
1874 "You always need to set the number of molecules [TT]-nmol[tt].",
1875 "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1876 "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]"
1877 "When the [TT]-viol[tt] option is set, the time averaged",
1878 "violations are plotted and the running time-averaged and",
1879 "instantaneous sum of violations are recalculated. Additionally",
1880 "running time-averaged and instantaneous distances between",
1881 "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1883 "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1884 "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1885 "The first two options plot the orientation, the last three the",
1886 "deviations of the orientations from the experimental values.",
1887 "The options that end on an 'a' plot the average over time",
1888 "as a function of restraint. The options that end on a 't'",
1889 "prompt the user for restraint label numbers and plot the data",
1890 "as a function of time. Option [TT]-odr[tt] plots the RMS",
1891 "deviation as a function of restraint.",
1892 "When the run used time or ensemble averaged orientation restraints,",
1893 "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1894 "not ensemble-averaged orientations and deviations instead of",
1895 "the time and ensemble averages.[PAR]",
1897 "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1898 "tensor for each orientation restraint experiment. With option",
1899 "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1901 "Option [TT]-odh[tt] extracts and plots the free energy data",
1902 "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1903 "from the [TT]ener.edr[tt] file.[PAR]",
1905 "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1906 "difference with an ideal gas state: [BR]",
1907 " [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]",
1908 " [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]",
1909 "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1910 "the average is over the ensemble (or time in a trajectory).",
1911 "Note that this is in principle",
1912 "only correct when averaging over the whole (Boltzmann) ensemble",
1913 "and using the potential energy. This also allows for an entropy",
1914 "estimate using:[BR]",
1915 " [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]",
1916 " [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",
1919 "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1920 "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1921 "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1922 "files, and the average is over the ensemble A. The running average",
1923 "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1924 "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1927 static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1928 static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1929 static int skip = 0, nmol = 1, nbmin = 5, nbmax = 5;
1930 static real reftemp = 300.0, ezero = 0;
1932 { "-fee", FALSE, etBOOL, {&bFee},
1933 "Do a free energy estimate" },
1934 { "-fetemp", FALSE, etREAL, {&reftemp},
1935 "Reference temperature for free energy calculation" },
1936 { "-zero", FALSE, etREAL, {&ezero},
1937 "Subtract a zero-point energy" },
1938 { "-sum", FALSE, etBOOL, {&bSum},
1939 "Sum the energy terms selected rather than display them all" },
1940 { "-dp", FALSE, etBOOL, {&bDp},
1941 "Print energies in high precision" },
1942 { "-nbmin", FALSE, etINT, {&nbmin},
1943 "Minimum number of blocks for error estimate" },
1944 { "-nbmax", FALSE, etINT, {&nbmax},
1945 "Maximum number of blocks for error estimate" },
1946 { "-mutot", FALSE, etBOOL, {&bMutot},
1947 "Compute the total dipole moment from the components" },
1948 { "-skip", FALSE, etINT, {&skip},
1949 "Skip number of frames between data points" },
1950 { "-aver", FALSE, etBOOL, {&bPrAll},
1951 "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1952 { "-nmol", FALSE, etINT, {&nmol},
1953 "Number of molecules in your sample: the energies are divided by this number" },
1954 { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1955 "Compute properties based on energy fluctuations, like heat capacity" },
1956 { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1957 "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1958 { "-fluc", FALSE, etBOOL, {&bFluct},
1959 "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1960 { "-orinst", FALSE, etBOOL, {&bOrinst},
1961 "Analyse instantaneous orientation data" },
1962 { "-ovec", FALSE, etBOOL, {&bOvec},
1963 "Also plot the eigenvectors with [TT]-oten[tt]" }
1965 const char * drleg[] = {
1969 static const char *setnm[] = {
1970 "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1971 "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1972 "Volume", "Pressure"
1975 FILE *out = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL;
1976 FILE *fp_dhdl = NULL;
1981 gmx_localtop_t *top = NULL;
1985 gmx_enxnm_t *enm = NULL;
1986 t_enxframe *frame, *fr = NULL;
1988 #define NEXT (1-cur)
1989 int nre, teller, teller_disre, nfr;
1990 gmx_int64_t start_step;
1991 int nor = 0, nex = 0, norfr = 0, enx_i = 0;
1993 real *bounds = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL;
1994 int *index = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL;
1995 int nbounds = 0, npairs;
1996 gmx_bool bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL;
1997 gmx_bool bFoundStart, bCont, bEDR, bVisco;
1998 double sum, sumaver, sumt, ener, dbl;
1999 double *time = NULL;
2001 int *set = NULL, i, j, k, nset, sss;
2002 gmx_bool *bIsEner = NULL;
2003 char **pairleg, **odtleg, **otenleg;
2006 char *anm_j, *anm_k, *resnm_j, *resnm_k;
2007 int resnr_j, resnr_k;
2008 const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
2011 t_enxblock *blk = NULL;
2012 t_enxblock *blk_disre = NULL;
2014 int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
2017 { efEDR, "-f", NULL, ffREAD },
2018 { efEDR, "-f2", NULL, ffOPTRD },
2019 { efTPX, "-s", NULL, ffOPTRD },
2020 { efXVG, "-o", "energy", ffWRITE },
2021 { efXVG, "-viol", "violaver", ffOPTWR },
2022 { efXVG, "-pairs", "pairs", ffOPTWR },
2023 { efXVG, "-ora", "orienta", ffOPTWR },
2024 { efXVG, "-ort", "orientt", ffOPTWR },
2025 { efXVG, "-oda", "orideva", ffOPTWR },
2026 { efXVG, "-odr", "oridevr", ffOPTWR },
2027 { efXVG, "-odt", "oridevt", ffOPTWR },
2028 { efXVG, "-oten", "oriten", ffOPTWR },
2029 { efXVG, "-corr", "enecorr", ffOPTWR },
2030 { efXVG, "-vis", "visco", ffOPTWR },
2031 { efXVG, "-ravg", "runavgdf", ffOPTWR },
2032 { efXVG, "-odh", "dhdl", ffOPTWR }
2034 #define NFILE asize(fnm)
2039 ppa = add_acf_pargs(&npargs, pa);
2040 if (!parse_common_args(&argc, argv,
2041 PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
2042 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
2047 bDRAll = opt2bSet("-pairs", NFILE, fnm);
2048 bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
2049 bORA = opt2bSet("-ora", NFILE, fnm);
2050 bORT = opt2bSet("-ort", NFILE, fnm);
2051 bODA = opt2bSet("-oda", NFILE, fnm);
2052 bODR = opt2bSet("-odr", NFILE, fnm);
2053 bODT = opt2bSet("-odt", NFILE, fnm);
2054 bORIRE = bORA || bORT || bODA || bODR || bODT;
2055 bOTEN = opt2bSet("-oten", NFILE, fnm);
2056 bDHDL = opt2bSet("-odh", NFILE, fnm);
2061 fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
2062 do_enxnms(fp, &nre, &enm);
2066 bVisco = opt2bSet("-vis", NFILE, fnm);
2068 if ((!bDisRe) && (!bDHDL))
2072 nset = asize(setnm);
2074 /* This is nasty code... To extract Pres tensor, Volume and Temperature */
2075 for (j = 0; j < nset; j++)
2077 for (i = 0; i < nre; i++)
2079 if (strstr(enm[i].name, setnm[j]))
2087 if (gmx_strcasecmp(setnm[j], "Volume") == 0)
2089 printf("Enter the box volume (" unit_volume "): ");
2090 if (1 != scanf("%lf", &dbl))
2092 gmx_fatal(FARGS, "Error reading user input");
2098 gmx_fatal(FARGS, "Could not find term %s for viscosity calculation",
2106 set = select_by_name(nre, enm, &nset);
2108 /* Print all the different units once */
2109 sprintf(buf, "(%s)", enm[set[0]].unit);
2110 for (i = 1; i < nset; i++)
2112 for (j = 0; j < i; j++)
2114 if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
2122 strcat(buf, enm[set[i]].unit);
2126 out = xvgropen(opt2fn("-o", NFILE, fnm), "Gromacs Energies", "Time (ps)", buf,
2130 for (i = 0; (i < nset); i++)
2132 leg[i] = enm[set[i]].name;
2136 leg[nset] = strdup("Sum");
2137 xvgr_legend(out, nset+1, (const char**)leg, oenv);
2141 xvgr_legend(out, nset, (const char**)leg, oenv);
2144 snew(bIsEner, nset);
2145 for (i = 0; (i < nset); i++)
2148 for (j = 0; (j <= F_ETOT); j++)
2150 bIsEner[i] = bIsEner[i] ||
2151 (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
2155 if (bPrAll && nset > 1)
2157 gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
2162 if (bORIRE || bOTEN)
2164 get_orires_parms(ftp2fn(efTPX, NFILE, fnm), &nor, &nex, &or_label, &oobs);
2188 fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n");
2189 fprintf(stderr, "End your selection with 0\n");
2196 if (1 != scanf("%d", &(orsel[j])))
2198 gmx_fatal(FARGS, "Error reading user input");
2201 while (orsel[j] > 0);
2204 fprintf(stderr, "Selecting all %d orientation restraints\n", nor);
2207 for (i = 0; i < nor; i++)
2214 /* Build the selection */
2216 for (i = 0; i < j; i++)
2218 for (k = 0; k < nor; k++)
2220 if (or_label[k] == orsel[i])
2229 fprintf(stderr, "Orientation restraint label %d not found\n",
2234 snew(odtleg, norsel);
2235 for (i = 0; i < norsel; i++)
2237 snew(odtleg[i], 256);
2238 sprintf(odtleg[i], "%d", or_label[orsel[i]]);
2242 fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
2243 "Time (ps)", "", oenv);
2246 fprintf(fort, "%s", orinst_sub);
2248 xvgr_legend(fort, norsel, (const char**)odtleg, oenv);
2252 fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
2253 "Orientation restraint deviation",
2254 "Time (ps)", "", oenv);
2257 fprintf(fodt, "%s", orinst_sub);
2259 xvgr_legend(fodt, norsel, (const char**)odtleg, oenv);
2265 foten = xvgropen(opt2fn("-oten", NFILE, fnm),
2266 "Order tensor", "Time (ps)", "", oenv);
2267 snew(otenleg, bOvec ? nex*12 : nex*3);
2268 for (i = 0; i < nex; i++)
2270 for (j = 0; j < 3; j++)
2272 sprintf(buf, "eig%d", j+1);
2273 otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
2277 for (j = 0; j < 9; j++)
2279 sprintf(buf, "vec%d%s", j/3+1, j%3 == 0 ? "x" : (j%3 == 1 ? "y" : "z"));
2280 otenleg[12*i+3+j] = strdup(buf);
2284 xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv);
2289 nbounds = get_bounds(ftp2fn(efTPX, NFILE, fnm), &bounds, &index, &pair, &npairs,
2291 snew(violaver, npairs);
2292 out = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations",
2293 "Time (ps)", "nm", oenv);
2294 xvgr_legend(out, 2, drleg, oenv);
2297 fp_pairs = xvgropen(opt2fn("-pairs", NFILE, fnm), "Pair Distances",
2298 "Time (ps)", "Distance (nm)", oenv);
2299 if (output_env_get_print_xvgr_codes(oenv))
2301 fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
2308 get_dhdl_parms(ftp2fn(efTPX, NFILE, fnm), &ir);
2311 /* Initiate energies and set them to zero */
2320 /* Initiate counters */
2323 bFoundStart = FALSE;
2328 /* This loop searches for the first frame (when -b option is given),
2329 * or when this has been found it reads just one energy frame
2333 bCont = do_enx(fp, &(frame[NEXT]));
2336 timecheck = check_times(frame[NEXT].t);
2339 while (bCont && (timecheck < 0));
2341 if ((timecheck == 0) && bCont)
2343 /* We read a valid frame, so we can use it */
2344 fr = &(frame[NEXT]);
2348 /* The frame contains energies, so update cur */
2351 if (edat.nframes % 1000 == 0)
2353 srenew(edat.step, edat.nframes+1000);
2354 memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
2355 srenew(edat.steps, edat.nframes+1000);
2356 memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
2357 srenew(edat.points, edat.nframes+1000);
2358 memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
2360 for (i = 0; i < nset; i++)
2362 srenew(edat.s[i].ener, edat.nframes+1000);
2363 memset(&(edat.s[i].ener[edat.nframes]), 0,
2364 1000*sizeof(edat.s[i].ener[0]));
2365 srenew(edat.s[i].es, edat.nframes+1000);
2366 memset(&(edat.s[i].es[edat.nframes]), 0,
2367 1000*sizeof(edat.s[i].es[0]));
2372 edat.step[nfr] = fr->step;
2377 /* Initiate the previous step data */
2378 start_step = fr->step;
2380 /* Initiate the energy sums */
2381 edat.steps[nfr] = 1;
2382 edat.points[nfr] = 1;
2383 for (i = 0; i < nset; i++)
2386 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2387 edat.s[i].es[nfr].sum2 = 0;
2394 edat.steps[nfr] = fr->nsteps;
2396 if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2400 edat.points[nfr] = 1;
2401 for (i = 0; i < nset; i++)
2404 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2405 edat.s[i].es[nfr].sum2 = 0;
2411 edat.points[nfr] = fr->nsum;
2412 for (i = 0; i < nset; i++)
2415 edat.s[i].es[nfr].sum = fr->ener[sss].esum;
2416 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2418 edat.npoints += fr->nsum;
2423 /* The interval does not match fr->nsteps:
2424 * can not do exact averages.
2428 edat.nsteps = fr->step - start_step + 1;
2431 for (i = 0; i < nset; i++)
2433 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2437 * Define distance restraint legends. Can only be done after
2438 * the first frame has been read... (Then we know how many there are)
2440 blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL);
2441 if (bDisRe && bDRAll && !leg && blk_disre)
2446 fa = top->idef.il[F_DISRES].iatoms;
2447 ip = top->idef.iparams;
2448 if (blk_disre->nsub != 2 ||
2449 (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2451 gmx_incons("Number of disre sub-blocks not equal to 2");
2454 ndisre = blk_disre->sub[0].nr;
2455 if (ndisre != top->idef.il[F_DISRES].nr/3)
2457 gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2458 ndisre, top->idef.il[F_DISRES].nr/3);
2460 snew(pairleg, ndisre);
2461 for (i = 0; i < ndisre; i++)
2463 snew(pairleg[i], 30);
2466 gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j);
2467 gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k);
2468 sprintf(pairleg[i], "%d %s %d %s (%d)",
2469 resnr_j, anm_j, resnr_k, anm_k,
2470 ip[fa[3*i]].disres.label);
2472 set = select_it(ndisre, pairleg, &nset);
2474 for (i = 0; (i < nset); i++)
2477 sprintf(leg[2*i], "a %s", pairleg[set[i]]);
2478 snew(leg[2*i+1], 32);
2479 sprintf(leg[2*i+1], "i %s", pairleg[set[i]]);
2481 xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv);
2485 * Store energies for analysis afterwards...
2487 if (!bDisRe && !bDHDL && (fr->nre > 0))
2489 if (edat.nframes % 1000 == 0)
2491 srenew(time, edat.nframes+1000);
2493 time[edat.nframes] = fr->t;
2497 * Printing time, only when we do not want to skip frames
2499 if (!skip || teller % skip == 0)
2503 /*******************************************
2504 * D I S T A N C E R E S T R A I N T S
2505 *******************************************/
2509 float *disre_rt = blk_disre->sub[0].fval;
2510 float *disre_rm3tav = blk_disre->sub[1].fval;
2512 double *disre_rt = blk_disre->sub[0].dval;
2513 double *disre_rm3tav = blk_disre->sub[1].dval;
2516 print_time(out, fr->t);
2517 if (violaver == NULL)
2519 snew(violaver, ndisre);
2522 /* Subtract bounds from distances, to calculate violations */
2523 calc_violations(disre_rt, disre_rm3tav,
2524 nbounds, pair, bounds, violaver, &sumt, &sumaver);
2526 fprintf(out, " %8.4f %8.4f\n", sumaver, sumt);
2529 print_time(fp_pairs, fr->t);
2530 for (i = 0; (i < nset); i++)
2533 fprintf(fp_pairs, " %8.4f", mypow(disre_rm3tav[sss], minthird));
2534 fprintf(fp_pairs, " %8.4f", disre_rt[sss]);
2536 fprintf(fp_pairs, "\n");
2543 do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2546 /*******************************************
2548 *******************************************/
2555 /* We skip frames with single points (usually only the first frame),
2556 * since they would result in an average plot with outliers.
2560 print_time(out, fr->t);
2561 print1(out, bDp, fr->ener[set[0]].e);
2562 print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
2563 print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum));
2569 print_time(out, fr->t);
2573 for (i = 0; i < nset; i++)
2575 sum += fr->ener[set[i]].e;
2577 print1(out, bDp, sum/nmol-ezero);
2581 for (i = 0; (i < nset); i++)
2585 print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero);
2589 print1(out, bDp, fr->ener[set[i]].e);
2596 blk = find_block_id_enxframe(fr, enx_i, NULL);
2600 xdr_datatype dt = xdr_datatype_float;
2602 xdr_datatype dt = xdr_datatype_double;
2606 if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2608 gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2611 vals = blk->sub[0].fval;
2613 vals = blk->sub[0].dval;
2616 if (blk->sub[0].nr != (size_t)nor)
2618 gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2622 for (i = 0; i < nor; i++)
2624 orient[i] += vals[i];
2629 for (i = 0; i < nor; i++)
2631 odrms[i] += sqr(vals[i]-oobs[i]);
2636 fprintf(fort, " %10f", fr->t);
2637 for (i = 0; i < norsel; i++)
2639 fprintf(fort, " %g", vals[orsel[i]]);
2641 fprintf(fort, "\n");
2645 fprintf(fodt, " %10f", fr->t);
2646 for (i = 0; i < norsel; i++)
2648 fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]);
2650 fprintf(fodt, "\n");
2654 blk = find_block_id_enxframe(fr, enxORT, NULL);
2658 xdr_datatype dt = xdr_datatype_float;
2660 xdr_datatype dt = xdr_datatype_double;
2664 if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2666 gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2669 vals = blk->sub[0].fval;
2671 vals = blk->sub[0].dval;
2674 if (blk->sub[0].nr != (size_t)(nex*12))
2676 gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2677 blk->sub[0].nr/12, nex);
2679 fprintf(foten, " %10f", fr->t);
2680 for (i = 0; i < nex; i++)
2682 for (j = 0; j < (bOvec ? 12 : 3); j++)
2684 fprintf(foten, " %g", vals[i*12+j]);
2687 fprintf(foten, "\n");
2694 while (bCont && (timecheck == 0));
2696 fprintf(stderr, "\n");
2718 out = xvgropen(opt2fn("-ora", NFILE, fnm),
2719 "Average calculated orientations",
2720 "Restraint label", "", oenv);
2723 fprintf(out, "%s", orinst_sub);
2725 for (i = 0; i < nor; i++)
2727 fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr);
2733 out = xvgropen(opt2fn("-oda", NFILE, fnm),
2734 "Average restraint deviation",
2735 "Restraint label", "", oenv);
2738 fprintf(out, "%s", orinst_sub);
2740 for (i = 0; i < nor; i++)
2742 fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr-oobs[i]);
2748 out = xvgropen(opt2fn("-odr", NFILE, fnm),
2749 "RMS orientation restraint deviations",
2750 "Restraint label", "", oenv);
2753 fprintf(out, "%s", orinst_sub);
2755 for (i = 0; i < nor; i++)
2757 fprintf(out, "%5d %g\n", or_label[i], sqrt(odrms[i]/norfr));
2768 analyse_disre(opt2fn("-viol", NFILE, fnm),
2769 teller_disre, violaver, bounds, index, pair, nbounds, oenv);
2776 printf("\n\nWrote %d lambda values with %d samples as ",
2777 dh_lambdas, dh_samples);
2780 printf("%d dH histograms ", dh_hists);
2784 printf("%d dH data blocks ", dh_blocks);
2786 printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2791 gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2797 double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2798 analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
2799 bFee, bSum, opt2parg_bSet("-nmol", npargs, ppa),
2800 bVisco, opt2fn("-vis", NFILE, fnm),
2802 start_step, start_t, frame[cur].step, frame[cur].t,
2803 time, reftemp, &edat,
2804 nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
2808 calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat,
2812 if (opt2bSet("-f2", NFILE, fnm))
2814 fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm),
2815 reftemp, nset, set, leg, &edat, time, oenv);
2819 const char *nxy = "-nxy";
2821 do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2822 do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2823 do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
2824 do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
2825 do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
2826 do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
2827 do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
2828 do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
2829 do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);