Moved another bunch of mdlib utils to C++
authorErik Lindahl <erik@kth.se>
Wed, 8 Jul 2015 14:32:06 +0000 (16:32 +0200)
committerDavid van der Spoel <spoel@xray.bmc.uu.se>
Wed, 29 Jul 2015 18:34:39 +0000 (20:34 +0200)
Changed calls to pow() where the second argument is 2 so we use
the sqr() function instead. Tables.c has been left as C since the
C++ version exposed a bug in icc-14 targeting MIC - this will be
addressed in a separate change.

Change-Id: Ie537519f4c7c7ab17a1aee575c00f83e3f71a068

src/gromacs/mdlib/adress.cpp [moved from src/gromacs/mdlib/adress.c with 94% similarity]
src/gromacs/mdlib/adress.h
src/gromacs/mdlib/expanded.cpp [moved from src/gromacs/mdlib/expanded.c with 88% similarity]
src/gromacs/mdlib/rf_util.cpp [moved from src/gromacs/mdlib/rf_util.c with 96% similarity]
src/gromacs/mdlib/sim_util.cpp
src/gromacs/mdlib/wall.cpp [moved from src/gromacs/mdlib/wall.c with 94% similarity]
src/gromacs/mdlib/wnblist.cpp [moved from src/gromacs/mdlib/wnblist.c with 89% similarity]

similarity index 94%
rename from src/gromacs/mdlib/adress.c
rename to src/gromacs/mdlib/adress.cpp
index 9e9db2908984f30cd6a842a4780d2979b194382e..4f5fb2a471dc8ee3d3af2e82491b337b40b63f4e 100644 (file)
@@ -2,7 +2,7 @@
  * This file is part of the GROMACS molecular simulation package.
  *
  * Copyright (c) 2009 Christoph Junghans, Brad Lambeth.
- * Copyright (c) 2011,2012,2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 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.
@@ -38,6 +38,8 @@
 
 #include "adress.h"
 
+#include <cmath>
+
 #include "gromacs/legacyheaders/typedefs.h"
 #include "gromacs/legacyheaders/types/simple.h"
 #include "gromacs/math/utilities.h"
@@ -95,7 +97,7 @@ adress_weight(rvec                 x,
             return 1;
     }
 
-    dl = sqrt(sqr_dl);
+    dl = std::sqrt(sqr_dl);
 
     /* molecule is coarse grained */
     if (dl > l2)
@@ -110,7 +112,7 @@ adress_weight(rvec                 x,
     /* hybrid region */
     else
     {
-        tmp = cos((dl-adressr)*M_PI/2/adressw);
+        tmp = std::cos((dl-adressr)*M_PI/2/adressw);
         return tmp*tmp;
     }
 }
