Merge origin/release-2019 into release-2020

[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_angle.cpp

diff --git a/src/gromacs/gmxana/gmx_angle.cpp b/src/gromacs/gmxana/gmx_angle.cpp
index 0a92c634d86bb4a2eb91a7e44024e3d772db76c5..e99bee5628e958d5ed9574e09c120adb03d3d753 100644 (file)
--- a/src/gromacs/gmxana/gmx_angle.cpp
+++ b/src/gromacs/gmxana/gmx_angle.cpp
@@ -46,6 +46,7 @@
 #include "gromacs/correlationfunctions/autocorr.h"
 #include "gromacs/fileio/trrio.h"
 #include "gromacs/fileio/xvgr.h"
+#include "gromacs/gmxana/angle_correction.h"
 #include "gromacs/gmxana/gmx_ana.h"
 #include "gromacs/gmxana/gstat.h"
 #include "gromacs/math/functions.h"
@@ -174,7 +175,6 @@ int gmx_g_angle(int argc, char* argv[])
     gmx_bool      bAver, bRb, bPeriodic, bFrac, /* calculate fraction too?  */
             bTrans,                             /* worry about transtions too? */
             bCorr;                              /* correlation function ? */
-    real     aver, aver2, aversig;              /* fraction trans dihedrals */
     double   tfrac = 0;
     char     title[256];
     real**   dih = nullptr; /* mega array with all dih. angles at all times*/
@@ -403,19 +403,40 @@ int gmx_g_angle(int argc, char* argv[])
     for (first = 0; (first < maxangstat - 1) && (angstat[first + 1] == 0); first++) {}
     for (last = maxangstat - 1; (last > 0) && (angstat[last - 1] == 0); last--) {}
 
-    aver = aver2 = 0;
-    for (i = 0; (i < nframes); i++)
-    {
-        aver += RAD2DEG * aver_angle[i];
-        aver2 += gmx::square(RAD2DEG * aver_angle[i]);
+    double aver = 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))
+    { /* It's better to re-calculate Std. Dev per sample */
+        real b_aver = aver_angle[0];
+        real b      = dih[0][0];
+        real delta;
+        for (int i = 0; (i < nframes); i++)
+        {
+            delta = correctRadianAngleRange(aver_angle[i] - b_aver);
+            b_aver += delta;
+            aver += b_aver;
+            for (int j = 0; (j < nangles); j++)
+            {
+                delta = correctRadianAngleRange(dih[j][i] - b);
+                b += delta;
+            }
+        }
+    }
+    else
+    { /* Incorrect  for Std. Dev. */
+        real delta, b_aver = aver_angle[0];
+        for (i = 0; (i < nframes); i++)
+        {
+            delta = correctRadianAngleRange(aver_angle[i] - b_aver);
+            b_aver += delta;
+            aver += b_aver;
+        }
     }
     aver /= nframes;
-    aver2 /= nframes;
-    aversig = std::sqrt(aver2 - gmx::square(aver));
-    printf("Found points in the range from %d to %d (max %d)\n", first, last, maxangstat);
-    printf(" < angle >  = %g\n", aver);
-    printf("< angle^2 > = %g\n", aver2);
-    printf("Std. Dev.   = %g\n", aversig);
+    double aversig = correctRadianAngleRange(aver);
+    aversig *= RAD2DEG;
+    aver *= RAD2DEG;
+    printf(" < angle >  = %g\n", aversig);
 
     if (mult == 3)
     {