Convert gmx_stochd_t to C++

Changed raw arrays to std::vector

Moved init functionality into constructor

Made owning gmx_update_t use a unique_ptr

Changed functions using struct to take constant references

Change-Id: I3664d9cc2b7e89eba57ad36a89e9c44b7895a7ef

diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp
index 71b7746cc2270658eda8755f8879c620f0fc5719..fbfd884ad7562e04cd5cbef8cbb5137ead505565 100644 (file)
--- a/src/gromacs/mdlib/coupling.cpp
+++ b/src/gromacs/mdlib/coupling.cpp
@@ -781,7 +781,8 @@ void berendsen_tcoupl(const t_inputrec *ir, const gmx_ekindata_t *ekind, real dt
 void andersen_tcoupl(const t_inputrec *ir, int64_t step,
                      const t_commrec *cr, const t_mdatoms *md,
                      gmx::ArrayRef<gmx::RVec> v,
-                     real rate, const gmx_bool *randomize, const real *boltzfac)
+                     real rate, const std::vector<bool> &randomize,
+                     gmx::ArrayRef<const real> boltzfac)
 {
     const int                                 *gatindex = (DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
     int                                        i;

diff --git a/src/gromacs/mdlib/update.cpp b/src/gromacs/mdlib/update.cpp
index c2a93a826ef305963e34a704049597dfef3a23a9..80aba49efa53b8d1f1d02cdd4309c62aba835495 100644 (file)
--- a/src/gromacs/mdlib/update.cpp
+++ b/src/gromacs/mdlib/update.cpp
@@ -43,6 +43,7 @@
 
 #include <algorithm>
 
+#include "gromacs/compat/make_unique.h"
 #include "gromacs/domdec/domdec_struct.h"
 #include "gromacs/fileio/confio.h"
 #include "gromacs/gmxlib/network.h"
@@ -92,22 +93,25 @@ typedef struct {
     real V;
 } gmx_sd_sigma_t;
 
-typedef struct {
+struct gmx_stochd_t
+{
     /* BD stuff */
-    real           *bd_rf;
+    std::vector<real>           bd_rf;
     /* SD stuff */
-    gmx_sd_const_t *sdc;
-    gmx_sd_sigma_t *sdsig;
+    std::vector<gmx_sd_const_t> sdc;
+    std::vector<gmx_sd_sigma_t> sdsig;
     /* andersen temperature control stuff */
-    gmx_bool       *randomize_group;
-    real           *boltzfac;
-} gmx_stochd_t;
+    std::vector<bool>           randomize_group;
+    std::vector<real>           boltzfac;
+
+    gmx_stochd_t(const t_inputrec *ir);
+};
 
 struct gmx_update_t
 {
-    gmx_stochd_t           *sd;
+    std::unique_ptr<gmx_stochd_t> sd;
     /* xprime for constraint algorithms */
-    PaddedVector<gmx::RVec> xp;
+    PaddedVector<gmx::RVec>       xp;
 
     /* Variables for the deform algorithm */
     int64_t           deformref_step;
@@ -763,25 +767,19 @@ static void do_update_vv_pos(int start, int nrend, real dt,
     }
 } /* do_update_vv_pos */
 
-static gmx_stochd_t *init_stochd(const t_inputrec *ir)
+gmx_stochd_t::gmx_stochd_t(const t_inputrec *ir)
 {
-    gmx_stochd_t   *sd;
-
-    snew(sd, 1);
-
     const t_grpopts *opts = &ir->opts;
-    int              ngtc = opts->ngtc;
+    const int        ngtc = opts->ngtc;
 
     if (ir->eI == eiBD)
     {
-        snew(sd->bd_rf, ngtc);
+        bd_rf.resize(ngtc);
     }
     else if (EI_SD(ir->eI))
     {
-        snew(sd->sdc, ngtc);
-        snew(sd->sdsig, ngtc);
-
-        gmx_sd_const_t *sdc = sd->sdc;
+        sdc.resize(ngtc);
+        sdsig.resize(ngtc);
 
         for (int gt = 0; gt < ngtc; gt++)
         {
@@ -798,8 +796,8 @@ static gmx_stochd_t *init_stochd(const t_inputrec *ir)
     }
     else if (ETC_ANDERSEN(ir->etc))
     {
-        snew(sd->randomize_group, ngtc);
-        snew(sd->boltzfac, ngtc);
+        randomize_group.resize(ngtc);
+        boltzfac.resize(ngtc);
 
         /* for now, assume that all groups, if randomized, are randomized at the same rate, i.e. tau_t is the same. */
         /* since constraint groups don't necessarily match up with temperature groups! This is checked in readir.c */
@@ -809,17 +807,15 @@ static gmx_stochd_t *init_stochd(const t_inputrec *ir)
             real reft = std::max<real>(0, opts->ref_t[gt]);
             if ((opts->tau_t[gt] > 0) && (reft > 0))  /* tau_t or ref_t = 0 means that no randomization is done */
             {
-                sd->randomize_group[gt] = TRUE;
-                sd->boltzfac[gt]        = BOLTZ*opts->ref_t[gt];
+                randomize_group[gt] = true;
+                boltzfac[gt]        = BOLTZ*opts->ref_t[gt];
             }
             else
             {
-                sd->randomize_group[gt] = FALSE;
+                randomize_group[gt] = false;
             }
         }
     }
-
-    return sd;
 }
 
 void update_temperature_constants(gmx_update_t *upd, const t_inputrec *ir)
@@ -857,10 +853,7 @@ gmx_update_t *init_update(const t_inputrec    *ir,
 {
     gmx_update_t *upd = new(gmx_update_t);
 
-    if (ir->eI == eiBD || EI_SD(ir->eI) || ir->etc == etcVRESCALE || ETC_ANDERSEN(ir->etc))
-    {
-        upd->sd    = init_stochd(ir);
-    }
+    upd->sd    = gmx::compat::make_unique<gmx_stochd_t>(ir);
 
     update_temperature_constants(upd, ir);
 
@@ -899,7 +892,7 @@ enum class SDUpdate : int
  * two with only one contribution, and one with both contributions. */
 template <SDUpdate updateType>
 static void
-doSDUpdateGeneral(gmx_stochd_t *sd,
+doSDUpdateGeneral(const gmx_stochd_t &sd,
                   int start, int nrend, real dt,
                   const rvec accel[], const ivec nFreeze[],
                   const real invmass[], const unsigned short ptype[],
@@ -930,9 +923,6 @@ doSDUpdateGeneral(gmx_stochd_t *sd,
     gmx::ThreeFry2x64<0> rng(seed, gmx::RandomDomain::UpdateCoordinates);
     gmx::TabulatedNormalDistribution<real, 14> dist;
 
-    gmx_sd_const_t *sdc = sd->sdc;
-    gmx_sd_sigma_t *sig = sd->sdsig;
-
     for (int n = start; n < nrend; n++)
     {
         int globalAtomIndex = gatindex ? gatindex[n] : n;
@@ -960,8 +950,8 @@ doSDUpdateGeneral(gmx_stochd_t *sd,
                 else if (updateType == SDUpdate::FrictionAndNoiseOnly)
                 {
                     real vn      = v[n][d];
-                    v[n][d]      = (vn*sdc[temperatureGroup].em +
-                                    invsqrtMass*sig[temperatureGroup].V*dist(rng));
+                    v[n][d]      = (vn*sd.sdc[temperatureGroup].em +
+                                    invsqrtMass*sd.sdsig[temperatureGroup].V*dist(rng));
                     // The previous phase already updated the
                     // positions with a full v*dt term that must
                     // now be half removed.
@@ -970,8 +960,8 @@ doSDUpdateGeneral(gmx_stochd_t *sd,
                 else
                 {
                     real vn      = v[n][d] + (inverseMass*f[n][d] + accel[accelerationGroup][d])*dt;
-                    v[n][d]      = (vn*sdc[temperatureGroup].em +
-                                    invsqrtMass*sig[temperatureGroup].V*dist(rng));
+                    v[n][d]      = (vn*sd.sdc[temperatureGroup].em +
+                                    invsqrtMass*sd.sdsig[temperatureGroup].V*dist(rng));
                     // Here we include half of the friction+noise
                     // update of v into the position update.
                     xprime[n][d] = x[n][d] + 0.5*(vn + v[n][d])*dt;
@@ -1568,7 +1558,7 @@ update_sd_second_half(int64_t                        step,
                 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
 
                 doSDUpdateGeneral<SDUpdate::FrictionAndNoiseOnly>
-                    (upd->sd,
+                    (*upd->sd,
                     start_th, end_th, dt,
                     inputrec->opts.acc, inputrec->opts.nFreeze,
                     md->invmass, md->ptype,
@@ -1846,7 +1836,7 @@ void update_coords(int64_t                             step,
                     {
                         // With constraints, the SD update is done in 2 parts
                         doSDUpdateGeneral<SDUpdate::ForcesOnly>
-                            (upd->sd,
+                            (*upd->sd,
                             start_th, end_th, dt,
                             inputrec->opts.acc, inputrec->opts.nFreeze,
                             md->invmass, md->ptype,
@@ -1857,7 +1847,7 @@ void update_coords(int64_t                             step,
                     else
                     {
                         doSDUpdateGeneral<SDUpdate::Combined>
-                            (upd->sd,
+                            (*upd->sd,
                             start_th, end_th, dt,
                             inputrec->opts.acc, inputrec->opts.nFreeze,
                             md->invmass, md->ptype,
@@ -1873,7 +1863,7 @@ void update_coords(int64_t                             step,
                                  md->cFREEZE, md->cTC,
                                  x_rvec, xp_rvec, v_rvec, f_rvec,
                                  inputrec->bd_fric,
-                                 upd->sd->bd_rf,
+                                 upd->sd->bd_rf.data(),
                                  step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
                     break;
                 case (eiVV):
@@ -1940,7 +1930,8 @@ extern gmx_bool update_randomize_velocities(const t_inputrec *ir, int64_t step,
     if ((ir->etc == etcANDERSEN) || do_per_step(step, roundToInt(1.0/rate)))
     {
         andersen_tcoupl(ir, step, cr, md, v, rate,
-                        upd->sd->randomize_group, upd->sd->boltzfac);
+                        upd->sd->randomize_group,
+                        upd->sd->boltzfac);
         return TRUE;
     }
     return FALSE;

diff --git a/src/gromacs/mdlib/update.h b/src/gromacs/mdlib/update.h
index be81265beac8f84f761c1718b1b4490afa16345d..da2a9b765b839c1e28d6f71e79984c7bc82dd16f 100644 (file)
--- a/src/gromacs/mdlib/update.h
+++ b/src/gromacs/mdlib/update.h
@@ -221,7 +221,8 @@ void berendsen_tcoupl(const t_inputrec *ir, const gmx_ekindata_t *ekind, real dt
 void andersen_tcoupl(const t_inputrec *ir, int64_t step,
                      const t_commrec *cr, const t_mdatoms *md,
                      gmx::ArrayRef<gmx::RVec> v,
-                     real rate, const gmx_bool *randomize, const real *boltzfac);
+                     real rate, const std::vector<bool> &randomize,
+                     gmx::ArrayRef<const real> boltzfac);
 
 void nosehoover_tcoupl(const t_grpopts *opts, const gmx_ekindata_t *ekind, real dt,
                        double xi[], double vxi[], const t_extmass *MassQ);