Merge branch 'master' into pygromacs
[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_sham.cpp
similarity index 94%
rename from src/gromacs/gmxana/gmx_sham.c
rename to src/gromacs/gmxana/gmx_sham.cpp
index 87b28e124d386a724bbdab8e63e377bab0e788d7..f24788810c2121678d1a8407a5e5264d34f180f4 100644 (file)
  */
 #include "gmxpre.h"
 
-#include <math.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+
+#include <algorithm>
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/fileio/matio.h"
@@ -229,7 +231,7 @@ static gmx_inline
 void print_minimum(FILE *fp, int num, const t_minimum *min)
 {
     fprintf(fp,
-            "Minimum %d at index " "%"GMX_PRId64 " energy %10.3f\n",
+            "Minimum %d at index " "%" GMX_PRId64 " energy %10.3f\n",
             num, min->index, min->ener);
 }
 
@@ -428,9 +430,9 @@ static void do_sham(const char *fn, const char *ndx,
     real       **PP, *W, *E, **WW, **EE, *S, **SS, *M, *bE;
     rvec         xxx;
     char        *buf;
-    double      *bfac, efac, bref, Pmax, Wmin, Wmax, Winf, Emin, Emax, Einf, Smin, Smax, Sinf, Mmin, Mmax;
+    double      *bfac, efac, bref, Pmax, Wmin, Wmax, Winf, Emin, Emax, Einf, Smin, Smax, Sinf;
     real        *delta;
-    int          i, j, k, imin, len, index, d, *nbin, *bindex, bi;
+    int          i, j, k, imin, len, index, *nbin, *bindex, bi;
     int         *nxyz, maxbox;
     t_blocka    *b;
     gmx_bool     bOutside;
@@ -451,8 +453,8 @@ static void do_sham(const char *fn, const char *ndx,
         min_eig[i] = max_eig[i] = eig[i][0];
         for (j = 0; (j < n); j++)
         {
-            min_eig[i] = min(min_eig[i], eig[i][j]);
-            max_eig[i] = max(max_eig[i], eig[i][j]);
+            min_eig[i] = std::min(min_eig[i], eig[i][j]);
+            max_eig[i] = std::max(max_eig[i], eig[i][j]);
             delta[i]   = (max_eig[i]-min_eig[i])/(2.0*ibox[i]);
         }
         /* Add some extra space, half a bin on each side, unless the
@@ -501,7 +503,7 @@ static void do_sham(const char *fn, const char *ndx,
             {
                 bE[j] = (bref - 1/(BOLTZ*enerT[1][j]))*enerT[0][j];
             }
-            Emin  = min(Emin, bE[j]);
+            Emin  = std::min(Emin, static_cast<double>(bE[j]));
         }
     }
     else
@@ -531,7 +533,7 @@ static void do_sham(const char *fn, const char *ndx,
         bOutside = FALSE;
         for (i = 0; (i < neig); i++)
         {
-            nxyz[i] = bfac[i]*(eig[i][j]-min_eig[i]);
+            nxyz[i] = static_cast<int>(bfac[i]*(eig[i][j]-min_eig[i]));
             if (nxyz[i] < 0 || nxyz[i] >= ibox[i])
             {
                 bOutside = TRUE;
@@ -544,7 +546,7 @@ static void do_sham(const char *fn, const char *ndx,
             /* Compute the exponential factor */
             if (enerT)
             {
-                efac = exp(-bE[j]+Emin);
+                efac = std::exp(-bE[j]+Emin);
             }
             else
             {
@@ -563,7 +565,7 @@ static void do_sham(const char *fn, const char *ndx,
                 }
                 else if (idim[i] == -1)
                 {
-                    efac /= sin(DEG2RAD*eig[i][j]);
+                    efac /= std::sin(DEG2RAD*eig[i][j]);
                 }
             }
             /* Update the probability */
@@ -592,16 +594,16 @@ static void do_sham(const char *fn, const char *ndx,
     {
         if (P[i] != 0)
         {
-            Pmax = max(P[i], Pmax);
-            W[i] = -BOLTZ*Tref*log(P[i]);
+            Pmax = std::max(P[i], Pmax);
+            W[i] = -BOLTZ*Tref*std::log(P[i]);
             if (W[i] < Wmin)
             {
                 Wmin = W[i];
                 imin = i;
             }
-            Emin = min(E[i], Emin);
-            Emax = max(E[i], Emax);
-            Wmax = max(W[i], Wmax);
+            Emin = std::min(static_cast<double>(E[i]), Emin);
+            Emax = std::max(static_cast<double>(E[i]), Emax);
+            Wmax = std::max(static_cast<double>(W[i]), Wmax);
         }
     }
     if (pmax > 0)
@@ -681,12 +683,12 @@ static void do_sham(const char *fn, const char *ndx,
     snew(axis_x, ibox[0]+1);
     snew(axis_y, ibox[1]+1);
     snew(axis_z, ibox[2]+1);
-    maxbox = max(ibox[0], max(ibox[1], ibox[2]));
+    maxbox = std::max(ibox[0], std::max(ibox[1], ibox[2]));
     snew(PP, maxbox*maxbox);
     snew(WW, maxbox*maxbox);
     snew(EE, maxbox*maxbox);
     snew(SS, maxbox*maxbox);
-    for (i = 0; (i < min(neig, 3)); i++)
+    for (i = 0; (i < std::min(neig, 3)); i++)
     {
         switch (i)
         {
@@ -776,9 +778,9 @@ static void do_sham(const char *fn, const char *ndx,
                 WW[i][j] = W[index3(ibox, i, j, nxyz[ZZ])];
             }
         }
-        snew(buf, strlen(xpm)+4);
+        snew(buf, std::strlen(xpm)+4);
         sprintf(buf, "%s", xpm);
-        sprintf(&buf[strlen(xpm)-4], "12.xpm");
+        sprintf(&buf[std::strlen(xpm)-4], "12.xpm");
         fp = gmx_ffopen(buf, "w");
         write_xpm(fp, flags, "Gibbs Energy Landscape", "W (kJ/mol)", "PC1", "PC2",
                   ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels);
@@ -790,7 +792,7 @@ static void do_sham(const char *fn, const char *ndx,
                 WW[i][j] = W[index3(ibox, i, nxyz[YY], j)];
             }
         }
-        sprintf(&buf[strlen(xpm)-4], "13.xpm");
+        sprintf(&buf[std::strlen(xpm)-4], "13.xpm");
         fp = gmx_ffopen(buf, "w");
         write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC1", "PC3",
                   ibox[0], ibox[2], axis_x, axis_z, WW, 0, gmax, rlo, rhi, &nlevels);
@@ -802,7 +804,7 @@ static void do_sham(const char *fn, const char *ndx,
                 WW[i][j] = W[index3(ibox, nxyz[XX], i, j)];
             }
         }
-        sprintf(&buf[strlen(xpm)-4], "23.xpm");
+        sprintf(&buf[std::strlen(xpm)-4], "23.xpm");
         fp = gmx_ffopen(buf, "w");
         write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC2", "PC3",
                   ibox[1], ibox[2], axis_y, axis_z, WW, 0, gmax, rlo, rhi, &nlevels);
@@ -840,11 +842,11 @@ static void ehisto(const char *fh, int n, real **enerT, const output_env_t oenv)
             T[nbin]   = enerT[1][j];
             nbin++;
         }
-        bmin = min(enerT[0][j], bmin);
-        bmax = max(enerT[0][j], bmax);
+        bmin = std::min(enerT[0][j], bmin);
+        bmax = std::max(enerT[0][j], bmax);
     }
     bwidth  = 1.0;
