Remove problematic output of gmx angle

It could happen that the calculation of the standard deviation
for angles caused a divide by zero error for empty populations.
Because this standard deviation of per frame averages over all
frames was anyhow meaningless, it has been removed.

Refs #3206

Change-Id: Icf541a8b3bef502da9a22b55a8e8154b2d4fc8f8

diff --git a/docs/release-notes/2019/2019.5.rst b/docs/release-notes/2019/2019.5.rst
index 4b5a496923c2777de2de50a3f3c53cde6102339b..57a041e63ac5668959d2969d73de70e29f7eef27 100644 (file)
--- a/docs/release-notes/2019/2019.5.rst
+++ b/docs/release-notes/2019/2019.5.rst
@@ -96,6 +96,16 @@ when close to the -180 or 180 degree boundary.
 :issue:`3225`
 
 
+Remove problematic output of gmx angle tool
+"""""""""""""""""""""""""""""""""""""""""""
+
+It could happen that the calculation of the standard deviation
+for angles caused a divide by zero error for empty populations.
+Because this standard deviation was meaningless, it has been
+removed.
+
+:issue:`3206`
+
 Fixes that affect portability
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

diff --git a/src/gromacs/gmxana/gmx_angle.cpp b/src/gromacs/gmxana/gmx_angle.cpp
index 4fa39ad4f129962be4121bc4b7ee152c71669575..b8f5da3b7a7bfb6f90666c2d2a34d0408601a2a9 100644 (file)
--- a/src/gromacs/gmxana/gmx_angle.cpp
+++ b/src/gromacs/gmxana/gmx_angle.cpp
@@ -165,11 +165,11 @@ int gmx_g_angle(int argc, char *argv[])
                        bFrac,               /* calculate fraction too?  */
                        bTrans,              /* worry about transtions too? */
                        bCorr;               /* correlation function ? */
-    double            tfrac = 0;
-    char              title[256];
-    real            **dih = nullptr; /* mega array with all dih. angles at all times*/
-    real             *time, *trans_frac, *aver_angle;
-    t_filenm          fnm[] = {
+    double             tfrac = 0;
+    char               title[256];
+    real             **dih = nullptr; /* mega array with all dih. angles at all times*/
+    real              *time, *trans_frac, *aver_angle;
+    t_filenm           fnm[] = {
         { efTRX, "-f", nullptr,  ffREAD  },
         { efNDX, nullptr, "angle",  ffREAD  },
         { efXVG, "-od", "angdist",  ffWRITE },
@@ -181,9 +181,9 @@ int gmx_g_angle(int argc, char *argv[])
         { efTRR, "-or", nullptr,       ffOPTWR }
     };
 #define NFILE asize(fnm)
-    int               npargs;
-    t_pargs          *ppa;
-    gmx_output_env_t *oenv;
+    int                npargs;
+    t_pargs           *ppa;
+    gmx_output_env_t  *oenv;
 
     npargs = asize(pa);
     ppa    = add_acf_pargs(&npargs, pa);
@@ -412,7 +412,6 @@ int gmx_g_angle(int argc, char *argv[])
     }
 
     double aver  = 0;
-    double aver2 = 0;
     printf("Found points in the range from %d to %d (max %d)\n",
            first, last, maxangstat);
     if (bTrans || bCorr  || bALL || opt2bSet("-or", NFILE, fnm))
@@ -429,10 +428,8 @@ int gmx_g_angle(int argc, char *argv[])
             {
                 delta  = correctRadianAngleRange(dih[j][i] - b);
                 b     += delta;
-                aver2 += gmx::square(b);
             }
         }
-        aver2 /= nangles;
     }
     else
     {   /* Incorrect  for Std. Dev. */
@@ -442,24 +439,13 @@ int gmx_g_angle(int argc, char *argv[])
             delta   = correctRadianAngleRange(aver_angle[i] - b_aver);
             b_aver += delta;
             aver   += b_aver;
-            aver2  += gmx::square(b_aver);
         }
     }
     aver  /= nframes;
-    aver2 /= nframes/(RAD2DEG * RAD2DEG);
     double aversig = correctRadianAngleRange(aver);
     aversig *= RAD2DEG;
     aver    *= RAD2DEG;
     printf(" < angle >  = %g\n", aversig);
-    printf("< angle^2 > = %g\n", aver2);
-    if ((aversig = aver2-gmx::square(aver)) > 0)
-    {
-        printf("Std. Dev.   = %g\n", std::sqrt(aversig));
-    }
-    else
-    {
-        printf("Std. Dev. - can't be calculated\n");
-    }
 
     if (mult == 3)
     {