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
int nbmin, int nbmax)
{
int i, j;
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
double NANO3;
enum {
eVol, eEnth, eTemp, eEtot, eNR
fprintf(fp,"Found %s data.\n",my_ener[i]);
*/ }
/* Compute it all! */
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;
/* Temperature */
tt = NOTSET;
/* Alpha, dcp */
if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
{
/* Alpha, dcp */
if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
{
+ 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++)
{
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);
{
fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol);
}
{
fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol);
}
- if (vvhh != NOTSET)
- {
- fprintf(fp, "vvhh = %10g (m^3 J)\n", vvhh);
- }