Moved second half of gmxana tools to C++
[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_msd.cpp
similarity index 97%
rename from src/gromacs/gmxana/gmx_msd.c
rename to src/gromacs/gmxana/gmx_msd.cpp
index 3a6361d8969d93c6a82784dfbd9462c91f6cda5d..90fe9725e2a78659d7fc99b5274633e1234d3e40 100644 (file)
@@ -36,8 +36,8 @@
  */
 #include "gmxpre.h"
 
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/fileio/confio.h"
@@ -55,6 +55,7 @@
 #include "gromacs/topology/index.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 
 #define FACTOR  1000.0  /* Convert nm^2/ps to 10e-5 cm^2/s */
@@ -229,7 +230,7 @@ static void calc_corr(t_corr *curr, int nr, int nx, atom_id index[], rvec xc[],
     {
         if ((thistime(curr) >= (curr->nlast*curr->delta_t)) && (nr == 0))
         {
-            memcpy(curr->x0[curr->nlast], xc, curr->ncoords*sizeof(xc[0]));
+            std::memcpy(curr->x0[curr->nlast], xc, curr->ncoords*sizeof(xc[0]));
             curr->n_offs[curr->nlast] = curr->nframes;
             copy_rvec(com, curr->com[curr->nlast]);
             curr->nlast++;
@@ -490,7 +491,6 @@ static void calc_com(gmx_bool bMol, int gnx, atom_id index[],
 
     clear_dvec(sx);
     tmass = 0;
-    mass  = 1;
 
     prep_data(bMol, gnx, index, xcur, xprev, box);
     for (i = 0; (i < gnx); i++)
@@ -595,7 +595,7 @@ void printmol(t_corr *curr, const char *fn,
         fprintf(out, "%10d  %10g\n", i, D);
         if (pdbinfo)
         {
-            sqrtD = sqrt(D);
+            sqrtD = std::sqrt(D);
             if (sqrtD > sqrtD_max)
             {
                 sqrtD_max = sqrtD;
@@ -614,7 +614,7 @@ void printmol(t_corr *curr, const char *fn,
     D2av /= curr->nmol;
     VarD  = D2av - sqr(Dav);
     printf("<D> = %.4f Std. Dev. = %.4f Error = %.4f\n",
-           Dav, sqrt(VarD), sqrt(VarD/curr->nmol));
+           Dav, std::sqrt(VarD), std::sqrt(VarD/curr->nmol));
 
     if (fn_pdb && x)
     {
@@ -627,6 +627,7 @@ void printmol(t_corr *curr, const char *fn,
         {
             scale *= 10;
         }
+        GMX_RELEASE_ASSERT(pdbinfo != NULL, "Internal error - pdbinfo not set for PDB input");
         for (i = 0; i < top->atoms.nr; i++)
         {
             pdbinfo[i].bfac *= scale;
@@ -778,7 +779,7 @@ int corr_loop(t_corr *curr, const char *fn, t_topology *top, int ePBC,
         /* for the first frame, the previous frame is a copy of the first frame */
         if (bFirst)
         {
-            memcpy(xa[prev], xa[cur], curr->ncoords*sizeof(xa[prev][0]));
+            std::memcpy(xa[prev], xa[cur], curr->ncoords*sizeof(xa[prev][0]));
             bFirst = FALSE;
         }
 
@@ -953,7 +954,7 @@ void do_corr(const char *trx_file, const char *ndx_file, const char *msd_file,
 
     if (beginfit == -1)
     {
-        i0       = (int)(0.1*(msd->nframes - 1) + 0.5);
+        i0       = static_cast<int>(0.1*(msd->nframes - 1) + 0.5);
         beginfit = msd->time[i0];
     }
     else
@@ -966,7 +967,7 @@ void do_corr(const char *trx_file, const char *ndx_file, const char *msd_file,
 
     if (endfit == -1)
     {
-        i1     = (int)(0.9*(msd->nframes - 1) + 0.5) + 1;
+        i1     = static_cast<int>(0.9*(msd->nframes - 1) + 0.5) + 1;
         endfit = msd->time[i1-1];
     }
     else
@@ -995,7 +996,7 @@ void do_corr(const char *trx_file, const char *ndx_file, const char *msd_file,
             {
                 lsq_y_ax_b(N/2, &(msd->time[i0]), &(msd->data[j][i0]), &a, &b, &r, &chi2);
                 lsq_y_ax_b(N/2, &(msd->time[i0+N/2]), &(msd->data[j][i0+N/2]), &a2, &b, &r, &chi2);
-                SigmaD[j] = fabs(a-a2);
+                SigmaD[j] = std::abs(a-a2);
             }
             else
             {
@@ -1160,6 +1161,9 @@ int gmx_msd(int argc, char *argv[])
         fprintf(stderr, "Calculating diffusion coefficients for molecules.\n");
     }
 
+    GMX_RELEASE_ASSERT(normtype[0] != 0, "Options inconsistency; normtype[0] is NULL");
+    GMX_RELEASE_ASSERT(axtitle[0] != 0, "Options inconsistency; axtitle[0] is NULL");
+
     if (normtype[0][0] != 'n')
     {
         type       = normtype[0][0] - 'x' + X; /* See defines above */