Merge origin/release-2019 into release-2020

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

diff --git a/src/gromacs/gmxana/anadih.cpp b/src/gromacs/gmxana/anadih.cpp
index b76f99dc1aed9e2bf31e3bcd5ccedf66b7595668..bfe16bede134eead6646469ce52e79c2d1819237 100644 (file)
--- a/src/gromacs/gmxana/anadih.cpp
+++ b/src/gromacs/gmxana/anadih.cpp
@@ -45,6 +45,7 @@
 #include "gromacs/fileio/confio.h"
 #include "gromacs/fileio/trxio.h"
 #include "gromacs/fileio/xvgr.h"
+#include "gromacs/gmxana/angle_correction.h"
 #include "gromacs/gmxana/gstat.h"
 #include "gromacs/listed_forces/bonded.h"
 #include "gromacs/math/functions.h"
@@ -833,7 +834,7 @@ void read_ang_dih(const char*             trj_fn,
     t_trxstatus*  status;
     int           i, angind, total, teller;
     int           nangles, n_alloc;
-    real          t, fraction, pifac, aa, angle;
+    real          t, fraction, pifac, angle;
     real*         angles[2];
     matrix        box;
     rvec*         x;
@@ -946,9 +947,16 @@ void read_ang_dih(const char*             trj_fn,
         }
 
         /* Average angles */
-        aa = 0;
+        double aa = 0;
         for (i = 0; (i < nangles); i++)
         {
+            if (!bAngles && i > 0)
+            {
+                real diffa     = angles[cur][i] - angles[cur][i - 1];
+                diffa          = correctRadianAngleRange(diffa);
+                angles[cur][i] = angles[cur][i - 1] + diffa;
+            }
+
             aa = aa + angles[cur][i];
 
             /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
@@ -960,15 +968,7 @@ void read_ang_dih(const char*             trj_fn,
             angle = angles[cur][i];
             if (!bAngles)
             {
-                while (angle < -M_PI)
-                {
-                    angle += 2 * M_PI;
-                }
-                while (angle >= M_PI)
-                {
-                    angle -= 2 * M_PI;
-                }
-
+                angle = correctRadianAngleRange(angle);
                 angle += M_PI;
             }
 
@@ -994,14 +994,22 @@ void read_ang_dih(const char*             trj_fn,
         }
 
         /* average over all angles */
-        (*aver_angle)[teller] = (aa / nangles);
+        aa                    = correctRadianAngleRange(aa / nangles);
+        (*aver_angle)[teller] = (aa);
 
         /* this copies all current dih. angles to dih[i], teller is frame */
         if (bSaveAll)
         {
             for (i = 0; i < nangles; i++)
             {
-                dih[i][teller] = angles[cur][i];
+                if (!bAngles)
+                {
+                    dih[i][teller] = correctRadianAngleRange(angles[cur][i]);
+                }
+                else
+                {
+                    dih[i][teller] = angles[cur][i];
+                }
             }
         }