Moved first batch of analysis tool source to C++
[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_pme_error.cpp
index fabca7c8d3cd497f6b728aabb869c917e2d74e26..d7589f8b4d288ad9b984920abb43a55dc63daba7 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
+#include "gmxpre.h"
+
+#include "config.h"
+
+#include <cmath>
+
 #include <algorithm>
 
 #include "gromacs/commandline/pargs.h"
-#include "typedefs.h"
-#include "types/commrec.h"
-#include "gromacs/utility/smalloc.h"
-#include "vec.h"
-#include "copyrite.h"
 #include "gromacs/fileio/tpxio.h"
-#include "readinp.h"
-#include "calcgrid.h"
-#include "checkpoint.h"
-#include "gmx_ana.h"
+#include "gromacs/gmxana/gmx_ana.h"
+#include "gromacs/legacyheaders/calcgrid.h"
+#include "gromacs/legacyheaders/checkpoint.h"
+#include "gromacs/legacyheaders/copyrite.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/main.h"
+#include "gromacs/legacyheaders/mdatoms.h"
+#include "gromacs/legacyheaders/network.h"
+#include "gromacs/legacyheaders/readinp.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/calculate-ewald-splitting-coefficient.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/vec.h"
 #include "gromacs/random/random.h"
-#include "physics.h"
-#include "mdatoms.h"
-#include "coulomb.h"
-#include "mtop_util.h"
-#include "network.h"
-#include "main.h"
-#include "macros.h"
-
-#include "gromacs/legacyheaders/gmx_fatal.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
 
-/* We use the same defines as in mvdata.c here */
+/* We use the same defines as in broadcaststructs.cpp here */
 #define  block_bc(cr,   d) gmx_bcast(     sizeof(d),     &(d), (cr))
 #define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
 #define   snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
@@ -208,19 +213,19 @@ static inline real eps_poly1(
     {
         tmp  = m / K + i;
         tmp *= 2.0*M_PI;
-        nom += pow( tmp, -n );
+        nom += std::pow( tmp, -n );
     }
 
     for (i = SUMORDER; i > 0; i--)
     {
         tmp  = m / K + i;
         tmp *= 2.0*M_PI;
-        nom += pow( tmp, -n );
+        nom += std::pow( tmp, -n );
     }
 
     tmp   = m / K;
     tmp  *= 2.0*M_PI;
-    denom = pow( tmp, -n )+nom;
+    denom = std::pow( tmp, -n )+nom;
 
     return -nom/denom;
 
@@ -245,21 +250,21 @@ static inline real eps_poly2(
     {
         tmp  = m / K + i;
         tmp *= 2.0*M_PI;
-        nom += pow( tmp, -2*n );
+        nom += std::pow( tmp, -2*n );
     }
 
     for (i = SUMORDER; i > 0; i--)
     {
         tmp  = m / K + i;
         tmp *= 2.0*M_PI;
-        nom += pow( tmp, -2*n );
+        nom += std::pow( tmp, -2*n );
     }
 
     for (i = -SUMORDER; i < SUMORDER+1; i++)
     {
         tmp    = m / K + i;
         tmp   *= 2.0*M_PI;
-        denom += pow( tmp, -n );
+        denom += std::pow( tmp, -n );
     }
     tmp = eps_poly1(m, K, n);
     return nom / denom / denom + tmp*tmp;
@@ -285,21 +290,21 @@ static inline real eps_poly3(
     {
         tmp  = m / K + i;
         tmp *= 2.0*M_PI;
-        nom += i * pow( tmp, -2*n );
+        nom += i * std::pow( tmp, -2*n );
     }
 
     for (i = SUMORDER; i > 0; i--)
     {
         tmp  = m / K + i;
         tmp *= 2.0*M_PI;
-        nom += i * pow( tmp, -2*n );
+        nom += i * std::pow( tmp, -2*n );
     }
 
     for (i = -SUMORDER; i < SUMORDER+1; i++)
     {
         tmp    = m / K + i;
         tmp   *= 2.0*M_PI;
-        denom += pow( tmp, -n );
+        denom += std::pow( tmp, -n );
     }
 
     return 2.0 * M_PI * nom / denom / denom;
@@ -325,21 +330,21 @@ static inline real eps_poly4(
     {
         tmp  = m / K + i;
         tmp *= 2.0*M_PI;
-        nom += i * i * pow( tmp, -2*n );
+        nom += i * i * std::pow( tmp, -2*n );
     }
 
     for (i = SUMORDER; i > 0; i--)
     {
         tmp  = m / K + i;
         tmp *= 2.0*M_PI;
-        nom += i * i * pow( tmp, -2*n );
+        nom += i * i * std::pow( tmp, -2*n );
     }
 
     for (i = -SUMORDER; i < SUMORDER+1; i++)
     {
         tmp    = m / K + i;
         tmp   *= 2.0*M_PI;
-        denom += pow( tmp, -n );
+        denom += std::pow( tmp, -n );
     }
 
     return 4.0 * M_PI * M_PI * nom / denom / denom;
@@ -372,18 +377,18 @@ static inline real eps_self(
 
     for (i = -SUMORDER; i < 0; i++)
     {
-        tmp    = -sin(2.0 * M_PI * i * K * rcoord);
+        tmp    = -std::sin(2.0 * M_PI * i * K * rcoord);
         tmp1   = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
-        tmp2   = pow(tmp1, -n);
+        tmp2   = std::pow(tmp1, -n);
         nom   += tmp * tmp2 * i;
         denom += tmp2;
     }
 
     for (i = SUMORDER; i > 0; i--)
     {
-        tmp    = -sin(2.0 * M_PI * i * K * rcoord);
+        tmp    = -std::sin(2.0 * M_PI * i * K * rcoord);
         tmp1   = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
-        tmp2   = pow(tmp1, -n);
+        tmp2   = std::pow(tmp1, -n);
         nom   += tmp * tmp2 * i;
         denom += tmp2;
     }
@@ -488,7 +493,7 @@ static real estimate_reciprocal(
     xtot        = stopglobal*2+1;
     if (PAR(cr))
     {
-        x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
+        x_per_core = static_cast<int>(std::ceil(static_cast<real>(xtot) / cr->nnodes));
         startlocal = startglobal + x_per_core*cr->nodeid;
         stoplocal  = startlocal + x_per_core -1;
         if (stoplocal > stopglobal)
@@ -540,7 +545,7 @@ static real estimate_reciprocal(
                 svmul(nz, info->recipbox[ZZ], tmpvec);
                 rvec_add(gridpxy, tmpvec, gridp);
                 tmp    = norm2(gridp);
-                coeff  = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
+                coeff  = std::exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
                 coeff /= 2.0 * M_PI * info->volume * tmp;
                 coeff2 = tmp;
 
@@ -634,14 +639,14 @@ static real estimate_reciprocal(
         /* Here xtot is the number of samples taken for the Monte Carlo calculation
          * of the average of term IV of equation 35 in Wang2010. Round up to a
          * number of samples that is divisible by the number of nodes */
-        x_per_core  = static_cast<int>(ceil(info->fracself * nr / cr->nnodes));
+        x_per_core  = static_cast<int>(std::ceil(info->fracself * nr / cr->nnodes));
         xtot        = x_per_core * cr->nnodes;
     }
     else
     {
         /* In this case we use all nr particle positions */
         xtot       = nr;
-        x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
+        x_per_core = static_cast<int>(std::ceil(static_cast<real>(xtot) / cr->nnodes));
     }
 
     startlocal = x_per_core *  cr->nodeid;
@@ -657,7 +662,7 @@ static real estimate_reciprocal(
         {
             for (i = 0; i < xtot; i++)
             {
-                numbers[i] = static_cast<int>(floor(gmx_rng_uniform_real(rng) * nr));
+                numbers[i] = static_cast<int>(std::floor(gmx_rng_uniform_real(rng) * nr));
             }
         }
         /* Broadcast the random number array to the other nodes */
@@ -717,7 +722,7 @@ static real estimate_reciprocal(
                     svmul(nz, info->recipbox[ZZ], tmpvec);
                     rvec_add(gridpxy, tmpvec, gridp);
                     tmp      = norm2(gridp);
-                    coeff    = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
+                    coeff    = std::exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
                     coeff   /= tmp;
                     e_rec3x += coeff*eps_self(nx, info->nkx[0], info->recipbox[XX], info->pme_order[0], x[ci]);
                     e_rec3y += coeff*eps_self(ny, info->nky[0], info->recipbox[YY], info->pme_order[0], x[ci]);
@@ -779,7 +784,7 @@ static real estimate_reciprocal(
        e_rec2*=  q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
        e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
      */
-    e_rec = sqrt(e_rec1+e_rec2+e_rec3);
+    e_rec = std::sqrt(e_rec1+e_rec2+e_rec3);
 
 
     return ONE_4PI_EPS0 * e_rec;
@@ -1020,7 +1025,7 @@ static void estimate_PME_error(t_inputinfo *info, t_state *state,
         edir = info->e_dir[0];
         erec = info->e_rec[0];
         derr = edir-erec;
-        while (fabs(derr/std::min(erec, edir)) > 1e-4)
+        while (std::abs(derr/std::min(erec, edir)) > 1e-4)
         {
 
             beta                = info->ewald_beta[0];
@@ -1045,7 +1050,7 @@ static void estimate_PME_error(t_inputinfo *info, t_state *state,
             if (MASTER(cr))
             {
                 i++;
-                fprintf(stderr, "difference between real and rec. space error (step %d): %g\n", i, fabs(derr));
+                fprintf(stderr, "difference between real and rec. space error (step %d): %g\n", i, std::abs(derr));
                 fprintf(stderr, "old beta: %f\n", beta0);
                 fprintf(stderr, "new beta: %f\n", beta);
             }
@@ -1102,16 +1107,16 @@ int gmx_pme_error(int argc, char *argv[])
 
 
     static t_filenm fnm[] = {
-        { efTPX, "-s",     NULL,    ffREAD },
+        { efTPR, "-s",     NULL,    ffREAD },
         { efOUT, "-o",    "error",  ffWRITE },
-        { efTPX, "-so",   "tuned",  ffOPTWR }
+        { efTPR, "-so",   "tuned",  ffOPTWR }
     };
 
     output_env_t    oenv = NULL;
 
     t_pargs         pa[] = {
         { "-beta",     FALSE, etREAL, {&user_beta},
-          "If positive, overwrite ewald_beta from [TT].tpr[tt] file with this value" },
+          "If positive, overwrite ewald_beta from [REF].tpr[ref] file with this value" },
         { "-tune",     FALSE, etBOOL, {&bTUNE},
           "Tune the splitting parameter such that the error is equally distributed between real and reciprocal space" },
         { "-self",     FALSE, etREAL, {&fracself},
@@ -1126,12 +1131,8 @@ int gmx_pme_error(int argc, char *argv[])
 #define NFILE asize(fnm)
 
     cr = init_commrec();
-#ifdef GMX_LIB_MPI
-    MPI_Barrier(MPI_COMM_WORLD);
-#endif
 
     PCA_Flags  = PCA_NOEXIT_ON_ARGS;
-    PCA_Flags |= (MASTER(cr) ? 0 : PCA_QUIET);
 
     if (!parse_common_args(&argc, argv, PCA_Flags,
                            NFILE, fnm, asize(pa), pa, asize(desc), desc,