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

diff --git a/src/tools/gmx_energy.c b/src/tools/gmx_energy.c
index 54e911015c67d685daa7e4ddc6e69d382848e48b..3323f35786b0e755d16dfabbfe5c436e1f278ddd 100644 (file)
--- 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)
         {