1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Green Red Orange Magenta Azure Cyan Skyblue
45 #include "gmx_fatal.h"
59 #include "mtop_util.h"
63 static real minthird = -1.0/3.0, minsixth = -1.0/6.0;
81 gmx_large_int_t nsteps;
82 gmx_large_int_t npoints;
90 static double mypow(double x, double y)
102 static int *select_it(int nre, char *nm[], int *nset)
107 gmx_bool bVerbose = TRUE;
109 if ((getenv("VERBOSE")) != NULL)
114 fprintf(stderr, "Select the terms you want from the following list\n");
115 fprintf(stderr, "End your selection with 0\n");
119 for (k = 0; (k < nre); )
121 for (j = 0; (j < 4) && (k < nre); j++, k++)
123 fprintf(stderr, " %3d=%14s", k+1, nm[k]);
125 fprintf(stderr, "\n");
132 if (1 != scanf("%d", &n))
134 gmx_fatal(FARGS, "Error reading user input");
136 if ((n > 0) && (n <= nre))
144 for (i = (*nset) = 0; (i < nre); i++)
157 static void chomp(char *buf)
159 int len = strlen(buf);
161 while ((len > 0) && (buf[len-1] == '\n'))
168 static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
171 int n, k, kk, j, i, nmatch, nind, nss;
173 gmx_bool bEOF, bVerbose = TRUE, bLong = FALSE;
174 char *ptr, buf[STRLEN];
175 const char *fm4 = "%3d %-14s";
176 const char *fm2 = "%3d %-34s";
179 if ((getenv("VERBOSE")) != NULL)
184 fprintf(stderr, "\n");
185 fprintf(stderr, "Select the terms you want from the following list by\n");
186 fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
187 fprintf(stderr, "End your selection with an empty line or a zero.\n");
188 fprintf(stderr, "-------------------------------------------------------------------\n");
192 for (k = 0; k < nre; k++)
194 newnm[k] = strdup(nm[k].name);
195 /* Insert dashes in all the names */
196 while ((ptr = strchr(newnm[k], ' ')) != NULL)
206 fprintf(stderr, "\n");
209 for (kk = k; kk < k+4; kk++)
211 if (kk < nre && strlen(nm[kk].name) > 14)
219 fprintf(stderr, " ");
223 fprintf(stderr, fm4, k+1, newnm[k]);
232 fprintf(stderr, fm2, k+1, newnm[k]);
243 fprintf(stderr, "\n\n");
249 while (!bEOF && (fgets2(buf, STRLEN-1, stdin)))
251 /* Remove newlines */
257 /* Empty line means end of input */
258 bEOF = (strlen(buf) == 0);
266 /* First try to read an integer */
267 nss = sscanf(ptr, "%d", &nind);
270 /* Zero means end of input */
275 else if ((1 <= nind) && (nind <= nre))
281 fprintf(stderr, "number %d is out of range\n", nind);
286 /* Now try to read a string */
289 for (nind = 0; nind < nre; nind++)
291 if (gmx_strcasecmp(newnm[nind], ptr) == 0)
301 for (nind = 0; nind < nre; nind++)
303 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
311 fprintf(stderr, "String '%s' does not match anything\n", ptr);
316 /* Look for the first space, and remove spaces from there */
317 if ((ptr = strchr(ptr, ' ')) != NULL)
322 while (!bEOF && (ptr && (strlen(ptr) > 0)));
327 for (i = (*nset) = 0; (i < nre); i++)
339 gmx_fatal(FARGS, "No energy terms selected");
342 for (i = 0; (i < nre); i++)
351 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
358 /* all we need is the ir to be able to write the label */
359 read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, &mtop);
362 static void get_orires_parms(const char *topnm,
363 int *nor, int *nex, int **label, real **obs)
374 read_tpx(topnm, &ir, box, &natoms, NULL, NULL, NULL, &mtop);
375 top = gmx_mtop_generate_local_top(&mtop, &ir);
377 ip = top->idef.iparams;
378 iatom = top->idef.il[F_ORIRES].iatoms;
380 /* Count how many distance restraint there are... */
381 nb = top->idef.il[F_ORIRES].nr;
384 gmx_fatal(FARGS, "No orientation restraints in topology!\n");
391 for (i = 0; i < nb; i += 3)
393 (*label)[i/3] = ip[iatom[i]].orires.label;
394 (*obs)[i/3] = ip[iatom[i]].orires.obs;
395 if (ip[iatom[i]].orires.ex >= *nex)
397 *nex = ip[iatom[i]].orires.ex+1;
400 fprintf(stderr, "Found %d orientation restraints with %d experiments",
404 static int get_bounds(const char *topnm,
405 real **bounds, int **index, int **dr_pair, int *npairs,
406 gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir)
409 t_functype *functype;
411 int natoms, i, j, k, type, ftype, natom;
419 read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, mtop);
421 top = gmx_mtop_generate_local_top(mtop, ir);
424 functype = top->idef.functype;
425 ip = top->idef.iparams;
427 /* Count how many distance restraint there are... */
428 nb = top->idef.il[F_DISRES].nr;
431 gmx_fatal(FARGS, "No distance restraints in topology!\n");
434 /* Allocate memory */
439 /* Fill the bound array */
441 for (i = 0; (i < top->idef.ntypes); i++)
444 if (ftype == F_DISRES)
447 label1 = ip[i].disres.label;
448 b[nb] = ip[i].disres.up1;
455 /* Fill the index array */
457 disres = &(top->idef.il[F_DISRES]);
458 iatom = disres->iatoms;
459 for (i = j = k = 0; (i < disres->nr); )
462 ftype = top->idef.functype[type];
463 natom = interaction_function[ftype].nratoms+1;
464 if (label1 != top->idef.iparams[type].disres.label)
467 label1 = top->idef.iparams[type].disres.label;
477 gmx_incons("get_bounds for distance restraints");
486 static void calc_violations(real rt[], real rav3[], int nb, int index[],
487 real bounds[], real *viol, double *st, double *sa)
489 const real sixth = 1.0/6.0;
491 double rsum, rav, sumaver, sumt;
495 for (i = 0; (i < nb); i++)
499 for (j = index[i]; (j < index[i+1]); j++)
503 viol[j] += mypow(rt[j], -3.0);
506 rsum += mypow(rt[j], -6);
508 rsum = max(0.0, mypow(rsum, -sixth)-bounds[i]);
509 rav = max(0.0, mypow(rav, -sixth)-bounds[i]);
518 static void analyse_disre(const char *voutfn, int nframes,
519 real violaver[], real bounds[], int index[],
520 int pair[], int nbounds,
521 const output_env_t oenv)
524 double sum, sumt, sumaver;
527 /* Subtract bounds from distances, to calculate violations */
528 calc_violations(violaver, violaver,
529 nbounds, pair, bounds, NULL, &sumt, &sumaver);
532 fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
534 fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n",
537 vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm",
541 for (i = 0; (i < nbounds); i++)
543 /* Do ensemble averaging */
545 for (j = pair[i]; (j < pair[i+1]); j++)
547 sumaver += sqr(violaver[j]/nframes);
549 sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]);
552 sum = max(sum, sumaver);
553 fprintf(vout, "%10d %10.5e\n", index[i], sumaver);
556 for (j = 0; (j < dr.ndr); j++)
558 fprintf(vout, "%10d %10.5e\n", j, mypow(violaver[j]/nframes, minthird));
563 fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
565 fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
567 do_view(oenv, voutfn, "-graphtype bar");
570 static void einstein_visco(const char *fn, const char *fni, int nsets,
571 int nframes, real **sum,
572 real V, real T, int nsteps, double time[],
573 const output_env_t oenv)
576 double av[4], avold[4];
585 dt = (time[1]-time[0]);
588 for (i = 0; i <= nsets; i++)
592 fp0 = xvgropen(fni, "Shear viscosity integral",
593 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
594 fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
595 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
596 for (i = 1; i < nf4; i++)
598 fac = dt*nframes/nsteps;
599 for (m = 0; m <= nsets; m++)
603 for (j = 0; j < nframes-i; j++)
605 for (m = 0; m < nsets; m++)
607 di = sqr(fac*(sum[m][j+i]-sum[m][j]));
610 av[nsets] += di/nsets;
613 /* Convert to SI for the viscosity */
614 fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i);
615 fprintf(fp0, "%10g", time[i]-time[0]);
616 for (m = 0; (m <= nsets); m++)
619 fprintf(fp0, " %10g", av[m]);
622 fprintf(fp1, "%10g", 0.5*(time[i]+time[i-1])-time[0]);
623 for (m = 0; (m <= nsets); m++)
625 fprintf(fp1, " %10g", (av[m]-avold[m])/dt);
645 gmx_large_int_t nst_min;
648 static void clear_ee_sum(ee_sum_t *ees)
656 static void add_ee_sum(ee_sum_t *ees, double sum, int np)
662 static void add_ee_av(ee_sum_t *ees)
666 av = ees->sum/ees->np;
673 static double calc_ee2(int nb, ee_sum_t *ees)
675 return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
678 static void set_ee_av(ener_ee_t *eee)
682 char buf[STEPSTRSIZE];
683 fprintf(debug, "Storing average for err.est.: %s steps\n",
684 gmx_step_str(eee->nst, buf));
686 add_ee_av(&eee->sum);
688 if (eee->b == 1 || eee->nst < eee->nst_min)
690 eee->nst_min = eee->nst;
695 static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
698 double sum, sum2, sump, see2;
699 gmx_large_int_t steps, np, p, bound_nb;
703 double x, sx, sy, sxx, sxy;
706 /* Check if we have exact statistics over all points */
707 for (i = 0; i < nset; i++)
710 ed->bExactStat = FALSE;
711 if (edat->npoints > 0)
713 /* All energy file sum entries 0 signals no exact sums.
714 * But if all energy values are 0, we still have exact sums.
717 for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
719 if (ed->ener[i] != 0)
723 ed->bExactStat = (ed->es[f].sum != 0);
727 ed->bExactStat = TRUE;
733 for (i = 0; i < nset; i++)
744 for (nb = nbmin; nb <= nbmax; nb++)
747 clear_ee_sum(&eee[nb].sum);
751 for (f = 0; f < edat->nframes; f++)
757 /* Add the sum and the sum of variances to the totals. */
763 sum2 += dsqr(sum/np - (sum + es->sum)/(np + p))
769 /* Add a single value to the sum and sum of squares. */
775 /* sum has to be increased after sum2 */
779 /* For the linear regression use variance 1/p.
780 * Note that sump is the sum, not the average, so we don't need p*.
782 x = edat->step[f] - 0.5*(edat->steps[f] - 1);
788 for (nb = nbmin; nb <= nbmax; nb++)
790 /* Check if the current end step is closer to the desired
791 * block boundary than the next end step.
793 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
794 if (eee[nb].nst > 0 &&
795 bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
805 eee[nb].nst += edat->step[f] - edat->step[f-1];
809 add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
813 add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
815 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
816 if (edat->step[f]*nb >= bound_nb)
823 edat->s[i].av = sum/np;
826 edat->s[i].rmsd = sqrt(sum2/np);
830 edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av));
833 if (edat->nframes > 1)
835 edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
839 edat->s[i].slope = 0;
844 for (nb = nbmin; nb <= nbmax; nb++)
846 /* Check if we actually got nb blocks and if the smallest
847 * block is not shorter than 80% of the average.
851 char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
852 fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
854 gmx_step_str(eee[nb].nst_min, buf1),
855 gmx_step_str(edat->nsteps, buf2));
857 if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
859 see2 += calc_ee2(nb, &eee[nb].sum);
865 edat->s[i].ee = sqrt(see2/nee);
875 static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax)
886 snew(s->ener, esum->nframes);
887 snew(s->es, esum->nframes);
889 s->bExactStat = TRUE;
891 for (i = 0; i < nset; i++)
893 if (!edat->s[i].bExactStat)
895 s->bExactStat = FALSE;
897 s->slope += edat->s[i].slope;
900 for (f = 0; f < edat->nframes; f++)
903 for (i = 0; i < nset; i++)
905 sum += edat->s[i].ener[f];
909 for (i = 0; i < nset; i++)
911 sum += edat->s[i].es[f].sum;
917 calc_averages(1, esum, nbmin, nbmax);
922 static char *ee_pr(double ee, char *buf)
929 sprintf(buf, "%s", "--");
933 /* Round to two decimals by printing. */
934 sprintf(tmp, "%.1e", ee);
935 sscanf(tmp, "%lf", &rnd);
936 sprintf(buf, "%g", rnd);
942 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat)
944 /* Remove the drift by performing a fit to y = ax+b.
945 Uses 5 iterations. */
947 double delta, d, sd, sd2;
949 edat->npoints = edat->nframes;
950 edat->nsteps = edat->nframes;
952 for (k = 0; (k < 5); k++)
954 for (i = 0; (i < nset); i++)
956 delta = edat->s[i].slope*dt;
960 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
963 for (j = 0; (j < edat->nframes); j++)
965 edat->s[i].ener[j] -= j*delta;
966 edat->s[i].es[j].sum = 0;
967 edat->s[i].es[j].sum2 = 0;
970 calc_averages(nset, edat, nbmin, nbmax);
974 static void calc_fluctuation_props(FILE *fp,
975 gmx_bool bDriftCorr, real dt,
977 char **leg, enerdata_t *edat,
978 int nbmin, int nbmax)
981 double vvhh, vv, v, h, hh2, vv2, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
984 eVol, eEnth, eTemp, eEtot, eNR
986 const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
989 NANO3 = NANO*NANO*NANO;
992 fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n"
993 "for spurious drift in the graphs. Note that this is not\n"
994 "a substitute for proper equilibration and sampling!\n");
998 remove_drift(nset, nbmin, nbmax, dt, edat);
1000 for (i = 0; (i < eNR); i++)
1002 for (ii[i] = 0; (ii[i] < nset &&
1003 (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++)
1007 /* if (ii[i] < nset)
1008 fprintf(fp,"Found %s data.\n",my_ener[i]);
1010 /* Compute it all! */
1011 vvhh = alpha = kappa = cp = dcp = cv = NOTSET;
1015 if (ii[eTemp] < nset)
1017 tt = edat->s[ii[eTemp]].av;
1021 if ((ii[eVol] < nset) && (ii[eTemp] < nset))
1023 vv = edat->s[ii[eVol]].av*NANO3;
1024 varv = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
1025 kappa = (varv/vv)/(BOLTZMANN*tt);
1029 if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
1031 hh = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
1032 varh = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
1033 cp = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
1036 et = varet = NOTSET;
1037 if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
1039 /* Only compute cv in constant volume runs, which we can test
1040 by checking whether the enthalpy was computed.
1042 et = edat->s[ii[eEtot]].av;
1043 varet = sqr(edat->s[ii[eEtot]].rmsd);
1044 cv = KILO*((varet/nmol)/(BOLTZ*tt*tt));
1047 if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
1050 for (j = 0; (j < edat->nframes); j++)
1052 v = edat->s[ii[eVol]].ener[j]*NANO3;
1053 h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
1056 vvhh /= edat->nframes;
1057 alpha = (vvhh-vv*hh)/(vv*BOLTZMANN*tt*tt);
1058 dcp = (vv*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
1065 fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n",
1068 fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
1069 fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
1070 fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
1071 fprintf(fp, "please use the g_dos program.\n\n");
1072 fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n"
1073 "a block-averaging error analysis (not implemented in g_energy yet)\n");
1079 fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol);
1083 fprintf(fp, "vvhh = %10g (m^3 J)\n", vvhh);
1088 fprintf(fp, "Volume = %10g m^3/mol\n",
1093 fprintf(fp, "Enthalpy = %10g kJ/mol\n",
1094 hh*AVOGADRO/(KILO*nmol));
1096 if (alpha != NOTSET)
1098 fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
1101 if (kappa != NOTSET)
1103 fprintf(fp, "Isothermal Compressibility Kappa = %10g (J/m^3)\n",
1105 fprintf(fp, "Adiabatic bulk modulus = %10g (m^3/J)\n",
1110 fprintf(fp, "Heat capacity at constant pressure Cp = %10g J/mol K\n",
1115 fprintf(fp, "Heat capacity at constant volume Cv = %10g J/mol K\n",
1120 fprintf(fp, "Cp-Cv = %10g J/mol K\n",
1123 please_cite(fp, "Allen1987a");
1127 fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
1131 static void analyse_ener(gmx_bool bCorr, const char *corrfn,
1132 gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct,
1133 gmx_bool bVisco, const char *visfn, int nmol,
1134 gmx_large_int_t start_step, double start_t,
1135 gmx_large_int_t step, double t,
1136 double time[], real reftemp,
1138 int nset, int set[], gmx_bool *bIsEner,
1139 char **leg, gmx_enxnm_t *enm,
1140 real Vaver, real ezero,
1141 int nbmin, int nbmax,
1142 const output_env_t oenv)
1145 /* Check out the printed manual for equations! */
1146 double Dt, aver, stddev, errest, delta_t, totaldrift;
1147 enerdata_t *esum = NULL;
1148 real xxx, integral, intBulk, Temp = 0, Pres = 0;
1149 real sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest;
1150 double beta = 0, expE, expEtot, *fee = NULL;
1151 gmx_large_int_t nsteps;
1152 int nexact, nnotexact;
1156 char buf[256], eebuf[100];
1158 nsteps = step - start_step + 1;
1161 fprintf(stdout, "Not enough steps (%s) for statistics\n",
1162 gmx_step_str(nsteps, buf));
1166 /* Calculate the time difference */
1167 delta_t = t - start_t;
1169 fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
1170 gmx_step_str(nsteps, buf), start_t, t, nset);
1172 calc_averages(nset, edat, nbmin, nbmax);
1176 esum = calc_sum(nset, edat, nbmin, nbmax);
1179 if (edat->npoints == 0)
1188 for (i = 0; (i < nset); i++)
1190 if (edat->s[i].bExactStat)
1203 fprintf(stdout, "All statistics are over %s points\n",
1204 gmx_step_str(edat->npoints, buf));
1206 else if (nexact == 0 || edat->npoints == edat->nframes)
1208 fprintf(stdout, "All statistics are over %d points (frames)\n",
1213 fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
1214 for (i = 0; (i < nset); i++)
1216 if (!edat->s[i].bExactStat)
1218 fprintf(stdout, " '%s'", leg[i]);
1221 fprintf(stdout, " %s has statistics over %d points (frames)\n",
1222 nnotexact == 1 ? "is" : "are", edat->nframes);
1223 fprintf(stdout, "All other statistics are over %s points\n",
1224 gmx_step_str(edat->npoints, buf));
1226 fprintf(stdout, "\n");
1228 fprintf(stdout, "%-24s %10s %10s %10s %10s",
1229 "Energy", "Average", "Err.Est.", "RMSD", "Tot-Drift");
1232 fprintf(stdout, " %10s\n", "-kT ln<e^(E/kT)>");
1236 fprintf(stdout, "\n");
1238 fprintf(stdout, "-------------------------------------------------------------------------------\n");
1240 /* Initiate locals, only used with -sum */
1244 beta = 1.0/(BOLTZ*reftemp);
1247 for (i = 0; (i < nset); i++)
1249 aver = edat->s[i].av;
1250 stddev = edat->s[i].rmsd;
1251 errest = edat->s[i].ee;
1256 for (j = 0; (j < edat->nframes); j++)
1258 expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
1262 expEtot += expE/edat->nframes;
1265 fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
1267 if (strstr(leg[i], "empera") != NULL)
1271 else if (strstr(leg[i], "olum") != NULL)
1275 else if (strstr(leg[i], "essure") != NULL)
1281 pr_aver = aver/nmol-ezero;
1282 pr_stddev = stddev/nmol;
1283 pr_errest = errest/nmol;
1292 /* Multiply the slope in steps with the number of steps taken */
1293 totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
1299 fprintf(stdout, "%-24s %10g %10s %10g %10g",
1300 leg[i], pr_aver, ee_pr(pr_errest, eebuf), pr_stddev, totaldrift);
1303 fprintf(stdout, " %10g", fee[i]);
1306 fprintf(stdout, " (%s)\n", enm[set[i]].unit);
1310 for (j = 0; (j < edat->nframes); j++)
1312 edat->s[i].ener[j] -= aver;
1318 totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1319 fprintf(stdout, "%-24s %10g %10s %10s %10g (%s)",
1320 "Total", esum->s[0].av/nmol, ee_pr(esum->s[0].ee/nmol, eebuf),
1321 "--", totaldrift/nmol, enm[set[0]].unit);
1322 /* pr_aver,pr_stddev,a,totaldrift */
1325 fprintf(stdout, " %10g %10g\n",
1326 log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta);
1330 fprintf(stdout, "\n");
1334 /* Do correlation function */
1335 if (edat->nframes > 1)
1337 Dt = delta_t/(edat->nframes - 1);
1345 const char* leg[] = { "Shear", "Bulk" };
1350 /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1352 /* Symmetrise tensor! (and store in first three elements)
1353 * And subtract average pressure!
1356 for (i = 0; i < 12; i++)
1358 snew(eneset[i], edat->nframes);
1361 for (i = 0; i < 3; i++)
1363 snew(enesum[i], edat->nframes);
1365 for (i = 0; (i < edat->nframes); i++)
1367 eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1368 eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1369 eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1370 for (j = 3; j <= 11; j++)
1372 eneset[j][i] = edat->s[j].ener[i];
1374 eneset[11][i] -= Pres;
1375 enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum);
1376 enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
1377 enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
1380 einstein_visco("evisco.xvg", "eviscoi.xvg",
1381 3, edat->nframes, enesum, Vaver, Temp, nsteps, time, oenv);
1383 /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1384 /* Do it for shear viscosity */
1385 strcpy(buf, "Shear Viscosity");
1386 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3,
1387 (edat->nframes+1)/2, eneset, Dt,
1388 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1390 /* Now for bulk viscosity */
1391 strcpy(buf, "Bulk Viscosity");
1392 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1,
1393 (edat->nframes+1)/2, &(eneset[11]), Dt,
1394 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1396 factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1397 fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1398 xvgr_legend(fp, asize(leg), leg, oenv);
1400 /* Use trapezium rule for integration */
1403 nout = get_acfnout();
1404 if ((nout < 2) || (nout >= edat->nframes/2))
1406 nout = edat->nframes/2;
1408 for (i = 1; (i < nout); i++)
1410 integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor;
1411 intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1412 fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk);
1420 strcpy(buf, "Autocorrelation of Energy Fluctuations");
1424 strcpy(buf, "Energy Autocorrelation");
1427 do_autocorr(corrfn, oenv, buf, edat->nframes,
1429 bSum ? &edat->s[nset-1].ener : eneset,
1430 (delta_t/edat->nframes), eacNormal, FALSE);
1436 static void print_time(FILE *fp, double t)
1438 fprintf(fp, "%12.6f", t);
1441 static void print1(FILE *fp, gmx_bool bDp, real e)
1445 fprintf(fp, " %16.12f", e);
1449 fprintf(fp, " %10.6f", e);
1453 static void fec(const char *ene2fn, const char *runavgfn,
1454 real reftemp, int nset, int set[], char *leg[],
1455 enerdata_t *edat, double time[],
1456 const output_env_t oenv)
1458 const char * ravgleg[] = {
1459 "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1460 "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N"
1464 int nre, timecheck, step, nenergy, nenergy2, maxenergy;
1470 gmx_enxnm_t *enm = NULL;
1474 /* read second energy file */
1477 enx = open_enx(ene2fn, "r");
1478 do_enxnms(enx, &(fr->nre), &enm);
1480 snew(eneset2, nset+1);
1486 /* This loop searches for the first frame (when -b option is given),
1487 * or when this has been found it reads just one energy frame
1491 bCont = do_enx(enx, fr);
1495 timecheck = check_times(fr->t);
1499 while (bCont && (timecheck < 0));
1501 /* Store energies for analysis afterwards... */
1502 if ((timecheck == 0) && bCont)
1506 if (nenergy2 >= maxenergy)
1509 for (i = 0; i <= nset; i++)
1511 srenew(eneset2[i], maxenergy);
1514 if (fr->t != time[nenergy2])
1516 fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n",
1517 fr->t, time[nenergy2], gmx_step_str(fr->step, buf));
1519 for (i = 0; i < nset; i++)
1521 eneset2[i][nenergy2] = fr->ener[set[i]].e;
1527 while (bCont && (timecheck == 0));
1530 if (edat->nframes != nenergy2)
1532 fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n",
1533 edat->nframes, nenergy2);
1535 nenergy = min(edat->nframes, nenergy2);
1537 /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1541 fp = xvgropen(runavgfn, "Running average free energy difference",
1542 "Time (" unit_time ")", "\\8D\\4E (" unit_energy ")", oenv);
1543 xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
1545 fprintf(stdout, "\n%-24s %10s\n",
1546 "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
1548 beta = 1.0/(BOLTZ*reftemp);
1549 for (i = 0; i < nset; i++)
1551 if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
1553 fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n",
1554 leg[i], enm[set[i]].name);
1556 for (j = 0; j < nenergy; j++)
1558 dE = eneset2[i][j] - edat->s[i].ener[j];
1559 sum += exp(-dE*beta);
1562 fprintf(fp, "%10g %10g %10g\n",
1563 time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
1566 aver = -BOLTZ*reftemp*log(sum/nenergy);
1567 fprintf(stdout, "%-24s %10g\n", leg[i], aver);
1577 static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl,
1578 const char *filename, gmx_bool bDp,
1579 int *blocks, int *hists, int *samples, int *nlambdas,
1580 const output_env_t oenv)
1582 const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1583 char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
1585 gmx_bool first = FALSE;
1586 int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
1589 double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0;
1590 static int setnr = 0;
1591 double *native_lambda_vec = NULL;
1592 const char **lambda_components = NULL;
1593 int n_lambda_vec = 0;
1594 gmx_bool changing_lambda = FALSE;
1595 int lambda_fep_state;
1597 /* now count the blocks & handle the global dh data */
1598 for (i = 0; i < fr->nblock; i++)
1600 if (fr->block[i].id == enxDHHIST)
1604 else if (fr->block[i].id == enxDH)
1608 else if (fr->block[i].id == enxDHCOLL)
1611 if ( (fr->block[i].nsub < 1) ||
1612 (fr->block[i].sub[0].type != xdr_datatype_double) ||
1613 (fr->block[i].sub[0].nr < 5))
1615 gmx_fatal(FARGS, "Unexpected block data");
1618 /* read the data from the DHCOLL block */
1619 temp = fr->block[i].sub[0].dval[0];
1620 start_time = fr->block[i].sub[0].dval[1];
1621 delta_time = fr->block[i].sub[0].dval[2];
1622 start_lambda = fr->block[i].sub[0].dval[3];
1623 delta_lambda = fr->block[i].sub[0].dval[4];
1624 changing_lambda = (delta_lambda != 0);
1625 if (fr->block[i].nsub > 1)
1627 lambda_fep_state = fr->block[i].sub[1].ival[0];
1628 if (n_lambda_vec == 0)
1630 n_lambda_vec = fr->block[i].sub[1].ival[1];
1634 if (n_lambda_vec != fr->block[i].sub[1].ival[1])
1637 "Unexpected change of basis set in lambda");
1640 if (lambda_components == NULL)
1642 snew(lambda_components, n_lambda_vec);
1644 if (native_lambda_vec == NULL)
1646 snew(native_lambda_vec, n_lambda_vec);
1648 for (j = 0; j < n_lambda_vec; j++)
1650 native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j];
1651 lambda_components[j] =
1652 efpt_singular_names[fr->block[i].sub[1].ival[2+j]];
1658 if (nblock_hist == 0 && nblock_dh == 0)
1660 /* don't do anything */
1663 if (nblock_hist > 0 && nblock_dh > 0)
1665 gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1671 /* we have standard, non-histogram data --
1672 call open_dhdl to open the file */
1673 /* TODO this is an ugly hack that needs to be fixed: this will only
1674 work if the order of data is always the same and if we're
1675 only using the g_energy compiled with the mdrun that produced
1677 *fp_dhdl = open_dhdl(filename, ir, oenv);
1681 sprintf(title, "N(%s)", deltag);
1682 sprintf(label_x, "%s (%s)", deltag, unit_energy);
1683 sprintf(label_y, "Samples");
1684 *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1685 sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
1686 xvgr_subtitle(*fp_dhdl, buf, oenv);
1690 (*hists) += nblock_hist;
1691 (*blocks) += nblock_dh;
1692 (*nlambdas) = nblock_hist+nblock_dh;
1694 /* write the data */
1695 if (nblock_hist > 0)
1697 gmx_large_int_t sum = 0;
1699 for (i = 0; i < fr->nblock; i++)
1701 t_enxblock *blk = &(fr->block[i]);
1702 if (blk->id == enxDHHIST)
1704 double foreign_lambda, dx;
1706 int nhist, derivative;
1708 /* check the block types etc. */
1709 if ( (blk->nsub < 2) ||
1710 (blk->sub[0].type != xdr_datatype_double) ||
1711 (blk->sub[1].type != xdr_datatype_large_int) ||
1712 (blk->sub[0].nr < 2) ||
1713 (blk->sub[1].nr < 2) )
1715 gmx_fatal(FARGS, "Unexpected block data in file");
1717 foreign_lambda = blk->sub[0].dval[0];
1718 dx = blk->sub[0].dval[1];
1719 nhist = blk->sub[1].lval[0];
1720 derivative = blk->sub[1].lval[1];
1721 for (j = 0; j < nhist; j++)
1724 x0 = blk->sub[1].lval[2+j];
1728 sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1729 deltag, lambda, foreign_lambda,
1730 lambda, start_lambda);
1734 sprintf(legend, "N(%s | %s=%g)",
1735 dhdl, lambda, start_lambda);
1739 xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1741 for (k = 0; k < blk->sub[j+2].nr; k++)
1746 hist = blk->sub[j+2].ival[k];
1749 fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist,
1753 /* multiple histogram data blocks in one histogram
1754 mean that the second one is the reverse of the first one:
1755 for dhdl derivatives, it's important to know both the
1756 maximum and minimum values */
1761 (*samples) += (int)(sum/nblock_hist);
1767 char **setnames = NULL;
1768 int nnames = nblock_dh;
1770 for (i = 0; i < fr->nblock; i++)
1772 t_enxblock *blk = &(fr->block[i]);
1773 if (blk->id == enxDH)
1777 len = blk->sub[2].nr;
1781 if (len != blk->sub[2].nr)
1783 gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1790 for (i = 0; i < len; i++)
1792 double time = start_time + delta_time*i;
1794 fprintf(*fp_dhdl, "%.4f ", time);
1796 for (j = 0; j < fr->nblock; j++)
1798 t_enxblock *blk = &(fr->block[j]);
1799 if (blk->id == enxDH)
1802 if (blk->sub[2].type == xdr_datatype_float)
1804 value = blk->sub[2].fval[i];
1808 value = blk->sub[2].dval[i];
1810 /* we need to decide which data type it is based on the count*/
1812 if (j == 1 && ir->bExpanded)
1814 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! */
1820 fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */
1824 fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */
1829 fprintf(*fp_dhdl, "\n");
1835 int gmx_energy(int argc, char *argv[])
1837 const char *desc[] = {
1839 "[TT]g_energy[tt] extracts energy components or distance restraint",
1840 "data from an energy file. The user is prompted to interactively",
1841 "select the desired energy terms.[PAR]",
1843 "Average, RMSD, and drift are calculated with full precision from the",
1844 "simulation (see printed manual). Drift is calculated by performing",
1845 "a least-squares fit of the data to a straight line. The reported total drift",
1846 "is the difference of the fit at the first and last point.",
1847 "An error estimate of the average is given based on a block averages",
1848 "over 5 blocks using the full-precision averages. The error estimate",
1849 "can be performed over multiple block lengths with the options",
1850 "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1851 "[BB]Note[bb] that in most cases the energy files contains averages over all",
1852 "MD steps, or over many more points than the number of frames in",
1853 "energy file. This makes the [TT]g_energy[tt] statistics output more accurate",
1854 "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1855 "file, the statistics mentioned above are simply over the single, per-frame",
1856 "energy values.[PAR]",
1858 "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1860 "Some fluctuation-dependent properties can be calculated provided",
1861 "the correct energy terms are selected, and that the command line option",
1862 "[TT]-fluct_props[tt] is given. The following properties",
1863 "will be computed:[BR]",
1864 "Property Energy terms needed[BR]",
1865 "---------------------------------------------------[BR]",
1866 "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp [BR]",
1867 "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp [BR]",
1868 "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1869 "Isothermal compressibility: Vol, Temp [BR]",
1870 "Adiabatic bulk modulus: Vol, Temp [BR]",
1871 "---------------------------------------------------[BR]",
1872 "You always need to set the number of molecules [TT]-nmol[tt].",
1873 "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1874 "for quantum effects. Use the [TT]g_dos[tt] program if you need that (and you do).[PAR]"
1875 "When the [TT]-viol[tt] option is set, the time averaged",
1876 "violations are plotted and the running time-averaged and",
1877 "instantaneous sum of violations are recalculated. Additionally",
1878 "running time-averaged and instantaneous distances between",
1879 "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1881 "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1882 "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1883 "The first two options plot the orientation, the last three the",
1884 "deviations of the orientations from the experimental values.",
1885 "The options that end on an 'a' plot the average over time",
1886 "as a function of restraint. The options that end on a 't'",
1887 "prompt the user for restraint label numbers and plot the data",
1888 "as a function of time. Option [TT]-odr[tt] plots the RMS",
1889 "deviation as a function of restraint.",
1890 "When the run used time or ensemble averaged orientation restraints,",
1891 "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1892 "not ensemble-averaged orientations and deviations instead of",
1893 "the time and ensemble averages.[PAR]",
1895 "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1896 "tensor for each orientation restraint experiment. With option",
1897 "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1899 "Option [TT]-odh[tt] extracts and plots the free energy data",
1900 "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1901 "from the [TT]ener.edr[tt] file.[PAR]",
1903 "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1904 "difference with an ideal gas state: [BR]",
1905 " [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]",
1906 " [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]",
1907 "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1908 "the average is over the ensemble (or time in a trajectory).",
1909 "Note that this is in principle",
1910 "only correct when averaging over the whole (Boltzmann) ensemble",
1911 "and using the potential energy. This also allows for an entropy",
1912 "estimate using:[BR]",
1913 " [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]",
1914 " [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",
1917 "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1918 "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1919 "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1920 "files, and the average is over the ensemble A. The running average",
1921 "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1922 "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1925 static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1926 static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1927 static int skip = 0, nmol = 1, nbmin = 5, nbmax = 5;
1928 static real reftemp = 300.0, ezero = 0;
1930 { "-fee", FALSE, etBOOL, {&bFee},
1931 "Do a free energy estimate" },
1932 { "-fetemp", FALSE, etREAL, {&reftemp},
1933 "Reference temperature for free energy calculation" },
1934 { "-zero", FALSE, etREAL, {&ezero},
1935 "Subtract a zero-point energy" },
1936 { "-sum", FALSE, etBOOL, {&bSum},
1937 "Sum the energy terms selected rather than display them all" },
1938 { "-dp", FALSE, etBOOL, {&bDp},
1939 "Print energies in high precision" },
1940 { "-nbmin", FALSE, etINT, {&nbmin},
1941 "Minimum number of blocks for error estimate" },
1942 { "-nbmax", FALSE, etINT, {&nbmax},
1943 "Maximum number of blocks for error estimate" },
1944 { "-mutot", FALSE, etBOOL, {&bMutot},
1945 "Compute the total dipole moment from the components" },
1946 { "-skip", FALSE, etINT, {&skip},
1947 "Skip number of frames between data points" },
1948 { "-aver", FALSE, etBOOL, {&bPrAll},
1949 "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1950 { "-nmol", FALSE, etINT, {&nmol},
1951 "Number of molecules in your sample: the energies are divided by this number" },
1952 { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1953 "Compute properties based on energy fluctuations, like heat capacity" },
1954 { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1955 "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1956 { "-fluc", FALSE, etBOOL, {&bFluct},
1957 "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1958 { "-orinst", FALSE, etBOOL, {&bOrinst},
1959 "Analyse instantaneous orientation data" },
1960 { "-ovec", FALSE, etBOOL, {&bOvec},
1961 "Also plot the eigenvectors with [TT]-oten[tt]" }
1963 const char * drleg[] = {
1967 static const char *setnm[] = {
1968 "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1969 "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1970 "Volume", "Pressure"
1973 FILE *out = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL;
1974 FILE *fp_dhdl = NULL;
1979 gmx_localtop_t *top = NULL;
1983 gmx_enxnm_t *enm = NULL;
1984 t_enxframe *frame, *fr = NULL;
1986 #define NEXT (1-cur)
1987 int nre, teller, teller_disre, nfr;
1988 gmx_large_int_t start_step;
1989 int nor = 0, nex = 0, norfr = 0, enx_i = 0;
1991 real *bounds = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL;
1992 int *index = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL;
1993 int nbounds = 0, npairs;
1994 gmx_bool bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL;
1995 gmx_bool bFoundStart, bCont, bEDR, bVisco;
1996 double sum, sumaver, sumt, ener, dbl;
1997 double *time = NULL;
1999 int *set = NULL, i, j, k, nset, sss;
2000 gmx_bool *bIsEner = NULL;
2001 char **pairleg, **odtleg, **otenleg;
2004 char *anm_j, *anm_k, *resnm_j, *resnm_k;
2005 int resnr_j, resnr_k;
2006 const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
2009 t_enxblock *blk = NULL;
2010 t_enxblock *blk_disre = NULL;
2012 int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
2015 { efEDR, "-f", NULL, ffREAD },
2016 { efEDR, "-f2", NULL, ffOPTRD },
2017 { efTPX, "-s", NULL, ffOPTRD },
2018 { efXVG, "-o", "energy", ffWRITE },
2019 { efXVG, "-viol", "violaver", ffOPTWR },
2020 { efXVG, "-pairs", "pairs", ffOPTWR },
2021 { efXVG, "-ora", "orienta", ffOPTWR },
2022 { efXVG, "-ort", "orientt", ffOPTWR },
2023 { efXVG, "-oda", "orideva", ffOPTWR },
2024 { efXVG, "-odr", "oridevr", ffOPTWR },
2025 { efXVG, "-odt", "oridevt", ffOPTWR },
2026 { efXVG, "-oten", "oriten", ffOPTWR },
2027 { efXVG, "-corr", "enecorr", ffOPTWR },
2028 { efXVG, "-vis", "visco", ffOPTWR },
2029 { efXVG, "-ravg", "runavgdf", ffOPTWR },
2030 { efXVG, "-odh", "dhdl", ffOPTWR }
2032 #define NFILE asize(fnm)
2037 ppa = add_acf_pargs(&npargs, pa);
2038 if (!parse_common_args(&argc, argv,
2039 PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
2040 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
2045 bDRAll = opt2bSet("-pairs", NFILE, fnm);
2046 bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
2047 bORA = opt2bSet("-ora", NFILE, fnm);
2048 bORT = opt2bSet("-ort", NFILE, fnm);
2049 bODA = opt2bSet("-oda", NFILE, fnm);
2050 bODR = opt2bSet("-odr", NFILE, fnm);
2051 bODT = opt2bSet("-odt", NFILE, fnm);
2052 bORIRE = bORA || bORT || bODA || bODR || bODT;
2053 bOTEN = opt2bSet("-oten", NFILE, fnm);
2054 bDHDL = opt2bSet("-odh", NFILE, fnm);
2059 fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
2060 do_enxnms(fp, &nre, &enm);
2064 bVisco = opt2bSet("-vis", NFILE, fnm);
2066 if ((!bDisRe) && (!bDHDL))
2070 nset = asize(setnm);
2072 /* This is nasty code... To extract Pres tensor, Volume and Temperature */
2073 for (j = 0; j < nset; j++)
2075 for (i = 0; i < nre; i++)
2077 if (strstr(enm[i].name, setnm[j]))
2085 if (gmx_strcasecmp(setnm[j], "Volume") == 0)
2087 printf("Enter the box volume (" unit_volume "): ");
2088 if (1 != scanf("%lf", &dbl))
2090 gmx_fatal(FARGS, "Error reading user input");
2096 gmx_fatal(FARGS, "Could not find term %s for viscosity calculation",
2104 set = select_by_name(nre, enm, &nset);
2106 /* Print all the different units once */
2107 sprintf(buf, "(%s)", enm[set[0]].unit);
2108 for (i = 1; i < nset; i++)
2110 for (j = 0; j < i; j++)
2112 if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
2120 strcat(buf, enm[set[i]].unit);
2124 out = xvgropen(opt2fn("-o", NFILE, fnm), "Gromacs Energies", "Time (ps)", buf,
2128 for (i = 0; (i < nset); i++)
2130 leg[i] = enm[set[i]].name;
2134 leg[nset] = strdup("Sum");
2135 xvgr_legend(out, nset+1, (const char**)leg, oenv);
2139 xvgr_legend(out, nset, (const char**)leg, oenv);
2142 snew(bIsEner, nset);
2143 for (i = 0; (i < nset); i++)
2146 for (j = 0; (j <= F_ETOT); j++)
2148 bIsEner[i] = bIsEner[i] ||
2149 (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
2153 if (bPrAll && nset > 1)
2155 gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
2160 if (bORIRE || bOTEN)
2162 get_orires_parms(ftp2fn(efTPX, NFILE, fnm), &nor, &nex, &or_label, &oobs);
2186 fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n");
2187 fprintf(stderr, "End your selection with 0\n");
2194 if (1 != scanf("%d", &(orsel[j])))
2196 gmx_fatal(FARGS, "Error reading user input");
2199 while (orsel[j] > 0);
2202 fprintf(stderr, "Selecting all %d orientation restraints\n", nor);
2205 for (i = 0; i < nor; i++)
2212 /* Build the selection */
2214 for (i = 0; i < j; i++)
2216 for (k = 0; k < nor; k++)
2218 if (or_label[k] == orsel[i])
2227 fprintf(stderr, "Orientation restraint label %d not found\n",
2232 snew(odtleg, norsel);
2233 for (i = 0; i < norsel; i++)
2235 snew(odtleg[i], 256);
2236 sprintf(odtleg[i], "%d", or_label[orsel[i]]);
2240 fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
2241 "Time (ps)", "", oenv);
2244 fprintf(fort, "%s", orinst_sub);
2246 xvgr_legend(fort, norsel, (const char**)odtleg, oenv);
2250 fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
2251 "Orientation restraint deviation",
2252 "Time (ps)", "", oenv);
2255 fprintf(fodt, "%s", orinst_sub);
2257 xvgr_legend(fodt, norsel, (const char**)odtleg, oenv);
2263 foten = xvgropen(opt2fn("-oten", NFILE, fnm),
2264 "Order tensor", "Time (ps)", "", oenv);
2265 snew(otenleg, bOvec ? nex*12 : nex*3);
2266 for (i = 0; i < nex; i++)
2268 for (j = 0; j < 3; j++)
2270 sprintf(buf, "eig%d", j+1);
2271 otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
2275 for (j = 0; j < 9; j++)
2277 sprintf(buf, "vec%d%s", j/3+1, j%3 == 0 ? "x" : (j%3 == 1 ? "y" : "z"));
2278 otenleg[12*i+3+j] = strdup(buf);
2282 xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv);
2287 nbounds = get_bounds(ftp2fn(efTPX, NFILE, fnm), &bounds, &index, &pair, &npairs,
2289 snew(violaver, npairs);
2290 out = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations",
2291 "Time (ps)", "nm", oenv);
2292 xvgr_legend(out, 2, drleg, oenv);
2295 fp_pairs = xvgropen(opt2fn("-pairs", NFILE, fnm), "Pair Distances",
2296 "Time (ps)", "Distance (nm)", oenv);
2297 if (output_env_get_print_xvgr_codes(oenv))
2299 fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
2306 get_dhdl_parms(ftp2fn(efTPX, NFILE, fnm), &ir);
2309 /* Initiate energies and set them to zero */
2318 /* Initiate counters */
2321 bFoundStart = FALSE;
2326 /* This loop searches for the first frame (when -b option is given),
2327 * or when this has been found it reads just one energy frame
2331 bCont = do_enx(fp, &(frame[NEXT]));
2334 timecheck = check_times(frame[NEXT].t);
2337 while (bCont && (timecheck < 0));
2339 if ((timecheck == 0) && bCont)
2341 /* We read a valid frame, so we can use it */
2342 fr = &(frame[NEXT]);
2346 /* The frame contains energies, so update cur */
2349 if (edat.nframes % 1000 == 0)
2351 srenew(edat.step, edat.nframes+1000);
2352 memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
2353 srenew(edat.steps, edat.nframes+1000);
2354 memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
2355 srenew(edat.points, edat.nframes+1000);
2356 memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
2358 for (i = 0; i < nset; i++)
2360 srenew(edat.s[i].ener, edat.nframes+1000);
2361 memset(&(edat.s[i].ener[edat.nframes]), 0,
2362 1000*sizeof(edat.s[i].ener[0]));
2363 srenew(edat.s[i].es, edat.nframes+1000);
2364 memset(&(edat.s[i].es[edat.nframes]), 0,
2365 1000*sizeof(edat.s[i].es[0]));
2370 edat.step[nfr] = fr->step;
2375 /* Initiate the previous step data */
2376 start_step = fr->step;
2378 /* Initiate the energy sums */
2379 edat.steps[nfr] = 1;
2380 edat.points[nfr] = 1;
2381 for (i = 0; i < nset; i++)
2384 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2385 edat.s[i].es[nfr].sum2 = 0;
2392 edat.steps[nfr] = fr->nsteps;
2394 if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2398 edat.points[nfr] = 1;
2399 for (i = 0; i < nset; i++)
2402 edat.s[i].es[nfr].sum = fr->ener[sss].e;
2403 edat.s[i].es[nfr].sum2 = 0;
2409 edat.points[nfr] = fr->nsum;
2410 for (i = 0; i < nset; i++)
2413 edat.s[i].es[nfr].sum = fr->ener[sss].esum;
2414 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2416 edat.npoints += fr->nsum;
2421 /* The interval does not match fr->nsteps:
2422 * can not do exact averages.
2426 edat.nsteps = fr->step - start_step + 1;
2429 for (i = 0; i < nset; i++)
2431 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2435 * Define distance restraint legends. Can only be done after
2436 * the first frame has been read... (Then we know how many there are)
2438 blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL);
2439 if (bDisRe && bDRAll && !leg && blk_disre)
2444 fa = top->idef.il[F_DISRES].iatoms;
2445 ip = top->idef.iparams;
2446 if (blk_disre->nsub != 2 ||
2447 (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2449 gmx_incons("Number of disre sub-blocks not equal to 2");
2452 ndisre = blk_disre->sub[0].nr;
2453 if (ndisre != top->idef.il[F_DISRES].nr/3)
2455 gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2456 ndisre, top->idef.il[F_DISRES].nr/3);
2458 snew(pairleg, ndisre);
2459 for (i = 0; i < ndisre; i++)
2461 snew(pairleg[i], 30);
2464 gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j);
2465 gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k);
2466 sprintf(pairleg[i], "%d %s %d %s (%d)",
2467 resnr_j, anm_j, resnr_k, anm_k,
2468 ip[fa[3*i]].disres.label);
2470 set = select_it(ndisre, pairleg, &nset);
2472 for (i = 0; (i < nset); i++)
2475 sprintf(leg[2*i], "a %s", pairleg[set[i]]);
2476 snew(leg[2*i+1], 32);
2477 sprintf(leg[2*i+1], "i %s", pairleg[set[i]]);
2479 xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv);
2483 * Store energies for analysis afterwards...
2485 if (!bDisRe && !bDHDL && (fr->nre > 0))
2487 if (edat.nframes % 1000 == 0)
2489 srenew(time, edat.nframes+1000);
2491 time[edat.nframes] = fr->t;
2495 * Printing time, only when we do not want to skip frames
2497 if (!skip || teller % skip == 0)
2501 /*******************************************
2502 * D I S T A N C E R E S T R A I N T S
2503 *******************************************/
2507 float *disre_rt = blk_disre->sub[0].fval;
2508 float *disre_rm3tav = blk_disre->sub[1].fval;
2510 double *disre_rt = blk_disre->sub[0].dval;
2511 double *disre_rm3tav = blk_disre->sub[1].dval;
2514 print_time(out, fr->t);
2515 if (violaver == NULL)
2517 snew(violaver, ndisre);
2520 /* Subtract bounds from distances, to calculate violations */
2521 calc_violations(disre_rt, disre_rm3tav,
2522 nbounds, pair, bounds, violaver, &sumt, &sumaver);
2524 fprintf(out, " %8.4f %8.4f\n", sumaver, sumt);
2527 print_time(fp_pairs, fr->t);
2528 for (i = 0; (i < nset); i++)
2531 fprintf(fp_pairs, " %8.4f", mypow(disre_rm3tav[sss], minthird));
2532 fprintf(fp_pairs, " %8.4f", disre_rt[sss]);
2534 fprintf(fp_pairs, "\n");
2541 do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2544 /*******************************************
2546 *******************************************/
2553 /* We skip frames with single points (usually only the first frame),
2554 * since they would result in an average plot with outliers.
2558 print_time(out, fr->t);
2559 print1(out, bDp, fr->ener[set[0]].e);
2560 print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
2561 print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum));
2567 print_time(out, fr->t);
2571 for (i = 0; i < nset; i++)
2573 sum += fr->ener[set[i]].e;
2575 print1(out, bDp, sum/nmol-ezero);
2579 for (i = 0; (i < nset); i++)
2583 print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero);
2587 print1(out, bDp, fr->ener[set[i]].e);
2594 blk = find_block_id_enxframe(fr, enx_i, NULL);
2598 xdr_datatype dt = xdr_datatype_float;
2600 xdr_datatype dt = xdr_datatype_double;
2604 if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2606 gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2609 vals = blk->sub[0].fval;
2611 vals = blk->sub[0].dval;
2614 if (blk->sub[0].nr != (size_t)nor)
2616 gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2620 for (i = 0; i < nor; i++)
2622 orient[i] += vals[i];
2627 for (i = 0; i < nor; i++)
2629 odrms[i] += sqr(vals[i]-oobs[i]);
2634 fprintf(fort, " %10f", fr->t);
2635 for (i = 0; i < norsel; i++)
2637 fprintf(fort, " %g", vals[orsel[i]]);
2639 fprintf(fort, "\n");
2643 fprintf(fodt, " %10f", fr->t);
2644 for (i = 0; i < norsel; i++)
2646 fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]);
2648 fprintf(fodt, "\n");
2652 blk = find_block_id_enxframe(fr, enxORT, NULL);
2656 xdr_datatype dt = xdr_datatype_float;
2658 xdr_datatype dt = xdr_datatype_double;
2662 if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2664 gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2667 vals = blk->sub[0].fval;
2669 vals = blk->sub[0].dval;
2672 if (blk->sub[0].nr != (size_t)(nex*12))
2674 gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2675 blk->sub[0].nr/12, nex);
2677 fprintf(foten, " %10f", fr->t);
2678 for (i = 0; i < nex; i++)
2680 for (j = 0; j < (bOvec ? 12 : 3); j++)
2682 fprintf(foten, " %g", vals[i*12+j]);
2685 fprintf(foten, "\n");
2692 while (bCont && (timecheck == 0));
2694 fprintf(stderr, "\n");
2716 out = xvgropen(opt2fn("-ora", NFILE, fnm),
2717 "Average calculated orientations",
2718 "Restraint label", "", oenv);
2721 fprintf(out, "%s", orinst_sub);
2723 for (i = 0; i < nor; i++)
2725 fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr);
2731 out = xvgropen(opt2fn("-oda", NFILE, fnm),
2732 "Average restraint deviation",
2733 "Restraint label", "", oenv);
2736 fprintf(out, "%s", orinst_sub);
2738 for (i = 0; i < nor; i++)
2740 fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr-oobs[i]);
2746 out = xvgropen(opt2fn("-odr", NFILE, fnm),
2747 "RMS orientation restraint deviations",
2748 "Restraint label", "", oenv);
2751 fprintf(out, "%s", orinst_sub);
2753 for (i = 0; i < nor; i++)
2755 fprintf(out, "%5d %g\n", or_label[i], sqrt(odrms[i]/norfr));
2766 analyse_disre(opt2fn("-viol", NFILE, fnm),
2767 teller_disre, violaver, bounds, index, pair, nbounds, oenv);
2774 printf("\n\nWrote %d lambda values with %d samples as ",
2775 dh_lambdas, dh_samples);
2778 printf("%d dH histograms ", dh_hists);
2782 printf("%d dH data blocks ", dh_blocks);
2784 printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2789 gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2795 double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2796 analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
2797 bFee, bSum, opt2parg_bSet("-nmol", npargs, ppa),
2798 bVisco, opt2fn("-vis", NFILE, fnm),
2800 start_step, start_t, frame[cur].step, frame[cur].t,
2801 time, reftemp, &edat,
2802 nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
2806 calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat,
2810 if (opt2bSet("-f2", NFILE, fnm))
2812 fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm),
2813 reftemp, nset, set, leg, &edat, time, oenv);
2817 const char *nxy = "-nxy";
2819 do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2820 do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2821 do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
2822 do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
2823 do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
2824 do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
2825 do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
2826 do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
2827 do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);