From 9696f6dd5b97b1ea92ec49c560fae551c5772912 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Mon, 2 Jun 2014 22:48:18 +0200 Subject: [PATCH] Fixed g_energy Einstein viscosity The g_energy -vis Einstein viscosity output was obviously incorrect. This bug has been introduced in version 4.5. Fixes #1516 Change-Id: I5ffc1a232f0c64769cc438c977b757e8d8b55b98 --- src/tools/gmx_energy.c | 67 +++++++++++++++++++++++++----------------- 1 file changed, 40 insertions(+), 27 deletions(-) diff --git a/src/tools/gmx_energy.c b/src/tools/gmx_energy.c index 54e911015c..3323f35786 100644 --- a/src/tools/gmx_energy.c +++ b/src/tools/gmx_energy.c @@ -570,22 +570,16 @@ static void analyse_disre(const char *voutfn, int nframes, } static void einstein_visco(const char *fn, const char *fni, int nsets, - int nframes, real **sum, - real V, real T, int nsteps, double time[], + int nint, real **eneint, + real V, real T, double dt, const output_env_t oenv) { FILE *fp0, *fp1; double av[4], avold[4]; - double fac, dt, di; + double fac, di; int i, j, m, nf4; - if (nframes < 1) - { - return; - } - - dt = (time[1]-time[0]); - nf4 = nframes/4+1; + nf4 = nint/4 + 1; for (i = 0; i <= nsets; i++) { @@ -595,33 +589,32 @@ static void einstein_visco(const char *fn, const char *fni, int nsets, "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv); fp1 = xvgropen(fn, "Shear viscosity using Einstein relation", "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv); - for (i = 1; i < nf4; i++) + for (i = 0; i < nf4; i++) { - fac = dt*nframes/nsteps; for (m = 0; m <= nsets; m++) { av[m] = 0; } - for (j = 0; j < nframes-i; j++) + for (j = 0; j < nint - i; j++) { for (m = 0; m < nsets; m++) { - di = sqr(fac*(sum[m][j+i]-sum[m][j])); + di = sqr(eneint[m][j+i] - eneint[m][j]); av[m] += di; av[nsets] += di/nsets; } } /* Convert to SI for the viscosity */ - fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i); - fprintf(fp0, "%10g", time[i]-time[0]); + fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nint - i); + fprintf(fp0, "%10g", i*dt); for (m = 0; (m <= nsets); m++) { av[m] = fac*av[m]; fprintf(fp0, " %10g", av[m]); } fprintf(fp0, "\n"); - fprintf(fp1, "%10g", 0.5*(time[i]+time[i-1])-time[0]); + fprintf(fp1, "%10g", (i + 0.5)*dt); for (m = 0; (m <= nsets); m++) { fprintf(fp1, " %10g", (av[m]-avold[m])/dt); @@ -1348,7 +1341,7 @@ static void analyse_ener(gmx_bool bCorr, const char *corrfn, const char* leg[] = { "Shear", "Bulk" }; real factor; real **eneset; - real **enesum; + real **eneint; /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */ @@ -1360,11 +1353,6 @@ static void analyse_ener(gmx_bool bCorr, const char *corrfn, { snew(eneset[i], edat->nframes); } - snew(enesum, 3); - for (i = 0; i < 3; i++) - { - snew(enesum[i], edat->nframes); - } for (i = 0; (i < edat->nframes); i++) { eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]); @@ -1375,13 +1363,32 @@ static void analyse_ener(gmx_bool bCorr, const char *corrfn, eneset[j][i] = edat->s[j].ener[i]; } eneset[11][i] -= Pres; - enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum); - enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum); - enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum); } + /* Determine integrals of the off-diagonal pressure elements */ + snew(eneint, 3); + for (i = 0; i < 3; i++) + { + snew(eneint[i], edat->nframes + 1); + } + eneint[0][0] = 0; + eneint[1][0] = 0; + eneint[2][0] = 0; + for (i = 0; i < edat->nframes; i++) + { + eneint[0][i+1] = eneint[0][i] + 0.5*(edat->s[1].es[i].sum + edat->s[3].es[i].sum)*Dt/edat->points[i]; + eneint[1][i+1] = eneint[1][i] + 0.5*(edat->s[2].es[i].sum + edat->s[6].es[i].sum)*Dt/edat->points[i]; + eneint[2][i+1] = eneint[2][i] + 0.5*(edat->s[5].es[i].sum + edat->s[7].es[i].sum)*Dt/edat->points[i]; + } + einstein_visco("evisco.xvg", "eviscoi.xvg", - 3, edat->nframes, enesum, Vaver, Temp, nsteps, time, oenv); + 3, edat->nframes+1, eneint, Vaver, Temp, Dt, oenv); + + for (i = 0; i < 3; i++) + { + sfree(eneint[i]); + } + sfree(eneint); /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/ /* Do it for shear viscosity */ @@ -1415,6 +1422,12 @@ static void analyse_ener(gmx_bool bCorr, const char *corrfn, fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk); } ffclose(fp); + + for (i = 0; i < 12; i++) + { + sfree(eneset[i]); + } + sfree(eneset); } else if (bCorr) { -- 2.22.0