Merge branch 'master' into pygromacs
[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_analyze.cpp
similarity index 94%
rename from src/gromacs/gmxana/gmx_analyze.c
rename to src/gromacs/gmxana/gmx_analyze.cpp
index 5d1c8c4a59d4a91a5955140c6e329d85a4415bde..c38228272764f6d0e18daf39be9e4767c2fa1b8f 100644 (file)
@@ -36,9 +36,9 @@
  */
 #include "gmxpre.h"
 
-#include <math.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/correlationfunctions/autocorr.h"
@@ -80,7 +80,7 @@ static void power_fit(int n, int nset, real **val, real *t)
         {
             if (t[0] > 0)
             {
-                x[i] = log(t[i]);
+                x[i] = std::log(t[i]);
             }
         }
     }
@@ -89,16 +89,15 @@ static void power_fit(int n, int nset, real **val, real *t)
         fprintf(stdout, "First time is not larger than 0, using index number as time for power fit\n");
         for (i = 0; i < n; i++)
         {
-            x[i] = gmx_log1p(i);
+            x[i] = gmx_log1p(static_cast<real>(i));
         }
     }
 
     for (s = 0; s < nset; s++)
     {
-        i = 0;
         for (i = 0; i < n && val[s][i] >= 0; i++)
         {
-            y[i] = log(val[s][i]);
+            y[i] = std::log(val[s][i]);
         }
         if (i < n)
         {
@@ -106,7 +105,7 @@ static void power_fit(int n, int nset, real **val, real *t)
         }
         lsq_y_ax_b(i, x, y, &a, &b, &r, &quality);
         fprintf(stdout, "Power fit set %3d:  error %.3f  a %g  b %g\n",
-                s+1, quality, a, exp(b));
+                s+1, quality, a, std::exp(b));
     }
 
     sfree(y);
@@ -130,7 +129,7 @@ static real cosine_content(int nhp, int n, real *y)
     yyint   = 0;
     for (i = 0; i < n; i++)
     {
-        cosyint += cos(fac*i)*y[i];
+        cosyint += std::cos(fac*i)*y[i];
         yyint   += y[i]*y[i];
     }
 