-    blength = (bmax - bmin)/bwidth + 2;
+    blength = static_cast<int>((bmax - bmin)/bwidth + 2);
     snew(histo, nbin);
     for (i = 0; (i < nbin); i++)
     {
@@ -852,7 +854,7 @@ static void ehisto(const char *fh, int n, real **enerT, const output_env_t oenv)
     }
     for (j = 0; (j < n); j++)
     {
-        k = (enerT[0][j]-bmin)/bwidth;
+        k = static_cast<int>((enerT[0][j]-bmin)/bwidth);
         histo[bindex[j]][k]++;
     }
     fp = xvgropen(fh, "Energy distribution", "E (kJ/mol)", "", oenv);
@@ -911,15 +913,14 @@ int gmx_sham(int argc, char *argv[])
         "is the natural quantity to use, as it will produce bins of the same",
         "volume."
     };
-    static real        tb        = -1, te = -1, frac = 0.5, filtlen = 0;
-    static gmx_bool    bHaveT    = TRUE, bDer = FALSE, bSubAv = TRUE, bAverCorr = FALSE, bXYdy = FALSE;
-    static gmx_bool    bEESEF    = FALSE, bEENLC = FALSE, bEeFitAc = FALSE, bPower = FALSE;
-    static gmx_bool    bShamEner = TRUE, bSham = TRUE;
+    static real        tb        = -1, te = -1;
+    static gmx_bool    bHaveT    = TRUE, bDer = FALSE;
+    static gmx_bool    bSham     = TRUE;
     static real        Tref      = 298.15, pmin = 0, ttol = 0, pmax = 0, gmax = 0, emin = 0, emax = 0;
     static rvec        nrdim     = {1, 1, 1};
     static rvec        nrbox     = {32, 32, 32};
     static rvec        xmin      = {0, 0, 0}, xmax = {1, 1, 1};
