Merge branch release-5-1
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / grompp.cpp
similarity index 97%
rename from src/gromacs/gmxpreprocess/grompp.c
rename to src/gromacs/gmxpreprocess/grompp.cpp
index 2f976b8d2eef1c8667d524a3ac287ca999f9aef2..d8a6737091a6bc1bbee23613862035b37f925ce3 100644 (file)
 
 #include "grompp.h"
 
-#include <assert.h>
 #include <errno.h>
 #include <limits.h>
-#include <math.h>
 #include <string.h>
 
+#include <cmath>
+
+#include <algorithm>
+
 #include <sys/types.h>
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/fileio/confio.h"
 #include "gromacs/fileio/enxio.h"
-#include "gromacs/fileio/gmxfio.h"
 #include "gromacs/fileio/tpxio.h"
-#include "gromacs/fileio/trnio.h"
 #include "gromacs/fileio/trxio.h"
 #include "gromacs/gmxpreprocess/add_par.h"
 #include "gromacs/gmxpreprocess/convparm.h"
@@ -84,6 +84,7 @@
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/snprintf.h"
 
@@ -188,7 +189,7 @@ static void check_cg_sizes(const char *topfn, t_block *cgs, warninp_t wi)
     maxsize = 0;
     for (cg = 0; cg < cgs->nr; cg++)
     {
-        maxsize = max(maxsize, cgs->index[cg+1]-cgs->index[cg]);
+        maxsize = std::max(maxsize, cgs->index[cg+1]-cgs->index[cg]);
     }
 
     if (maxsize > MAX_CHARGEGROUP_SIZE)
@@ -290,7 +291,7 @@ static void check_bonds_timestep(gmx_mtop_t *mtop, double dt, warninp_t wi)
                 if (debug)
                 {
                     fprintf(debug, "fc %g m1 %g m2 %g period %g\n",
-                            fc, m1, m2, sqrt(period2));
+                            fc, m1, m2, std::sqrt(period2));
                 }
                 if (period2 < limit2)
                 {
@@ -335,7 +336,7 @@ static void check_bonds_timestep(gmx_mtop_t *mtop, double dt, warninp_t wi)
                 *w_moltype->name,
                 w_a1+1, *w_moltype->atoms.atomname[w_a1],
                 w_a2+1, *w_moltype->atoms.atomname[w_a2],
-                sqrt(w_period2), bWarn ? min_steps_warn : min_steps_note, dt,
+                std::sqrt(w_period2), bWarn ? min_steps_warn : min_steps_note, dt,
                 bWater ?
                 "Maybe you asked for fexible water." :
                 "Maybe you forgot to change the constraints mdp option.");
