Moved first batch of analysis tool source to C++

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

diff --git a/src/gromacs/gmxana/gmx_current.c b/src/gromacs/gmxana/gmx_current.cpp
similarity index 94%
rename from src/gromacs/gmxana/gmx_current.c
rename to src/gromacs/gmxana/gmx_current.cpp
index 107d87223055f1bb335e96f62d8c06d7623d18bc..833e5a5ef306b27f89c145ac744b531a9532cf76 100644 (file)
--- a/src/gromacs/gmxana/gmx_current.c
+++ b/src/gromacs/gmxana/gmx_current.cpp
@@ -34,8 +34,8 @@
  */
 #include "gmxpre.h"
 
-#include <assert.h>
-#include <stdlib.h>
+#include <cmath>
+#include <cstdlib>
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/fileio/confio.h"
@@ -52,9 +52,9 @@
 #include "gromacs/topology/index.h"
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 
-#define SQR(x) (pow(x, 2.0))
 #define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
 
 static void index_atom2mol(int *n, int *index, t_block *mols)
@@ -139,7 +139,7 @@ static gmx_bool precalc(t_topology top, real mass2[], real qmol[])
         }
     }
 
-    if (fabs(qall) > 0.01)
+    if (std::abs(qall) > 0.01)
     {
         printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
         bNEU = FALSE;
@@ -275,12 +275,10 @@ static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int
 {
 
     int  i;
-    real corint;
     real deltat = 0.0;
     real rfr;
     real sigma     = 0.0;
     real sigma_ret = 0.0;
-    corint = 0.0;
 
     if (nfr > 1)
     {
@@ -288,8 +286,9 @@ static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int
 
         while (i < nfr)
         {
+            real corint;
 
-            rfr      = (real) (nfr/nshift-i/nshift);
+            rfr      = static_cast<real>(nfr+i)/nshift;
             cacf[i] /= rfr;
 
             if (time[vfr[i]] <= time[vfr[ei]])
@@ -327,14 +326,9 @@ static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], i
 {
 
     int     i;
-    real    rtmp;
-    real    rfr;
-
 
     fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
 
-
-
     for (i = 0; i < nfr; i++)
     {
 
@@ -353,18 +347,17 @@ static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], i
 static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
                        FILE *fmjdsp, gmx_bool bNoJump, gmx_bool bACF, gmx_bool bINT,
                        int ePBC, t_topology top, t_trxframe fr, real temp,
-                       real trust, real bfit, real efit, real bvit, real evit,
+                       real bfit, real efit, real bvit, real evit,
                        t_trxstatus *status, int isize, int nmols, int nshift,
                        atom_id *index0, int indexm[], real mass2[],
                        real qmol[], real eps_rf, const output_env_t oenv)
 {
-    int       i, j, k, l, f;
-    int       valloc, nalloc, nfr, nvfr, m, itrust = 0;
+    int       i, j;
+    int       valloc, nalloc, nfr, nvfr;
     int       vshfr;
     real     *xshfr       = NULL;
     int      *vfr         = NULL;
     real      refr        = 0.0;
-    real      rfr         = 0.0;
     real     *cacf        = NULL;
     real     *time        = NULL;
     real     *djc         = NULL;
@@ -374,7 +367,6 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
     real      volume;
     real      volume_av = 0.0;
     real      dk_s, dk_d, dk_f;
-    real      dm_s, dm_d;
     real      mj    = 0.0;
     real      mj2   = 0.0;
     real      mjd   = 0.0;
@@ -392,12 +384,9 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
     real     *dsp2  = NULL;
     real      t0;
     real      rtmp;
-    real      qtmp;
-
 
-
-    rvec   tmp;
-    rvec  *mtrans = NULL;
+    rvec      tmp;
+    rvec     *mtrans = NULL;
 
     /*
      * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
@@ -449,7 +438,7 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
     do
     {
 
-        refr = (real)(nfr+1);
+        refr = (nfr+1);
 
         if (nfr >= nalloc)
         {
@@ -470,8 +459,7 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
                 xshfr[i] = 0.0;
             }
         }
-        assert(time != NULL);
-
+        GMX_RELEASE_ASSERT(time != NULL, "Memory not allocated correctly - time array is NULL");
 
         if (nfr == 0)
         {
@@ -652,9 +640,6 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
 
     if (v0 != NULL)
     {
-        trust *= (real) nvfr;
-        itrust = (int) trust;
-
         if (bINT)
         {
 
@@ -678,15 +663,15 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
                 for (i = ii; i <= ie; i++)
                 {
 
-                    xfit[i-ii] = log(time[vfr[i]]);
-                    rtmp       = fabs(cacf[i]);
-                    yfit[i-ii] = log(rtmp);
+                    xfit[i-ii] = std::log(time[vfr[i]]);
+                    rtmp       = std::abs(cacf[i]);
+                    yfit[i-ii] = std::log(rtmp);
 
                 }
 
                 lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
 
-                malt = exp(malt);
+                malt = std::exp(malt);
 
                 sigma += 1.0;
 
@@ -729,7 +714,7 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
     if (bACF && (ii < nvfr))
     {
         fprintf(stderr, "Integral and integrated fit to the current acf yields at t=%f:\n", time[vfr[ii]]);
-        fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk-malt*pow(time[vfr[ii]], sigma), sgk);
+        fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk-malt*std::pow(time[vfr[ii]], sigma), sgk);
     }
 
     if (ei > bi)
@@ -792,7 +777,6 @@ int gmx_current(int argc, char *argv[])
 
     static int             nshift  = 1000;
     static real            temp    = 300.0;
-    static real            trust   = 0.25;
     static real            eps_rf  = 0.0;
     static gmx_bool        bNoJump = TRUE;
     static real            bfit    = 100.0;
@@ -814,8 +798,6 @@ int gmx_current(int argc, char *argv[])
           "Begin of the fit of the current autocorrelation function to a*t^b."},
         { "-evit", FALSE, etREAL, {&evit},
           "End of the fit of the current autocorrelation function to a*t^b."},
-        { "-tr", FALSE, etREAL, {&trust},
-          "Fraction of the trajectory taken into account for the integral."},
         { "-temp", FALSE, etREAL, {&temp},
           "Temperature for calculating epsilon."}
     };
@@ -834,21 +816,13 @@ int gmx_current(int argc, char *argv[])
     int                    isize;
     t_trxstatus           *status;
     int                    flags = 0;
-    gmx_bool               bTop;
-    gmx_bool               bNEU;
     gmx_bool               bACF;
     gmx_bool               bINT;
     int                    ePBC = -1;
-    int                    natoms;
     int                    nmols;
-    int                    i, j, k = 0, l;
-    int                    step;
-    real                   t;
-    real                   lambda;
+    int                    i;
     real                  *qmol;
     FILE                  *outf   = NULL;
-    FILE                  *outi   = NULL;
-    FILE                  *tfile  = NULL;
     FILE                  *mcor   = NULL;
     FILE                  *fmj    = NULL;
     FILE                  *fmd    = NULL;
@@ -915,7 +889,7 @@ int gmx_current(int argc, char *argv[])
     bACF = opt2bSet("-caf", NFILE, fnm);
     bINT = opt2bSet("-mc", NFILE, fnm);
 
-    bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, &vtop, box, TRUE);
+    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, &vtop, box, TRUE);
 
 
 
@@ -933,7 +907,7 @@ int gmx_current(int argc, char *argv[])
     snew(mass2, top.atoms.nr);
     snew(qmol, top.atoms.nr);
 
-    bNEU = precalc(top, mass2, qmol);
+    precalc(top, mass2, qmol);
 
 
     snew(indexm, isize);
@@ -990,7 +964,7 @@ int gmx_current(int argc, char *argv[])
      * and calculates the requested quantities */
 
     dielectric(fmj, fmd, outf, fcur, mcor, fmjdsp, bNoJump, bACF, bINT, ePBC, top, fr,
-               temp, trust, bfit, efit, bvit, evit, status, isize, nmols, nshift,
+               temp, bfit, efit, bvit, evit, status, isize, nmols, nshift,
                index0, indexm, mass2, qmol, eps_rf, oenv);
 
     xvgrclose(fmj);