@@ -244,36 +243,36 @@ void histogram(const char *distfile, real binwidth, int n, int nset, real **val,
 {
     FILE          *fp;
     int            i, s;
-    double         min, max;
+    double         minval, maxval;
     int            nbin;
     gmx_int64_t   *histo;
 
-    min = val[0][0];
-    max = val[0][0];
+    minval = val[0][0];
+    maxval = val[0][0];
     for (s = 0; s < nset; s++)
     {
         for (i = 0; i < n; i++)
         {
-            if (val[s][i] < min)
+            if (val[s][i] < minval)
             {
-                min = val[s][i];
+                minval = val[s][i];
             }
-            else if (val[s][i] > max)
+            else if (val[s][i] > maxval)
             {
-                max = val[s][i];
+                maxval = val[s][i];
             }
         }
     }
 
-    min = binwidth*floor(min/binwidth);
-    max = binwidth*ceil(max/binwidth);
-    if (min != 0)
+    minval = binwidth*std::floor(minval/binwidth);
+    maxval = binwidth*std::ceil(maxval/binwidth);
+    if (minval != 0)
     {
-        min -= binwidth;
+        minval -= binwidth;
     }
-    max += binwidth;
+    maxval += binwidth;
 
-    nbin = (int)((max - min)/binwidth + 0.5) + 1;
+    nbin = static_cast<int>(((maxval - minval)/binwidth + 0.5) + 1);
     fprintf(stderr, "Making distributions with %d bins\n", nbin);
     snew(histo, nbin);
     fp = xvgropen(distfile, "Distribution", "", "", oenv);
@@ -285,11 +284,11 @@ void histogram(const char *distfile, real binwidth, int n, int nset, real **val,
         }
         for (i = 0; i < n; i++)
         {
-            histo[(int)((val[s][i] - min)/binwidth + 0.5)]++;
+            histo[static_cast<int>((val[s][i] - minval)/binwidth + 0.5)]++;
         }
         for (i = 0; i < nbin; i++)
         {
-            fprintf(fp, " %g  %g\n", min+i*binwidth, (double)histo[i]/(n*binwidth));
+            fprintf(fp, " %g  %g\n", minval+i*binwidth, static_cast<double>(histo[i])/(n*binwidth));
         }
         if (s < nset-1)
         {
@@ -336,9 +335,9 @@ static void average(const char *avfile, int avbar_opt,
         {
             snew(tmp, nset);
             fprintf(fp, "@TYPE xydydy\n");
-            edge = (int)(nset*0.05+0.5);
+            edge = static_cast<int>(nset*0.05+0.5);
             fprintf(stdout, "Errorbars: discarding %d points on both sides: %d%%"
-                    " interval\n", edge, (int)(100*(nset-2*edge)/nset+0.5));
+                    " interval\n", edge, static_cast<int>(100*(nset-2*edge)/nset+0.5));
         }
         else
         {
@@ -375,11 +374,11 @@ static void average(const char *avfile, int avbar_opt,
                 }
                 if (avbar_opt == avbarSTDDEV)
                 {
-                    err = sqrt(var/nset);
+                    err = std::sqrt(var/nset);
                 }
                 else
                 {
-                    err = sqrt(var/(nset*(nset-1)));
+                    err = std::sqrt(var/(nset*(nset-1)));
                 }
                 fprintf(fp, " %g", err);
             }
@@ -408,7 +407,7 @@ static real optimal_error_estimate(double sigma, double fitparm[], real tTotal)
         return 0;
     }
 
-    return sigma*sqrt(2*ss/tTotal);
+    return sigma*std::sqrt(2*ss/tTotal);
 }
 
 static void estimate_error(const char *eefile, int nb_min, int resol, int n,
@@ -445,7 +444,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n,
     xvgr_legend(fp, 2*nset, (const char**)leg, oenv);
     sfree(leg);
 
-    spacing = pow(2, 1.0/resol);
+    spacing = std::pow(2.0, 1.0/resol);
     snew(tbs, n);
     snew(ybs, n);
     snew(fitsig, n);
@@ -456,7 +455,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n,
         nbr     = nb_min;
         while (nbr <= n)
         {
-            bs = n/(int)nbr;
+            bs = n/static_cast<int>(nbr);
             if (bs != prev_bs)
             {
                 nb  = n/bs;
@@ -482,7 +481,6 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n,
                 nbs++;
             }
             nbr    *= spacing;
-            nb      = (int)(nbr+0.5);
             prev_bs = bs;
         }
         if (sig[s] == 0)
@@ -507,7 +505,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n,
              * For a single exponential autocorrelation: ybs(tau1) = 2/e tau1
              * From this we take our initial guess for tau1.
              */
-            twooe = 2/exp(1);
+            twooe = 2.0/std::exp(1.0);
             i     = -1;
             do
             {
@@ -553,7 +551,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n,
                 {
                     dens = 0.5*(tbs[i+1]/tbs[i-1] - 1);
                 }
-                fitsig[i] = sqrt((tau_sig + tbs[i])/dens);
+                fitsig[i] = std::sqrt((tau_sig + tbs[i])/dens);
             }
 
             if (!bSingleExpFit)
@@ -563,7 +561,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n,
                 /* We set the initial guess for tau2
                  * to halfway between tau1_est and the total time (on log scale).
                  */
-                fitparm[2] = sqrt(tau1_est*(n-1)*dt);
+                fitparm[2] = std::sqrt(tau1_est*(n-1)*dt);
                 do_lmfit(nbs, ybs, fitsig, 0, tbs, 0, dt*n, oenv,
                          bDebugMode(), effnERREST, fitparm, 0,
                          NULL);
@@ -637,8 +635,8 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n,
         }
         for (i = 0; i < nbs; i++)
         {
-            fprintf(fp, "%g %g %g\n", tbs[i], sig[s]*sqrt(ybs[i]/(n*dt)),
-                    sig[s]*sqrt(fit_function(effnERREST, fitparm, tbs[i])/(n*dt)));
+            fprintf(fp, "%g %g %g\n", tbs[i], sig[s]*std::sqrt(ybs[i]/(n*dt)),
+                    sig[s]*std::sqrt(fit_function(effnERREST, fitparm, tbs[i])/(n*dt)));
         }
 
         if (bFitAc)
@@ -653,7 +651,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n,
                 ac[i] = val[s][i] - av[s];
                 if (i > 0)
                 {
-                    fitsig[i] = sqrt(i);
+                    fitsig[i] = std::sqrt(static_cast<real>(i));
                 }
                 else
                 {
@@ -677,7 +675,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n,
             /* Generate more or less appropriate sigma's */
             for (i = 0; i <= fitlen; i++)
             {
-                fitsig[i] = sqrt(acint + dt*i);
+                fitsig[i] = std::sqrt(acint + dt*i);
             }
 
             ac_fit[0] = 0.5*acint;
@@ -695,7 +693,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n,
             for (i = 0; i < nbs; i++)
             {
                 fprintf(fp, "%g %g\n", tbs[i],
-                        sig[s]*sqrt(fit_function(effnERREST, ac_fit, tbs[i]))/(n*dt));
+                        sig[s]*std::sqrt(fit_function(effnERREST, ac_fit, tbs[i]))/(n*dt));
             }
 
             sfree(ac);
@@ -714,9 +712,8 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n,
 static void luzar_correl(int nn, real *time, int nset, real **val, real temp,
                          gmx_bool bError, real fit_start)
 {
-    const real tol = 1e-8;
     real      *kt;
-    real       weight, d2;
+    real       d2;
     int        j;
 
     please_cite(stdout, "Spoel2006b");
@@ -737,7 +734,7 @@ static void luzar_correl(int nn, real *time, int nset, real **val, real temp,
             {
                 d2 += sqr(kt[j] - val[3][j]);
             }
-            fprintf(debug, "RMS difference in derivatives is %g\n", sqrt(d2/nn));
+            fprintf(debug, "RMS difference in derivatives is %g\n", std::sqrt(d2/nn));
         }
         analyse_corr(nn, time, val[0], val[2], kt, NULL, NULL, NULL, fit_start,
                      temp);
@@ -760,13 +757,13 @@ static void filter(real flen, int n, int nset, real **val, real dt)
     int     f, s, i, j;
     double *filt, sum, vf, fluc, fluctot;
 
-    f = (int)(flen/(2*dt));
+    f = static_cast<int>(flen/(2*dt));
     snew(filt, f+1);
     filt[0] = 1;
     sum     = 1;
     for (i = 1; i <= f; i++)
     {
-        filt[i] = cos(M_PI*dt*i/flen);
+        filt[i] = std::cos(M_PI*dt*i/flen);
         sum    += 2*filt[i];
     }
     for (i = 0; i <= f; i++)
@@ -790,9 +787,9 @@ static void filter(real flen, int n, int nset, real **val, real dt)
         }
         fluc    /= n - 2*f;
         fluctot += fluc;
-        fprintf(stdout, "Set %3d filtered fluctuation: %12.6e\n", s+1, sqrt(fluc));
+        fprintf(stdout, "Set %3d filtered fluctuation: %12.6e\n", s+1, std::sqrt(fluc));
     }
