From: Berk Hess <hess@kth.se>
Date: Mon, 22 Jun 2015 14:50:53 +0000 (+0200)
Subject: Fix g_energy average/RMSD bug
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=83283df8937e8419452929b76dbedac9e0e90dd5;p=alexxy%2Fgromacs.git

Fix g_energy average/RMSD bug

Made g_energy produce correct output for energy files from continued
and appended runs with nstcalcenergy=nstenergy. In that case a count
in the energy file is/was incorrect, but that entry is actually not
necessary for determining the average and RMSD of energy terms. The
RMSD would be NaN and the average would be off in the last decimal.

Fixes #1342.

Change-Id: I82007bfe508023e1c1e17366e10f76bc4470d238
---

diff --git a/src/tools/gmx_energy.c b/src/tools/gmx_energy.c
index 7341ddbd34..df155c2a0d 100644
--- a/src/tools/gmx_energy.c
+++ b/src/tools/gmx_energy.c
@@ -87,6 +87,7 @@ typedef struct {
     int             *steps;
     int             *points;
     enerdat_t       *s;
+    gmx_bool         bHaveSums;
 } enerdata_t;
 
 static double mypow(double x, double y)
@@ -703,7 +704,7 @@ static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
     {
         ed             = &edat->s[i];
         ed->bExactStat = FALSE;
-        if (edat->npoints > 0)
+        if (edat->bHaveSums)
         {
             /* All energy file sum entries 0 signals no exact sums.
              * But if all energy values are 0, we still have exact sums.
@@ -1172,7 +1173,7 @@ static void analyse_ener(gmx_bool bCorr, const char *corrfn,
             esum = calc_sum(nset, edat, nbmin, nbmax);
         }
 
-        if (edat->npoints == 0)
+        if (!edat->bHaveSums)
         {
             nexact    = 0;
             nnotexact = nset;
@@ -2321,12 +2322,13 @@ int gmx_energy(int argc, char *argv[])
     }
 
     /* Initiate energies and set them to zero */
-    edat.nsteps  = 0;
-    edat.npoints = 0;
-    edat.nframes = 0;
-    edat.step    = NULL;
-    edat.steps   = NULL;
-    edat.points  = NULL;
+    edat.nsteps    = 0;
+    edat.npoints   = 0;
+    edat.nframes   = 0;
+    edat.step      = NULL;
+    edat.steps     = NULL;
+    edat.points    = NULL;
+    edat.bHaveSums = TRUE;
     snew(edat.s, nset);
 
     /* Initiate counters */
@@ -2404,41 +2406,43 @@ int gmx_energy(int argc, char *argv[])
                 else
                 {
                     edat.steps[nfr] = fr->nsteps;
+
+                    if (fr->nsum <= 1)
                     {
-                        if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
+                        /* mdrun only calculated the energy at energy output
+                         * steps. We don't need to check step intervals.
+                         */
+                        edat.points[nfr] = 1;
+                        for (i = 0; i < nset; i++)
                         {
-                            if (fr->nsum <= 1)
-                            {
-                                edat.points[nfr] = 1;
-                                for (i = 0; i < nset; i++)
-                                {
-                                    sss                    = set[i];
-                                    edat.s[i].es[nfr].sum  = fr->ener[sss].e;
-                                    edat.s[i].es[nfr].sum2 = 0;
-                                }
-                                edat.npoints += 1;
-                            }
-                            else
-                            {
-                                edat.points[nfr] = fr->nsum;
-                                for (i = 0; i < nset; i++)
-                                {
-                                    sss                    = set[i];
-                                    edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
-                                    edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
-                                }
-                                edat.npoints += fr->nsum;
-                            }
+                            sss                    = set[i];
+                            edat.s[i].es[nfr].sum  = fr->ener[sss].e;
+                            edat.s[i].es[nfr].sum2 = 0;
                         }
-                        else
+                        edat.npoints  += 1;
+                        edat.bHaveSums = FALSE;
+                    }
+                    else if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
+                    {
+                        /* We have statistics  to the previous frame */
+                        edat.points[nfr] = fr->nsum;
+                        for (i = 0; i < nset; i++)
                         {
-                            /* The interval does not match fr->nsteps:
-                             * can not do exact averages.
-                             */
-                            edat.npoints = 0;
+                            sss                    = set[i];
+                            edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
+                            edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
                         }
-                        edat.nsteps = fr->step - start_step + 1;
+                        edat.npoints += fr->nsum;
                     }
+                    else
+                    {
+                        /* The interval does not match fr->nsteps:
+                         * can not do exact averages.
+                         */
+                        edat.bHaveSums = FALSE;
+                    }
+
+                    edat.nsteps = fr->step - start_step + 1;
                 }
                 for (i = 0; i < nset; i++)
                 {