Merge branch 'master' into pygromacs

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

diff --git a/src/gromacs/gmxana/gmx_dos.c b/src/gromacs/gmxana/gmx_dos.cpp
similarity index 92%
rename from src/gromacs/gmxana/gmx_dos.c
rename to src/gromacs/gmxana/gmx_dos.cpp
index d4423adc47ca1e50d943d761b0b704f1f4d2d1c6..04abdb911b6b00f1755ac0cd1f7270d5fac4b888 100644 (file)
--- a/src/gromacs/gmxana/gmx_dos.c
+++ b/src/gromacs/gmxana/gmx_dos.cpp
@@ -34,10 +34,10 @@
  */
 #include "gmxpre.h"
 
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/correlationfunctions/autocorr.h"
@@ -102,16 +102,16 @@ static int calcMoleculesInIndexGroup(t_block *mols, int natoms, atom_id *index,
 
 static double FD(double Delta, double f)
 {
-    return (2*pow(Delta, -4.5)*pow(f, 7.5) -
-            6*pow(Delta, -3)*pow(f, 5) -
-            pow(Delta, -1.5)*pow(f, 3.5) +
-            6*pow(Delta, -1.5)*pow(f, 2.5) +
+    return (2*std::pow(Delta, -4.5)*std::pow(f, 7.5) -
+            6*std::pow(Delta, -3.0)*std::pow(f, 5.0) -
+            std::pow(Delta, -1.5)*std::pow(f, 3.5) +
+            6*std::pow(Delta, -1.5)*std::pow(f, 2.5) +
             2*f - 2);
 }
 
 static double YYY(double f, double y)
 {
-    return (2*pow(y*f, 3) - sqr(f)*y*(1+6*y) +
+    return (2*std::pow(y*f, 3.0) - sqr(f)*y*(1+6*y) +
             (2+6*y)*f - 2);
 }
 
@@ -121,14 +121,14 @@ static double calc_compress(double y)
     {
         return 0;
     }
-    return ((1+y+sqr(y)-pow(y, 3))/(pow(1-y, 3)));
+    return ((1+y+sqr(y)-std::pow(y, 3.0))/(std::pow(1-y, 3.0)));
 }
 
 static double bisector(double Delta, double tol,
                        double ff0, double ff1,
                        double ff(double, double))
 {
-    double fd0, fd, fd1, f, f0, f1;
+    double fd, f, f0, f1;
     double tolmin = 1e-8;
 
     f0 = ff0;
@@ -141,8 +141,6 @@ static double bisector(double Delta, double tol,
 
     do
     {
-        fd0 = ff(Delta, f0);
-        fd1 = ff(Delta, f1);
         f   = (f0+f1)*0.5;
         fd  = ff(Delta, f);
         if (fd < 0)
@@ -172,9 +170,9 @@ static double calc_y(double f, double Delta, double toler)
 {
     double y1, y2;
 
-    y1 = pow(f/Delta, 1.5);
+    y1 = std::pow(f/Delta, 1.5);
     y2 = bisector(f, toler, 0, 10000, YYY);
-    if (fabs((y1-y2)/(y1+y2)) > 100*toler)
+    if (std::abs((y1-y2)/(y1+y2)) > 100*toler)
     {
         fprintf(stderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n",
                 y1, y2);
@@ -187,7 +185,7 @@ static double calc_Shs(double f, double y)
 {
     double fy  = f*y;
 
-    return BOLTZ*(log(calc_compress(fy)) + fy*(3*fy-4)/sqr(1-fy));
+    return BOLTZ*(std::log(calc_compress(fy)) + fy*(3*fy-4)/sqr(1-fy));
 }
 
 static real wCsolid(real nu, real beta)
@@ -201,7 +199,7 @@ static real wCsolid(real nu, real beta)
     }
     else
     {
-        ebn  = exp(bhn);
+        ebn  = std::exp(bhn);
         koko = sqr(1-ebn);
         return sqr(bhn)*ebn/koko;
     }
@@ -217,7 +215,7 @@ static real wSsolid(real nu, real beta)
     }
     else
     {
-        return bhn/gmx_expm1(bhn) - gmx_log1p(-exp(-bhn));
+        return bhn/gmx_expm1(bhn) - gmx_log1p(-std::exp(-bhn));
     }
 }
 
@@ -231,7 +229,7 @@ static real wAsolid(real nu, real beta)
     }
     else
     {
-        return log((1-exp(-bhn))/(exp(-bhn/2))) - log(bhn);
+        return std::log((1-std::exp(-bhn))/(std::exp(-bhn/2))) - std::log(bhn);
     }
 }
 
@@ -274,14 +272,14 @@ int gmx_dos(int argc, char *argv[])
     matrix              box;
     int                 gnx;
     char                title[256];
-    real                t0, t1, m;
+    real                t0, t1;
     t_trxstatus        *status;
-    int                 nV, nframes, n_alloc, i, j, k, l, fftcode, Nmol, Natom;
-    double              rho, dt, V2sum, Vsum, V, tmass, dostot, dos2, dosabs;
+    int                 nV, nframes, n_alloc, i, j, fftcode, Nmol, Natom;
+    double              rho, dt, Vsum, V, tmass, dostot, dos2;
     real              **c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
     output_env_t        oenv;
     gmx_fft_t           fft;
-    double              cP, S, A, E, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
+    double              cP, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
     double              wCdiff, wSdiff, wAdiff, wEdiff;
     int                 grpNatoms;
     atom_id            *index;
@@ -368,14 +366,13 @@ int gmx_dos(int argc, char *argv[])
 
     n_alloc = 0;
     nframes = 0;
-    Vsum    = V2sum = 0;
+    Vsum    = 0;
     nV      = 0;
     do
     {
         if (fr.bBox)
         {
             V      = det(fr.box);
-            V2sum += V*V;
             Vsum  += V;
             nV++;
         }
@@ -496,7 +493,7 @@ int gmx_dos(int argc, char *argv[])
         dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
         if (bAbsolute)
         {
-            dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
+            dos[DOS][j] = bfac*std::sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
         }
         else
         {
@@ -517,15 +514,15 @@ int gmx_dos(int argc, char *argv[])
     DoS0 = dos[DOS][0];
 
     /* Note this eqn. is incorrect in Pascal2011a! */
-    Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
-             pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));
+    Delta = ((2*DoS0/(9*Natom))*std::sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
+             std::pow((Natom/V), 1.0/3.0)*std::pow(6.0/M_PI, 2.0/3.0));
     f     = calc_fluidicity(Delta, toler);
     y     = calc_y(f, Delta, toler);
     z     = calc_compress(y);
-    Sig   = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
+    Sig   = BOLTZ*(5.0/2.0+std::log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
     Shs   = Sig+calc_Shs(f, y);
     rho   = (tmass*AMU)/(V*NANO*NANO*NANO);
-    sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);
+    sigHS = std::pow(6*y*V/(M_PI*Natom), 1.0/3.0);
 
     fprintf(fplog, "System = \"%s\"\n", title);
     fprintf(fplog, "Nmol = %d\n", Nmol);