Fix norm calculation in gmx_spatial

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

diff --git a/src/gromacs/gmxana/gmx_spatial.cpp b/src/gromacs/gmxana/gmx_spatial.cpp
index 0ec43ffafda9f4a68fdc4bbe948b17d00c0755d8..90eda00f9f476496105fd0c7f04c78c8f5097516 100644 (file)
--- a/src/gromacs/gmxana/gmx_spatial.cpp
+++ b/src/gromacs/gmxana/gmx_spatial.cpp
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 2007,2008,2009,2010,2011 by the GROMACS development team.
  * Copyright (c) 2012,2013,2014,2015,2017 by the GROMACS development team.
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -52,6 +52,7 @@
 #include "gromacs/utility/arraysize.h"
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 
 static const double bohr =
@@ -190,8 +191,8 @@ int gmx_spatial(int argc, char* argv[])
 
     /* This is the routine responsible for adding default options,
      * calling the X/motif interface, etc. */
-    if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa,
-                           asize(desc), desc, asize(bugs), bugs, &oenv))
+    if (!parse_common_args(
+                &argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs, &oenv))
     {
         return 0;
     }
@@ -284,10 +285,17 @@ int gmx_spatial(int argc, char* argv[])
             {
                 printf("There was an item outside of the allocated memory. Increase the value "
                        "given with the -nab option.\n");
-                printf("Memory was allocated for [%f,%f,%f]\tto\t[%f,%f,%f]\n", MINBIN[XX],
-                       MINBIN[YY], MINBIN[ZZ], MAXBIN[XX], MAXBIN[YY], MAXBIN[ZZ]);
-                printf("Memory was required for [%f,%f,%f]\n", fr.x[index[i]][XX],
-                       fr.x[index[i]][YY], fr.x[index[i]][ZZ]);
+                printf("Memory was allocated for [%f,%f,%f]\tto\t[%f,%f,%f]\n",
+                       MINBIN[XX],
+                       MINBIN[YY],
+                       MINBIN[ZZ],
+                       MAXBIN[XX],
+                       MAXBIN[YY],
+                       MAXBIN[ZZ]);
+                printf("Memory was required for [%f,%f,%f]\n",
+                       fr.x[index[i]][XX],
+                       fr.x[index[i]][YY],
+                       fr.x[index[i]][ZZ]);
                 exit(1);
             }
             x = static_cast<int>(std::ceil((fr.x[index[i]][XX] - MINBIN[XX]) / rBINWIDTH));
@@ -341,16 +349,15 @@ int gmx_spatial(int argc, char* argv[])
     flp = gmx_ffopen("grid.cube", "w");
     fprintf(flp, "Spatial Distribution Function\n");
     fprintf(flp, "test\n");
-    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", nidxp,
+    fprintf(flp,
+            "%5d%12.6f%12.6f%12.6f\n",
+            nidxp,
             (MINBIN[XX] + (minx + iIGNOREOUTER) * rBINWIDTH) * 10. / bohr,
             (MINBIN[YY] + (miny + iIGNOREOUTER) * rBINWIDTH) * 10. / bohr,
             (MINBIN[ZZ] + (minz + iIGNOREOUTER) * rBINWIDTH) * 10. / bohr);
-    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", maxx - minx + 1 - (2 * iIGNOREOUTER),
-            rBINWIDTH * 10. / bohr, 0., 0.);
-    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", maxy - miny + 1 - (2 * iIGNOREOUTER), 0.,
-            rBINWIDTH * 10. / bohr, 0.);
-    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", maxz - minz + 1 - (2 * iIGNOREOUTER), 0., 0.,
-            rBINWIDTH * 10. / bohr);
+    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", maxx - minx + 1 - (2 * iIGNOREOUTER), rBINWIDTH * 10. / bohr, 0., 0.);
+    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", maxy - miny + 1 - (2 * iIGNOREOUTER), 0., rBINWIDTH * 10. / bohr, 0.);
+    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", maxz - minz + 1 - (2 * iIGNOREOUTER), 0., 0., rBINWIDTH * 10. / bohr);
     for (i = 0; i < nidxp; i++)
     {
         v = 2;
@@ -374,8 +381,13 @@ int gmx_spatial(int argc, char* argv[])
         {
             v = 16;
         }
-        fprintf(flp, "%5d%12.6f%12.6f%12.6f%12.6f\n", v, 0., fr.x[indexp[i]][XX] * 10.0 / bohr,
-                fr.x[indexp[i]][YY] * 10.0 / bohr, fr.x[indexp[i]][ZZ] * 10.0 / bohr);
+        fprintf(flp,
+                "%5d%12.6f%12.6f%12.6f%12.6f\n",
+                v,
+                0.,
+                fr.x[indexp[i]][XX] * 10.0 / bohr,
+                fr.x[indexp[i]][YY] * 10.0 / bohr,
+                fr.x[indexp[i]][ZZ] * 10.0 / bohr);
     }
 
     tot = 0;
@@ -445,7 +457,8 @@ int gmx_spatial(int argc, char* argv[])
             * (maxz - minz + 1 - (2 * iIGNOREOUTER));
     if (bCALCDIV)
     {
-        norm = static_cast<double>(numcu * numfr) / tot;
+        norm = double(numcu) * numfr / tot;
+        GMX_ASSERT(norm >= 0, "The norm should be non-negative.");
     }
     else
     {
@@ -481,14 +494,18 @@ int gmx_spatial(int argc, char* argv[])
     if (bCALCDIV)
     {
         printf("Counts per frame in all %d cubes divided by %le\n", numcu, 1.0 / norm);
-        printf("Normalized data: average %le, min %le, max %le\n", 1.0, minval * norm / numfr,
+        printf("Normalized data: average %le, min %le, max %le\n",
+               1.0,
+               minval * norm / numfr,
                maxval * norm / numfr);
     }
     else
     {
         printf("grid.cube contains counts per frame in all %d cubes\n", numcu);
-        printf("Raw data: average %le, min %le, max %le\n", 1.0 / norm,
-               static_cast<double>(minval) / numfr, static_cast<double>(maxval) / numfr);
+        printf("Raw data: average %le, min %le, max %le\n",
+               1.0 / norm,
+               static_cast<double>(minval) / numfr,
+               static_cast<double>(maxval) / numfr);
     }
 
     return 0;