-    static int         nsets_in  = 1, nb_min = 4, resol = 10, nlevels = 25;
+    static int         nsets_in  = 1, nlevels = 25;
     t_pargs            pa[]      = {
         { "-time",    FALSE, etBOOL, {&bHaveT},
           "Expect a time in the input" },
@@ -960,11 +961,9 @@ int gmx_sham(int argc, char *argv[])
     };
 #define NPA asize(pa)
 
-    FILE           *out;
-    int             n, e_n, nlast, s, nset, e_nset, d_nset, i, j = 0, *idim, *ibox;
-    real          **val, **et_val, *t, *e_t, e_dt, d_dt, dt, tot, error;
+    int             n, e_n, nset, e_nset = 0, i, *idim, *ibox;
+    real          **val, **et_val, *t, *e_t, e_dt, dt;
     real           *rmin, *rmax;
-    double         *av, *sig, cum1, cum2, cum3, cum4, db;
     const char     *fn_ge, *fn_ene;
     output_env_t    oenv;
     gmx_int64_t     num_grid_points;
@@ -1045,21 +1044,21 @@ int gmx_sham(int argc, char *argv[])
         ehisto(opt2fn("-histo", NFILE, fnm), e_n, et_val, oenv);
     }
 
-    snew(idim, max(3, nset));
-    snew(ibox, max(3, nset));
-    snew(rmin, max(3, nset));
-    snew(rmax, max(3, nset));
-    for (i = 0; (i < min(3, nset)); i++)
+    snew(idim, std::max(3, nset));
+    snew(ibox, std::max(3, nset));
+    snew(rmin, std::max(3, nset));
+    snew(rmax, std::max(3, nset));
+    for (i = 0; (i < std::min(3, nset)); i++)
     {
-        idim[i] = nrdim[i];
-        ibox[i] = nrbox[i];
+        idim[i] = static_cast<int>(nrdim[i]);
+        ibox[i] = static_cast<int>(nrbox[i]);
         rmin[i] = xmin[i];
         rmax[i] = xmax[i];
     }
     for (; (i < nset); i++)
     {
-        idim[i] = nrdim[2];
-        ibox[i] = nrbox[2];
+        idim[i] = static_cast<int>(nrdim[2]);
+        ibox[i] = static_cast<int>(nrbox[2]);
         rmin[i] = xmin[2];
         rmax[i] = xmax[2];
     }