@@ -811,7 +812,7 @@ static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB,
                         rvec com,
                         warninp_t wi)
 {
-    gmx_bool        bFirst = TRUE, *hadAtom;
+    gmx_bool       *hadAtom;
     rvec           *x, *v, *xp;
     dvec            sum;
     double          totmass;
@@ -827,7 +828,7 @@ static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB,
     get_stx_coordnum(fn, &natoms);
     if (natoms != mtop->natoms)
     {
-        sprintf(warn_buf, "The number of atoms in %s (%d) does not match the number of atoms in the topology (%d). Will assume that the first %d atoms in the topology and %s match.", fn, natoms, mtop->natoms, min(mtop->natoms, natoms), fn);
+        sprintf(warn_buf, "The number of atoms in %s (%d) does not match the number of atoms in the topology (%d). Will assume that the first %d atoms in the topology and %s match.", fn, natoms, mtop->natoms, std::min(mtop->natoms, natoms), fn);
         warning(wi, warn_buf);
     }
     snew(x, natoms);
@@ -936,7 +937,7 @@ static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB,
 
     if (rc_scaling != erscNO)
     {
-        assert(npbcdim <= DIM);
+        GMX_ASSERT(npbcdim <= DIM, "Only DIM dimensions can have PBC");
 
         for (mb = 0; mb < mtop->nmolblock; mb++)
         {
@@ -994,8 +995,6 @@ static void gen_posres(gmx_mtop_t *mtop, t_molinfo *mi,
                        rvec com, rvec comB,
                        warninp_t wi)
 {
-    int i, j;
-
     read_posres  (mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi);
     /* It is safer to simply read the b-state posres rather than trying
      * to be smart and copy the positions.
@@ -1092,7 +1091,7 @@ interpolate1d( double     xmin,
     int    ix;
     double a, b;
 
-    ix = (x-xmin)/dx;
+    ix = static_cast<int>((x-xmin)/dx);
 
     a = (xmin+(ix+1)*dx-x)/dx;
     b = (x-xmin-ix*dx)/dx;
@@ -1182,7 +1181,7 @@ setup_cmap (int              grid_spacing,
 
 void init_cmap_grid(gmx_cmap_t *cmap_grid, int ngrid, int grid_spacing)
 {
-    int i, k, nelem;
+    int i, nelem;
 
     cmap_grid->ngrid        = ngrid;
     cmap_grid->grid_spacing = grid_spacing;
@@ -1309,22 +1308,19 @@ static real calc_temp(const gmx_mtop_t *mtop,
                       const t_inputrec *ir,
                       rvec             *v)
 {
-    double                  sum_mv2;
     gmx_mtop_atomloop_all_t aloop;
     t_atom                 *atom;
     int                     a;
-    int                     nrdf, g;
-
-    sum_mv2 = 0;
 
+    double                  sum_mv2 = 0;
     aloop = gmx_mtop_atomloop_all_init(mtop);
     while (gmx_mtop_atomloop_all_next(aloop, &a, &atom))
     {
         sum_mv2 += atom->m*norm2(v[a]);
     }
 
-    nrdf = 0;
-    for (g = 0; g < ir->opts.ngtc; g++)
+    double nrdf = 0;
+    for (int g = 0; g < ir->opts.ngtc; g++)
     {
         nrdf += ir->opts.nrdf[g];
     }
@@ -1349,7 +1345,7 @@ static real get_max_reference_temp(const t_inputrec *ir,
         }
         else
         {
-            ref_t = max(ref_t, ir->opts.ref_t[i]);
+            ref_t = std::max(ref_t, ir->opts.ref_t[i]);
         }
     }
 
@@ -1371,7 +1367,6 @@ static void set_verlet_buffer(const gmx_mtop_t *mtop,
                               matrix            box,
                               warninp_t         wi)
 {
-    int                    i;
     verletbuf_list_setup_t ls;
     real                   rlist_1x1;
     int                    n_nonlin_vsite;
@@ -1397,18 +1392,18 @@ static void set_verlet_buffer(const gmx_mtop_t *mtop,
     }
 
     printf("Calculated rlist for %dx%d atom pair-list as %.3f nm, buffer size %.3f nm\n",
-           1, 1, rlist_1x1, rlist_1x1-max(ir->rvdw, ir->rcoulomb));
+           1, 1, rlist_1x1, rlist_1x1-std::max(ir->rvdw, ir->rcoulomb));
 
     ir->rlistlong = ir->rlist;
     printf("Set rlist, assuming %dx%d atom pair-list, to %.3f nm, buffer size %.3f nm\n",
            ls.cluster_size_i, ls.cluster_size_j,
-           ir->rlist, ir->rlist-max(ir->rvdw, ir->rcoulomb));
+           ir->rlist, ir->rlist-std::max(ir->rvdw, ir->rcoulomb));
 
     printf("Note that mdrun will redetermine rlist based on the actual pair-list setup\n");
 
     if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
     {
-        gmx_fatal(FARGS, "The pair-list cut-off (%g nm) is longer than half the shortest box vector or longer than the smallest box diagonal element (%g nm). Increase the box size or decrease nstlist or increase verlet-buffer-tolerance.", ir->rlistlong, sqrt(max_cutoff2(ir->ePBC, box)));
+        gmx_fatal(FARGS, "The pair-list cut-off (%g nm) is longer than half the shortest box vector or longer than the smallest box diagonal element (%g nm). Increase the box size or decrease nstlist or increase verlet-buffer-tolerance.", ir->rlistlong, std::sqrt(max_cutoff2(ir->ePBC, box)));
     }
 }
 
@@ -1519,26 +1514,21 @@ int gmx_grompp(int argc, char *argv[])
     t_molinfo         *mi, *intermolecular_interactions;
     gpp_atomtype_t     atype;
     t_inputrec        *ir;
-    int                natoms, nvsite, comb, mt;
+    int                nvsite, comb, mt;
     t_params          *plist;
     t_state           *state;
     matrix             box;
-    real               max_spacing, fudgeQQ;
+    real               fudgeQQ;
     double             reppow;
     char               fn[STRLEN], fnB[STRLEN];
     const char        *mdparin;
     int                ntype;
     gmx_bool           bNeedVel, bGenVel;
     gmx_bool           have_atomnumber;
-    int                n12, n13, n14;
-    t_params          *gb_plist = NULL;
-    gmx_genborn_t     *born     = NULL;
     output_env_t       oenv;
     gmx_bool           bVerbose = FALSE;
     warninp_t          wi;
     char               warn_buf[STRLEN];
-    unsigned int       useed;
-    t_atoms            IMDatoms;   /* Atoms to be operated on interactively (IMD) */
 
     t_filenm           fnm[] = {
         { efMDP, NULL,  NULL,        ffREAD  },
@@ -1606,7 +1596,7 @@ int gmx_grompp(int argc, char *argv[])
     if (ir->ld_seed == -1)
     {
         ir->ld_seed = (gmx_int64_t)gmx_rng_make_seed();
-        fprintf(stderr, "Setting the LD random seed to %"GMX_PRId64 "\n", ir->ld_seed);
+        fprintf(stderr, "Setting the LD random seed to %" GMX_PRId64 "\n", ir->ld_seed);
     }
 
     if (ir->expandedvals->lmc_seed == -1)
@@ -1785,7 +1775,7 @@ int gmx_grompp(int argc, char *argv[])
     if (bRenum)
     {
         renum_atype(plist, sys, ir->wall_atomtype, atype, bVerbose);
-        ntype = get_atomtype_ntypes(atype);
+        get_atomtype_ntypes(atype);
     }
 
     if (ir->implicit_solvent != eisNO)
@@ -1865,7 +1855,6 @@ int gmx_grompp(int argc, char *argv[])
     }
     do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm),
              sys, bVerbose, ir,
-             bGenVel ? state->v : NULL,
              wi);
 
     if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0 &&
@@ -2012,8 +2001,8 @@ int gmx_grompp(int argc, char *argv[])
             set_warning_line(wi, mdparin, -1);
             warning_error(wi, "Some of the Fourier grid sizes are set, but all of them need to be set.");
         }
-        max_spacing = calc_grid(stdout, box, ir->fourier_spacing,
-                                &(ir->nkx), &(ir->nky), &(ir->nkz));
+        calc_grid(stdout, box, ir->fourier_spacing,
+                  &(ir->nkx), &(ir->nky), &(ir->nkz));
     }
 
     /* MRS: eventually figure out better logic for initializing the fep