@@ -261,14 +263,11 @@ void update_adress_weights_atom_per_atom(
         t_mdatoms *          mdatoms,
         t_pbc *              pbc)
 {
-    int            icg, k, k0, k1, d;
-    real           nrcg, inv_ncg, mtot, inv_mtot;
+    int            icg, k, k0, k1;
     atom_id *      cgindex;
-    rvec           ix;
     int            adresstype;
     real           adressr, adressw;
     rvec *         ref;
-    real *         massT;
     real *         wf;
 
 
@@ -281,7 +280,6 @@ void update_adress_weights_atom_per_atom(
     adresstype         = fr->adress_type;
     adressr            = fr->adress_ex_width;
     adressw            = fr->adress_hy_width;
-    massT              = mdatoms->massT;
     wf                 = mdatoms->wf;
     ref                = &(fr->adress_refs);
 
@@ -296,7 +294,6 @@ void update_adress_weights_atom_per_atom(
     {
         k0   = cgindex[icg];
         k1   = cgindex[icg + 1];
-        nrcg = k1 - k0;
 
         for (k = (k0); (k < k1); k++)
         {
@@ -326,7 +323,7 @@ update_adress_weights_cog(t_iparams            ip[],
                           t_mdatoms *          mdatoms,
                           t_pbc *              pbc)
 {
-    int            i, j, k, nr, nra, inc;
+    int            i, j, nr, nra, inc;
     int            ftype, adresstype;
     t_iatom        avsite, ai, aj, ak, al;
     t_iatom *      ia;
@@ -464,13 +461,11 @@ update_adress_weights_atom(int                  cg0,
     int            adresstype;
     real           adressr, adressw;
     rvec *         ref;
-    real *         massT;
     real *         wf;
 
     adresstype         = fr->adress_type;
     adressr            = fr->adress_ex_width;
     adressw            = fr->adress_hy_width;
-    massT              = mdatoms->massT;
     wf                 = mdatoms->wf;
     ref                = &(fr->adress_refs);
     cgindex            = cgs->index;
@@ -497,33 +492,25 @@ update_adress_weights_atom(int                  cg0,
 void
 adress_thermo_force(int                  start,
                     int                  homenr,
-                    t_block *            cgs,
                     rvec                 x[],
                     rvec                 f[],
                     t_forcerec *         fr,
                     t_mdatoms *          mdatoms,
                     t_pbc *              pbc)
 {
-    int              iatom, n0, nnn, nrcg, i;
+    int              iatom, n0, nnn, i;
     int              adresstype;
-    real             adressw, adressr;
-    atom_id *        cgindex;
     unsigned short * ptype;
     rvec *           ref;
-    real *           wf;
     real             tabscale;
     real *           ATFtab;
     rvec             dr;
-    real             w, wsq, wmin1, wmin1sq, wp, wt, rinv, sqr_dl, dl;
-    real             eps, eps2, F, Geps, Heps2, Fp, dmu_dwp, dwp_dr, fscal;
+    real             wt, rinv, sqr_dl, dl;
+    real             eps, eps2, F, Geps, Heps2, Fp, fscal;
 
     adresstype        = fr->adress_type;
-    adressw           = fr->adress_hy_width;
-    adressr           = fr->adress_ex_width;
-    cgindex           = cgs->index;
     ptype             = mdatoms->ptype;
     ref               = &(fr->adress_refs);
-    wf                = mdatoms->wf;
 
     for (iatom = start; (iatom < start+homenr); iatom++)
     {
@@ -531,7 +518,6 @@ adress_thermo_force(int                  start,
         {
             if (ptype[iatom] == eptVSite)
             {
-                w    = wf[iatom];
                 /* is it hybrid or apply the thermodynamics force everywhere?*/
                 if (mdatoms->tf_table_index[iatom] != NO_TF_TABLE)
                 {
@@ -548,7 +534,6 @@ adress_thermo_force(int                  start,
                         tabscale = fr->atf_tabs[DEFAULT_TF_TABLE].scale;
                     }
 
-                    fscal            = 0;
                     if (pbc)
                     {
                         pbc_dx(pbc, (*ref), x[iatom], dr);
@@ -583,10 +568,10 @@ adress_thermo_force(int                  start,
                             rinv = 0.0;
                     }
 
-                    dl = sqrt(sqr_dl);
+                    dl = std::sqrt(sqr_dl);
                     /* table origin is adress center */
                     wt               = dl*tabscale;
-                    n0               = wt;
+                    n0               = static_cast<int>(wt);
                     eps              = wt-n0;
                     eps2             = eps*eps;
                     nnn              = 4*n0;
index a25514d8b18e6c0b70129568f550e7203fa28a20..be88ab1d86a727fb8490dc41075a2f511e33a0b1 100644 (file)
@@ -2,7 +2,7 @@
  * This file is part of the GROMACS molecular simulation package.
  *
  * Copyright (c) 2009 Christoph Junghans, Brad Lambeth.
- * Copyright (c) 2011,2014, by the GROMACS development team, led by
+ * Copyright (c) 2011,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.
@@ -158,7 +158,6 @@ update_adress_weights_atom_per_atom(int                  cg0,
  *
  * \param[in] cg0 first charge group to update
  * \param[in] cg1 last+1 charge group to update
- * \param[in] cgs block containing the cg index
  * \param[in] x array with all the particle positions
  * \param[in,out] f the force array pointing at f_novirsum from sim_util.c
  * \param[in] fr the forcerec containing all the parameters
@@ -168,7 +167,6 @@ update_adress_weights_atom_per_atom(int                  cg0,
 void
 adress_thermo_force(int                  cg0,
                     int                  cg1,
-                    t_block *            cgs,
                     rvec                 x[],
                     rvec                 f[],
                     t_forcerec *         fr,
similarity index 88%
rename from src/gromacs/mdlib/expanded.c
rename to src/gromacs/mdlib/expanded.cpp
index 402662de4c31989a0eb6c13f3a657c3fe60cd7cd..6c53c76b1bc8c99340164f71edf6081d7368fb40 100644 (file)
  */
 #include "gmxpre.h"
 
-#include <math.h>
 #include <stdio.h>
 
+#include <cmath>
+
+#include <algorithm>
+
 #include "gromacs/domdec/domdec.h"
 #include "gromacs/fileio/confio.h"
 #include "gromacs/fileio/gmxfio.h"
@@ -104,12 +107,12 @@ static void GenerateGibbsProbabilities(real *ene, double *p_k, double *pks, int
     /* find the denominator */
     for (i = minfep; i <= maxfep; i++)
     {
-        *pks += exp(ene[i]-maxene);
+        *pks += std::exp(ene[i]-maxene);
     }
     /*numerators*/
     for (i = minfep; i <= maxfep; i++)
     {
-        p_k[i] = exp(ene[i]-maxene) / *pks;
+        p_k[i] = std::exp(ene[i]-maxene) / *pks;
     }
 }
 
@@ -128,11 +131,11 @@ static void GenerateWeightedGibbsProbabilities(real *ene, double *p_k, double *p
         {
             /* add the delta, since we need to make sure it's greater than zero, and
                we need a non-arbitrary number? */
-            nene[i] = ene[i] + log(nvals[i]+delta);
+            nene[i] = ene[i] + std::log(nvals[i]+delta);
         }
         else
         {
-            nene[i] = ene[i] + log(nvals[i]);
+            nene[i] = ene[i] + std::log(nvals[i]);
         }
     }
 
@@ -155,13 +158,13 @@ static void GenerateWeightedGibbsProbabilities(real *ene, double *p_k, double *p
     /* find the denominator */
     for (i = 0; i < nlim; i++)
     {
-        *pks += exp(nene[i]);
+        *pks += std::exp(nene[i]);
     }
 
     /*numerators*/
     for (i = 0; i < nlim; i++)
     {
-        p_k[i] = exp(nene[i]) / *pks;
+        p_k[i] = std::exp(nene[i]) / *pks;
     }
     sfree(nene);
 }
@@ -180,18 +183,18 @@ real do_logsum(int N, real *a_n)
     maxarg = a_n[0];
     for (i = 1; i < N; i++)
     {
-        maxarg = max(maxarg, a_n[i]);
+        maxarg = std::max(maxarg, a_n[i]);
     }
 
     /* compute sum of exp(a_n - maxarg) */
     sum = 0.0;
     for (i = 0; i < N; i++)
     {
-        sum = sum + exp(a_n[i] - maxarg);
+        sum = sum + std::exp(a_n[i] - maxarg);
     }
 
     /*     compute log sum */
-    logsum = log(sum) + maxarg;
+    logsum = std::log(sum) + maxarg;
     return logsum;
 }
 
@@ -256,95 +259,102 @@ static gmx_bool CheckIfDoneEquilibrating(int nlim, t_expanded *expand, df_histor
     gmx_bool bDoneEquilibrating = TRUE;
     gmx_bool bIfFlat;
 
-    /* assume we have equilibrated the weights, then check to see if any of the conditions are not met */
-
-    /* calculate the total number of samples */
-    switch (expand->elmceq)
+    /* If we are doing slow growth to get initial values, we haven't finished equilibrating */
+    if (expand->lmc_forced_nstart > 0)
     {
-        case elmceqNO:
-            /* We have not equilibrated, and won't, ever. */
-            return FALSE;
-        case elmceqYES:
-            /* we have equilibrated -- we're done */
-            return TRUE;
-        case elmceqSTEPS:
-            /* first, check if we are equilibrating by steps, if we're still under */
-            if (step < expand->equil_steps)
+        for (i = 0; i < nlim; i++)
+        {
+            if (dfhist->n_at_lam[i] < expand->lmc_forced_nstart) /* we are still doing the initial sweep, so we're definitely not
+                                                                    done equilibrating*/
             {
                 bDoneEquilibrating = FALSE;
+                break;
             }
-            break;
-        case elmceqSAMPLES:
-            totalsamples = 0;
-            for (i = 0; i < nlim; i++)
-            {
-                totalsamples += dfhist->n_at_lam[i];
-            }
-            if (totalsamples < expand->equil_samples)
-            {
+        }
+    }
+    else
+    {
+        /* assume we have equilibrated the weights, then check to see if any of the conditions are not met */
+        bDoneEquilibrating = TRUE;
+
+        /* calculate the total number of samples */
+        switch (expand->elmceq)
+        {
+            case elmceqNO:
+                /* We have not equilibrated, and won't, ever. */
                 bDoneEquilibrating = FALSE;
-            }
-            break;
-        case elmceqNUMATLAM:
-            for (i = 0; i < nlim; i++)
-            {
-                if (dfhist->n_at_lam[i] < expand->equil_n_at_lam) /* we are still doing the initial sweep, so we're definitely not
-                                                                     done equilibrating*/
+                break;
+            case elmceqYES:
+                /* we have equilibrated -- we're done */
+                bDoneEquilibrating = TRUE;
+                break;
+            case elmceqSTEPS:
+                /* first, check if we are equilibrating by steps, if we're still under */
+                if (step < expand->equil_steps)
                 {
-                    bDoneEquilibrating  = FALSE;
-                    break;
+                    bDoneEquilibrating = FALSE;
                 }
-            }
-            break;
-        case elmceqWLDELTA:
-            if (EWL(expand->elamstats)) /* This check is in readir as well, but
-                                           just to be sure */
-            {
-                if (dfhist->wl_delta > expand->equil_wl_delta)
+                break;
+            case elmceqSAMPLES:
+                totalsamples = 0;
+                for (i = 0; i < nlim; i++)
+                {
+                    totalsamples += dfhist->n_at_lam[i];
+                }
+                if (totalsamples < expand->equil_samples)
                 {
                     bDoneEquilibrating = FALSE;
                 }
-            }
-            break;
-        case elmceqRATIO:
-            /* we can use the flatness as a judge of good weights, as long as
-               we're not doing minvar, or Wang-Landau.
-               But turn off for now until we figure out exactly how we do this.
-             */
-
-            if (!(EWL(expand->elamstats) || expand->elamstats == elamstatsMINVAR))
-            {
-                /* we want to use flatness -avoiding- the forced-through samples.  Plus, we need to convert to
-                   floats for this histogram function. */
-
-                real *modhisto;
-                snew(modhisto, nlim);
+                break;
+            case elmceqNUMATLAM:
                 for (i = 0; i < nlim; i++)
                 {
-                    modhisto[i] = 1.0*(dfhist->n_at_lam[i]-expand->lmc_forced_nstart);
+                    if (dfhist->n_at_lam[i] < expand->equil_n_at_lam) /* we are still doing the initial sweep, so we're definitely not
+                                                                         done equilibrating*/
+                    {
+                        bDoneEquilibrating  = FALSE;
+                        break;
+                    }
                 }
-                bIfFlat = CheckHistogramRatios(nlim, modhisto, expand->equil_ratio);
-                sfree(modhisto);
-                if (!bIfFlat)
+                break;
+            case elmceqWLDELTA:
+                if (EWL(expand->elamstats)) /* This check is in readir as well, but
+                                               just to be sure */
                 {
-                    bDoneEquilibrating = FALSE;
+                    if (dfhist->wl_delta > expand->equil_wl_delta)
+                    {
+                        bDoneEquilibrating = FALSE;
+                    }
                 }
-            }
-        default:
-            bDoneEquilibrating = TRUE;
-    }
-    /* one last case to go though, if we are doing slow growth to get initial values, we haven't finished equilibrating */
+                break;
+            case elmceqRATIO:
+                /* we can use the flatness as a judge of good weights, as long as
+                   we're not doing minvar, or Wang-Landau.
+                   But turn off for now until we figure out exactly how we do this.
+                 */
 
-    if (expand->lmc_forced_nstart > 0)
-    {
-        for (i = 0; i < nlim; i++)
-        {
-            if (dfhist->n_at_lam[i] < expand->lmc_forced_nstart) /* we are still doing the initial sweep, so we're definitely not
-                                                                    done equilibrating*/
-            {
-                bDoneEquilibrating = FALSE;
+                if (!(EWL(expand->elamstats) || expand->elamstats == elamstatsMINVAR))
+                {
+                    /* we want to use flatness -avoiding- the forced-through samples.  Plus, we need to convert to
+                       floats for this histogram function. */
+
+                    real *modhisto;
+                    snew(modhisto, nlim);
+                    for (i = 0; i < nlim; i++)
+                    {
+                        modhisto[i] = 1.0*(dfhist->n_at_lam[i]-expand->lmc_forced_nstart);
+                    }
+                    bIfFlat = CheckHistogramRatios(nlim, modhisto, expand->equil_ratio);
+                    sfree(modhisto);
+                    if (!bIfFlat)
+                    {
+                        bDoneEquilibrating = FALSE;
+                    }
+                }
+                break;
+            default:
+                bDoneEquilibrating = TRUE;
                 break;
-            }
         }
     }
     return bDoneEquilibrating;
@@ -353,21 +363,18 @@ static gmx_bool CheckIfDoneEquilibrating(int nlim, t_expanded *expand, df_histor
 static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist,
                               int fep_state, real *scaled_lamee, real *weighted_lamee, gmx_int64_t step)
 {
-    real     maxdiff = 0.000000001;
-    gmx_bool bSufficientSamples;
-    int      i, k, n, nz, indexi, indexk, min_n, max_n, totali;
-    int      n0, np1, nm1, nval, min_nvalm, min_nvalp, maxc;
-    real     omega_m1_0, omega_p1_m1, omega_m1_p1, omega_p1_0, clam_osum;
-    real     de, de_function, dr, denom, maxdr;
-    real     min_val, cnval, zero_sum_weights;
-    real    *omegam_array, *weightsm_array, *omegap_array, *weightsp_array, *varm_array, *varp_array, *dwp_array, *dwm_array;
-    real     clam_varm, clam_varp, clam_weightsm, clam_weightsp, clam_minvar;
-    real    *lam_weights, *lam_minvar_corr, *lam_variance, *lam_dg;
-    double  *p_k;
-    double   pks = 0;
-    real    *numweighted_lamee, *logfrac;
-    int     *nonzero;
-    real     chi_m1_0, chi_p1_0, chi_m2_0, chi_p2_0, chi_p1_m1, chi_p2_m1, chi_m1_p1, chi_m2_p1;
+    gmx_bool  bSufficientSamples;
+    int       i;
+    int       n0, np1, nm1, nval, min_nvalm, min_nvalp, maxc;
+    real      omega_m1_0, omega_p1_m1, omega_m1_p1, omega_p1_0, clam_osum;
+    real      de, de_function;
+    real      cnval, zero_sum_weights;
+    real     *omegam_array, *weightsm_array, *omegap_array, *weightsp_array, *varm_array, *varp_array, *dwp_array, *dwm_array;
+    real      clam_varm, clam_varp, clam_weightsm, clam_weightsp, clam_minvar;
+    real     *lam_variance, *lam_dg;
+    double   *p_k;
+    double    pks = 0;
+    real      chi_m1_0, chi_p1_0, chi_m2_0, chi_p2_0, chi_p1_m1, chi_p2_m1, chi_m1_p1, chi_m2_p1;
 
     /* if we have equilibrated the weights, exit now */
     if (dfhist->bEquil)
@@ -459,7 +466,7 @@ static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist
         for (i = 0; i < nlim-1; i++)
         {   /* only through the second to last */
             lam_dg[i]       = dfhist->sum_dg[i+1] - dfhist->sum_dg[i];
-            lam_variance[i] = pow(dfhist->sum_variance[i+1], 2) - pow(dfhist->sum_variance[i], 2);
+            lam_variance[i] = sqr(dfhist->sum_variance[i+1]) - sqr(dfhist->sum_variance[i]);
         }
 
         /* accumulate running averages */
@@ -470,7 +477,7 @@ static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist
             /* actually, should be able to rewrite it w/o exponential, for better numerical stability */
             if (fep_state > 0)
             {
-                de = exp(cnval - (scaled_lamee[fep_state]-scaled_lamee[fep_state-1]));
+                de = std::exp(cnval - (scaled_lamee[fep_state]-scaled_lamee[fep_state-1]));
                 if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMINVAR)
                 {
                     de_function = 1.0/(1.0+de);
@@ -492,7 +499,7 @@ static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist
 
             if (fep_state < nlim-1)
             {
-                de = exp(-cnval + (scaled_lamee[fep_state+1]-scaled_lamee[fep_state]));
+                de = std::exp(-cnval + (scaled_lamee[fep_state+1]-scaled_lamee[fep_state]));
                 if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMINVAR)
                 {
                     de_function = 1.0/(1.0+de);
@@ -574,7 +581,7 @@ static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist
                 }
                 if ((n0 > 0) && (nm1 > 0))
                 {
-                    clam_weightsm = (log(chi_m1_0) - log(chi_p1_m1)) + cnval;
+                    clam_weightsm = (std::log(chi_m1_0) - std::log(chi_p1_m1)) + cnval;
                     clam_varm     = (1.0/n0)*(omega_m1_0) + (1.0/nm1)*(omega_p1_m1);
                 }
             }
@@ -591,7 +598,7 @@ static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist
                 }
                 if ((n0 > 0) && (np1 > 0))
                 {
-                    clam_weightsp = (log(chi_m1_p1) - log(chi_p1_0)) + cnval;
+                    clam_weightsp = (std::log(chi_m1_p1) - std::log(chi_p1_0)) + cnval;
                     clam_varp     = (1.0/np1)*(omega_m1_p1) + (1.0/n0)*(omega_p1_0);
                 }
             }
@@ -608,7 +615,7 @@ static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist
             varm_array[nval]               = clam_varm;
             if (nm1 > 0)
             {
-                dwm_array[nval]  = fabs( (cnval + log((1.0*n0)/nm1)) - lam_dg[fep_state-1] );
+                dwm_array[nval]  = fabs( (cnval + std::log((1.0*n0)/nm1)) - lam_dg[fep_state-1] );
             }
             else
             {
@@ -627,7 +634,7 @@ static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist
             varp_array[nval]               = clam_varp;
             if ((np1 > 0) && (n0 > 0))
             {
-                dwp_array[nval]  = fabs( (cnval + log((1.0*np1)/n0)) - lam_dg[fep_state] );
+                dwp_array[nval]  = fabs( (cnval + std::log((1.0*np1)/n0)) - lam_dg[fep_state] );
             }
             else
             {
@@ -652,7 +659,7 @@ static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist
         clam_minvar = 0;
         if (clam_osum > 0)
         {
-            clam_minvar = 0.5*log(clam_osum);
+            clam_minvar = 0.5*std::log(clam_osum);
         }
 
         if (fep_state > 0)
@@ -705,7 +712,7 @@ static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist
         for (i = 1; i < nlim; i++)
         {
             dfhist->sum_dg[i]       = lam_dg[i-1] + dfhist->sum_dg[i-1];
-            dfhist->sum_variance[i] = sqrt(lam_variance[i-1] + pow(dfhist->sum_variance[i-1], 2));
+            dfhist->sum_variance[i] = std::sqrt(lam_variance[i-1] + sqr(dfhist->sum_variance[i-1]));
             dfhist->sum_weights[i]  = dfhist->sum_dg[i] + dfhist->sum_minvar[i];
         }
 
@@ -730,13 +737,11 @@ static int ChooseNewLambda(int nlim, t_expanded *expand, df_history_t *dfhist, i
 {
     /* Choose new lambda value, and update transition matrix */
 
-    int      i, ifep, jfep, minfep, maxfep, lamnew, lamtrial, starting_fep_state;
-    real     r1, r2, de_old, de_new, de, trialprob, tprob = 0;
-    real   **Tij;
+    int      i, ifep, minfep, maxfep, lamnew, lamtrial, starting_fep_state;
+    real     r1, r2, de, trialprob, tprob = 0;
     double  *propose, *accept, *remainder;
     double   pks;
-    real     sum, pnorm;
-    gmx_bool bRestricted;
+    real     pnorm;
 
     starting_fep_state = fep_state;
     lamnew             = fep_state; /* so that there is a default setting -- stays the same */
@@ -784,13 +789,11 @@ static int ChooseNewLambda(int nlim, t_expanded *expand, df_history_t *dfhist, i
 
         if ((expand->elmcmove == elmcmoveGIBBS) || (expand->elmcmove == elmcmoveMETGIBBS))
         {
-            bRestricted = TRUE;
             /* use the Gibbs sampler, with restricted range */
             if (expand->gibbsdeltalam < 0)
             {
                 minfep      = 0;
                 maxfep      = nlim-1;
-                bRestricted = FALSE;
             }
             else
             {
@@ -968,7 +971,7 @@ static int ChooseNewLambda(int nlim, t_expanded *expand, df_history_t *dfhist, i
             if (expand->elmcmove == elmcmoveMETROPOLIS)
             {
                 tprob     = 1.0;
-                trialprob = exp(de);
+                trialprob = std::exp(de);
                 if (trialprob < tprob)
                 {
                     tprob = trialprob;
@@ -981,7 +984,7 @@ static int ChooseNewLambda(int nlim, t_expanded *expand, df_history_t *dfhist, i
             }
             else if (expand->elmcmove == elmcmoveBARKER)
             {
-                tprob = 1.0/(1.0+exp(-de));
+                tprob = 1.0/(1.0+std::exp(-de));
 
                 propose[fep_state] = (1-tprob);
                 propose[lamtrial] += tprob; /* we add, to account for the fact that at the end, they might be the same point */
@@ -1022,8 +1025,7 @@ extern void PrintFreeEnergyInfoToFile(FILE *outfile, t_lambda *fep, t_expanded *
                                       int fep_state, int frequency, gmx_int64_t step)
 {
     int         nlim, i, ifep, jfep;
-    real        dw, dg, dv, dm, Tprint;
-    real       *temps;
+    real        dw, dg, dv, Tprint;
     const char *print_names[efptNR] = {" FEPL", "MassL", "CoulL", " VdwL", "BondL", "RestT", "Temp.(K)"};
     gmx_bool    bSimTemp            = FALSE;
 
@@ -1033,7 +1035,7 @@ extern void PrintFreeEnergyInfoToFile(FILE *outfile, t_lambda *fep, t_expanded *
         bSimTemp = TRUE;
     }
 
-    if (mod(step, frequency) == 0)
+    if (step % frequency == 0)
     {
         fprintf(outfile, "             MC-lambda information\n");
         if (EWL(expand->elamstats) && (!(dfhist->bEquil)))
@@ -1068,15 +1070,12 @@ extern void PrintFreeEnergyInfoToFile(FILE *outfile, t_lambda *fep, t_expanded *
                 dw = 0.0;
                 dg = 0.0;
                 dv = 0.0;
-                dm = 0.0;
             }
             else
             {
                 dw = dfhist->sum_weights[ifep+1] - dfhist->sum_weights[ifep];
                 dg = dfhist->sum_dg[ifep+1] - dfhist->sum_dg[ifep];
-                dv = sqrt(pow(dfhist->sum_variance[ifep+1], 2) - pow(dfhist->sum_variance[ifep], 2));
-                dm = dfhist->sum_minvar[ifep+1] - dfhist->sum_minvar[ifep];
-
+                dv = std::sqrt(sqr(dfhist->sum_variance[ifep+1]) - sqr(dfhist->sum_variance[ifep]));
             }
             fprintf(outfile, "%3d", (ifep+1));
             for (i = 0; i < efptNR; i++)
@@ -1124,7 +1123,7 @@ extern void PrintFreeEnergyInfoToFile(FILE *outfile, t_lambda *fep, t_expanded *
         }
         fprintf(outfile, "\n");
 
-        if ((mod(step, expand->nstTij) == 0) && (expand->nstTij > 0) && (step > 0))
+        if ((step % expand->nstTij == 0) && (expand->nstTij > 0) && (step > 0))
         {
             fprintf(outfile, "                     Transition Matrix\n");
             for (ifep = 0; ifep < nlim; ifep++)
@@ -1202,7 +1201,6 @@ extern int ExpandedEnsembleDynamics(FILE *log, t_inputrec *ir, gmx_enerdata_t *e
     real        oneovert, maxscaled = 0, maxweighted = 0;
     t_expanded *expand;
     t_simtemp  *simtemp;
-    double     *temperature_lambdas;
     gmx_bool    bIfReset, bSwitchtoOneOverT, bDoneEquilibrating = FALSE;
 
     expand  = ir->expandedvals;
@@ -1322,7 +1320,7 @@ extern int ExpandedEnsembleDynamics(FILE *log, t_inputrec *ir, gmx_enerdata_t *e
             {
                 told              = ir->opts.ref_t[i];
                 ir->opts.ref_t[i] =  simtemp->temperatures[lamnew];
-                buf_ngtc[i]       = sqrt(ir->opts.ref_t[i]/told); /* using the buffer as temperature scaling */
+                buf_ngtc[i]       = std::sqrt(ir->opts.ref_t[i]/told); /* using the buffer as temperature scaling */
             }
         }
 
similarity index 96%
rename from src/gromacs/mdlib/rf_util.c
rename to src/gromacs/mdlib/rf_util.cpp
index a4703cca24ebd8af1d5be38ae523119e038dd3f0..0058ed7f6c0989e4d845d3285ef03553958dc53e 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 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.
@@ -36,6 +36,8 @@
  */
 #include "gmxpre.h"
 
+#include <cmath>
+
 #include "gromacs/legacyheaders/copyrite.h"
 #include "gromacs/legacyheaders/force.h"
 #include "gromacs/legacyheaders/names.h"
@@ -55,7 +57,7 @@ real RF_excl_correction(const t_forcerec *fr, t_graph *g,
      * epsfac q_i q_j (k_rf r_ij^2 - c_rf)
      * and force correction for all excluded pairs, including self pairs.
      */
-    int         top, i, j, j1, j2, k, ki;
+    int         i, j, j1, j2, k, ki;
     double      q2sumA, q2sumB, ener;
     const real *chargeA, *chargeB;
     real        ek, ec, L1, qiA, qiB, qqA, qqB, qqL, v;
@@ -221,7 +223,7 @@ void calc_rffac(FILE *fplog, int eel, real eps_r, real eps_rf, real Rc, real Tem
             }
             /* Ionic strength (only needed for eelGRF */
             I       = 0.5*zsq/vol;
-            *kappa  = sqrt(2*I/(EPSILON0*eps_rf*BOLTZ*Temp));
+            *kappa  = std::sqrt(2*I/(EPSILON0*eps_rf*BOLTZ*Temp));
         }
         else
         {
@@ -242,7 +244,8 @@ void calc_rffac(FILE *fplog, int eel, real eps_r, real eps_rf, real Rc, real Tem
             *krf = ((eps_rf - eps_r)*k1 + 0.5*k2)/((2*eps_rf + eps_r)*k1 + k2)/(Rc*Rc*Rc);
         }
         *crf   = 1/Rc + *krf*Rc*Rc;
-        rmin   = pow(*krf*2.0, -1.0/3.0);
+        // Make sure we don't lose resolution in pow() by casting real arg to double
+        rmin   = std::pow(static_cast<double>(*krf*2.0), -1.0/3.0);
 
         if (fplog)
         {
index b2d98a11a017e75e3b4021aab4b9b7b85371b075..357819c491f2cb6c5c6aee427808dea04782c17f 100644 (file)
@@ -1871,7 +1871,7 @@ void do_force_cutsGROUP(FILE *fplog, t_commrec *cr,
         if (bDoAdressWF && fr->adress_icor == eAdressICThermoForce)
         {
             /* Compute thermodynamic force in hybrid AdResS region */
-            adress_thermo_force(start, homenr, &(top->cgs), x, fr->f_novirsum, fr, mdatoms,
+            adress_thermo_force(start, homenr, x, fr->f_novirsum, fr, mdatoms,
                                 inputrec->ePBC == epbcNONE ? NULL : &pbc);
         }
 
similarity index 94%
rename from src/gromacs/mdlib/wall.c
rename to src/gromacs/mdlib/wall.cpp
index bdc5c90e404b06e96cd38a06e86a28b205c7f669..c8aee134c6338fe60cef51c52a5f243e61450f31 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2008, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 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.
 
 #include "gmxpre.h"
 
-#include <math.h>
 #include <string.h>
 
+#include <algorithm>
+
 #include "gromacs/fileio/filenm.h"
 #include "gromacs/legacyheaders/force.h"
 #include "gromacs/legacyheaders/macros.h"
@@ -55,7 +56,7 @@ void make_wall_tables(FILE *fplog, const output_env_t oenv,
                       const gmx_groups_t *groups,
                       t_forcerec *fr)
 {
-    int           w, negp_pp, egp, i, j;
+    int           negp_pp;
     int          *nm_ind;
     char          buf[STRLEN];
     t_forcetable *tab;
@@ -70,10 +71,10 @@ void make_wall_tables(FILE *fplog, const output_env_t oenv,
     }
 
     snew(fr->wall_tab, ir->nwall);
-    for (w = 0; w < ir->nwall; w++)
+    for (int w = 0; w < ir->nwall; w++)
     {
         snew(fr->wall_tab[w], negp_pp);
-        for (egp = 0; egp < negp_pp; egp++)
+        for (int egp = 0; egp < negp_pp; egp++)
         {
             /* If the energy group pair is excluded, we don't need a table */
             if (!(fr->egp_flags[egp*ir->opts.ngener+negp_pp+w] & EGP_EXCL))
@@ -86,9 +87,9 @@ void make_wall_tables(FILE *fplog, const output_env_t oenv,
                         ftp2ext(efXVG));
                 *tab = make_tables(fplog, oenv, fr, FALSE, buf, 0, GMX_MAKETABLES_FORCEUSER);
                 /* Since wall have no charge, we can compress the table */
-                for (i = 0; i <= tab->n; i++)
+                for (int i = 0; i <= tab->n; i++)
                 {
-                    for (j = 0; j < 8; j++)
+                    for (int j = 0; j < 8; j++)
                     {
                         tab->data[8*i+j] = tab->data[12*i+4+j];
                     }
@@ -109,13 +110,13 @@ static void wall_error(int a, rvec *x, real r)
 real do_walls(t_inputrec *ir, t_forcerec *fr, matrix box, t_mdatoms *md,
               rvec x[], rvec f[], real lambda, real Vlj[], t_nrnb *nrnb)
 {
-    int             nwall, w, lam, i;
+    int             nwall;
     int             ntw[2], at, ntype, ngid, ggid, *egp_flags, *type;
-    real           *nbfp, lamfac, fac_d[2], fac_r[2], Cd, Cr, Vtot, Fwall[2];
+    real           *nbfp, lamfac, fac_d[2], fac_r[2], Cd, Cr, Vtot;
     real            wall_z[2], r, mr, r1, r2, r4, Vd, Vr, V = 0, Fd, Fr, F = 0, dvdlambda;
     dvec            xf_z;
     int             n0, nnn;
-    real            tabscale, *VFtab, rt, eps, eps2, Yt, Ft, Geps, Heps, Heps2, Fp, VV, FF;
+    real            tabscale, *VFtab, rt, eps, eps2, Yt, Ft, Geps, Heps2, Fp, VV, FF;
     unsigned short *gid = md->cENER;
     t_forcetable   *tab;
 
@@ -125,7 +126,7 @@ real do_walls(t_inputrec *ir, t_forcerec *fr, matrix box, t_mdatoms *md,
     nbfp      = fr->nbfp;
     egp_flags = fr->egp_flags;
 
-    for (w = 0; w < nwall; w++)
+    for (int w = 0; w < nwall; w++)
     {
         ntw[w] = 2*ntype*ir->wall_atomtype[w];
         switch (ir->wall_type)
@@ -141,7 +142,6 @@ real do_walls(t_inputrec *ir, t_forcerec *fr, matrix box, t_mdatoms *md,
             default:
                 break;
         }
-        Fwall[w] = 0;
     }
     wall_z[0] = 0;
     wall_z[1] = box[ZZ][ZZ];
@@ -149,7 +149,7 @@ real do_walls(t_inputrec *ir, t_forcerec *fr, matrix box, t_mdatoms *md,
     Vtot      = 0;
     dvdlambda = 0;
     clear_dvec(xf_z);
-    for (lam = 0; lam < (md->nPerturbed ? 2 : 1); lam++)
+    for (int lam = 0; lam < (md->nPerturbed ? 2 : 1); lam++)
     {
         if (md->nPerturbed)
         {
@@ -169,9 +169,9 @@ real do_walls(t_inputrec *ir, t_forcerec *fr, matrix box, t_mdatoms *md,
             lamfac = 1;
             type   = md->typeA;
         }
-        for (i = 0; i < md->homenr; i++)
+        for (int i = 0; i < md->homenr; i++)
         {
-            for (w = 0; w < nwall; w++)
+            for (int w = 0; w < std::min(nwall, 2); w++)
             {
                 /* The wall energy groups are always at the end of the list */
                 ggid = gid[i]*ngid + ngid - nwall + w;
@@ -210,7 +210,7 @@ real do_walls(t_inputrec *ir, t_forcerec *fr, matrix box, t_mdatoms *md,
                             VFtab    = tab->data;
 
                             rt    = r*tabscale;
-                            n0    = rt;
+                            n0    = static_cast<int>(rt);
                             if (n0 >= tab->n)
                             {
                                 /* Beyond the table range, set V and F to zero */
@@ -323,7 +323,7 @@ real do_walls(t_inputrec *ir, t_forcerec *fr, matrix box, t_mdatoms *md,
         inc_nrnb(nrnb, eNR_WALLS, md->homenr);
     }
 
-    for (i = 0; i < DIM; i++)
+    for (int i = 0; i < DIM; i++)
     {
         fr->vir_wall_z[i] = -0.5*xf_z[i];
     }
similarity index 89%
rename from src/gromacs/mdlib/wnblist.c
rename to src/gromacs/mdlib/wnblist.cpp
index 0e72450f564f0229d872adcc4a345210ddc6d48a..89c030c71736aa0c9d95130ee77bfc09fb962898 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 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.
@@ -40,6 +40,8 @@
 #include <stdio.h>
 #include <string.h>
 
+#include <algorithm>
+
 #include "gromacs/domdec/domdec.h"
 #include "gromacs/fileio/gmxfio.h"
 #include "gromacs/legacyheaders/force.h"
@@ -75,7 +77,7 @@ static void write_nblist(FILE *out, gmx_domdec_t *dd, t_nblist *nblist, int nDNL
                 ca1[zi] = dd->cgindex[dd_zones->cg_range[zi+1]];
             }
             i = 0;
-            for (zi = 0; zi < dd_zones->nizone; zi++)
+            for (zi = 0; zi < dd_zones->nizone && zi < dd_zones->n; zi++)
             {
                 zj0 = dd_zones->izone[zi].j0;
                 zj1 = dd_zones->izone[zi].j1;
@@ -160,34 +162,6 @@ static void set_mat(FILE *fp, int **mat, int i0, int ni, int j0, int nj,
 
 void dump_nblist(FILE *out, t_commrec *cr, t_forcerec *fr, int nDNL)
 {
-#if 0
-    static FILE *fp = NULL;
-    char         buf[STRLEN];
-    int          n, i;
-
-    if (fp == NULL)
-    {
-        if (PAR(cr))
-        {
-            sprintf(buf, "nlist_n%d.txt", cr->nodeid);
-        }
-        else
-        {
-            sprintf(buf, "nlist.txt");
-        }
-        fp = gmx_fio_fopen(buf, "w");
-    }
-    fprintf(fp, "%s\n", header);
-
-    for (n = 0; (n < fr->nnblists); n++)
-    {
-        for (i = 0; (i < eNL_NR); i++)
-        {
-            write_nblist(fp, cr->dd, &fr->nblists[n].nlist_sr[i], nDNL);
-        }
-    }
-#endif
-    char buf[STRLEN];
     int  n, i;
 
     fprintf(out, "%s\n", header);