From a052b75ca8405fc167228fa409dfdc8337ba2aea Mon Sep 17 00:00:00 2001 From: David van der Spoel Date: Mon, 23 Dec 2013 21:42:08 +0100 Subject: [PATCH] Fixed precision in thermal expansion coefficient calc. Loss of accuracy was caused by different sampling of volume and enthalpy and as a result alpha was computed incorrectly. With the present "fix" the volume and enthalpy are both downsampled to what is written in the .edr file. The real fix would be to store the product of H and V in the .edr file, but that falls outside the 4.6 branch policy. Change-Id: I1be06d689002d7c9d6be92bf1e377912f0be1efd --- src/tools/gmx_energy.c | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/tools/gmx_energy.c b/src/tools/gmx_energy.c index cec4125cfd..54e911015c 100644 --- a/src/tools/gmx_energy.c +++ b/src/tools/gmx_energy.c @@ -980,7 +980,7 @@ static void calc_fluctuation_props(FILE *fp, int nbmin, int nbmax) { int i, j; - double vvhh, vv, v, h, hh2, vv2, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet; + double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet; double NANO3; enum { eVol, eEnth, eTemp, eEtot, eNR @@ -1010,7 +1010,7 @@ static void calc_fluctuation_props(FILE *fp, fprintf(fp,"Found %s data.\n",my_ener[i]); */ } /* Compute it all! */ - vvhh = alpha = kappa = cp = dcp = cv = NOTSET; + alpha = kappa = cp = dcp = cv = NOTSET; /* Temperature */ tt = NOTSET; @@ -1048,16 +1048,21 @@ static void calc_fluctuation_props(FILE *fp, /* Alpha, dcp */ if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset)) { - vvhh = 0; + double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver; + vh_sum = v_sum = h_sum = 0; for (j = 0; (j < edat->nframes); j++) { - v = edat->s[ii[eVol]].ener[j]*NANO3; - h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO; - vvhh += (v*h); + v = edat->s[ii[eVol]].ener[j]*NANO3; + h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO; + v_sum += v; + h_sum += h; + vh_sum += (v*h); } - vvhh /= edat->nframes; - alpha = (vvhh-vv*hh)/(vv*BOLTZMANN*tt*tt); - dcp = (vv*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa); + vh_aver = vh_sum / edat->nframes; + v_aver = v_sum / edat->nframes; + h_aver = h_sum / edat->nframes; + alpha = (vh_aver-v_aver*h_aver)/(v_aver*BOLTZMANN*tt*tt); + dcp = (v_aver*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa); } if (tt != NOTSET) @@ -1080,10 +1085,6 @@ static void calc_fluctuation_props(FILE *fp, { fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol); } - if (vvhh != NOTSET) - { - fprintf(fp, "vvhh = %10g (m^3 J)\n", vvhh); - } } if (vv != NOTSET) { -- 2.22.0