-    fprintf(stdout, "Overall filtered fluctuation: %12.6e\n", sqrt(fluctot/nset));
+    fprintf(stdout, "Overall filtered fluctuation: %12.6e\n", std::sqrt(fluctot/nset));
     fprintf(stdout, "\n");
 
     sfree(filt);
@@ -926,10 +923,10 @@ static void print_fitted_function(const char   *fitfile,
         int   s, buflen = 0;
         if (NULL != fn_fitted)
         {
-            buflen = strlen(fn_fitted)+32;
+            buflen = std::strlen(fn_fitted)+32;
             snew(buf2, buflen);
-            strncpy(buf2, fn_fitted, buflen);
-            buf2[strlen(buf2)-4] = '\0';
+            std::strncpy(buf2, fn_fitted, buflen);
+            buf2[std::strlen(buf2)-4] = '\0';
         }
         for (s = 0; s < nset; s++)
         {
@@ -1039,9 +1036,9 @@ int gmx_analyze(int argc, char *argv[])
     static real        tb         = -1, te = -1, frac = 0.5, filtlen = 0, binwidth = 0.1, aver_start = 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    bIntegrate = FALSE, bRegression = FALSE, bLuzar = FALSE, bLuzarError = FALSE;
-    static int         nsets_in   = 1, d = 1, nb_min = 4, resol = 10, nFitPoints = 100;
-    static real        temp       = 298.15, fit_start = 1, fit_end = 60, diffusion = 5e-5, rcut = 0.35;
+    static gmx_bool    bIntegrate = FALSE, bRegression = FALSE, bLuzar = FALSE;
+    static int         nsets_in   = 1, d = 1, nb_min = 4, resol = 10;
+    static real        temp       = 298.15, fit_start = 1, fit_end = 60;
 
     /* must correspond to enum avbar* declared at beginning of file */
     static const char *avbar_opt[avbarNR+1] = {
@@ -1106,7 +1103,7 @@ int gmx_analyze(int argc, char *argv[])
     };
 #define NPA asize(pa)
 
-    FILE           *out, *out_fit;
+    FILE           *out;
     int             n, nlast, s, nset, i, j = 0;
     real          **val, *t, dt, tot, error;
     double         *av, *sig, cum1, cum2, cum3, cum4, db;
@@ -1232,10 +1229,10 @@ int gmx_analyze(int argc, char *argv[])
         cum3  /= n;
         cum4  /= n;
         av[s]  = cum1;
-        sig[s] = sqrt(cum2);
+        sig[s] = std::sqrt(cum2);
         if (n > 1)
         {
-            error = sqrt(cum2/(n-1));
+            error = std::sqrt(cum2/(n-1));
         }
         else
         {
@@ -1243,7 +1240,7 @@ int gmx_analyze(int argc, char *argv[])
         }
         printf("SS%d  %13.6e   %12.6e   %12.6e      %6.3f   %6.3f\n",
                s+1, av[s], sig[s], error,
-               sig[s] ? cum3/(sig[s]*sig[s]*sig[s]*sqrt(8/M_PI)) : 0,
+               sig[s] ? cum3/(sig[s]*sig[s]*sig[s]*std::sqrt(8/M_PI)) : 0,
                sig[s] ? cum4/(sig[s]*sig[s]*sig[s]*sig[s]*3)-1 : 0);
     }
     printf("\n");
@@ -1257,7 +1254,7 @@ int gmx_analyze(int argc, char *argv[])
     {
         out = xvgropen(msdfile, "Mean square displacement",
                        "time", "MSD (nm\\S2\\N)", oenv);
-        nlast = (int)(n*frac);
+        nlast = static_cast<int>(n*frac);
         for (s = 0; s < nset; s++)
         {
             for (j = 0; j <= nlast; j++)
@@ -1271,7 +1268,7 @@ int gmx_analyze(int argc, char *argv[])
                 {
                     tot += sqr(val[s][i]-val[s][i+j]);
                 }
-                tot /= (real)(n-j);
+                tot /= (n-j);
                 fprintf(out, " %g %8g\n", dt*j, tot);
             }
             if (s < nset-1)