From 62b11995e20adc5688224bd144e89d4ef30ebd3a Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Thu, 20 Dec 2018 09:43:30 -0500 Subject: [PATCH] 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 --- src/gromacs/mdlib/coupling.cpp | 3 +- src/gromacs/mdlib/update.cpp | 79 +++++++++++++++++++----------------------- src/gromacs/mdlib/update.h | 3 +- 3 files changed, 39 insertions(+), 46 deletions(-) diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp index 71b7746cc2..fbfd884ad7 100644 --- 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 v, - real rate, const gmx_bool *randomize, const real *boltzfac) + real rate, const std::vector &randomize, + gmx::ArrayRef 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 c2a93a826e..80aba49efa 100644 --- a/src/gromacs/mdlib/update.cpp +++ b/src/gromacs/mdlib/update.cpp @@ -43,6 +43,7 @@ #include +#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 bd_rf; /* SD stuff */ - gmx_sd_const_t *sdc; - gmx_sd_sigma_t *sdsig; + std::vector sdc; + std::vector sdsig; /* andersen temperature control stuff */ - gmx_bool *randomize_group; - real *boltzfac; -} gmx_stochd_t; + std::vector randomize_group; + std::vector boltzfac; + + gmx_stochd_t(const t_inputrec *ir); +}; struct gmx_update_t { - gmx_stochd_t *sd; + std::unique_ptr sd; /* xprime for constraint algorithms */ - PaddedVector xp; + PaddedVector 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(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(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 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 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 - (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 - (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 - (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 be81265bea..da2a9b765b 100644 --- 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 v, - real rate, const gmx_bool *randomize, const real *boltzfac); + real rate, const std::vector &randomize, + gmx::ArrayRef boltzfac); void nosehoover_tcoupl(const t_grpopts *opts, const gmx_ekindata_t *ekind, real dt, double xi[], double vxi[], const t_extmass *MassQ); -- 2.11.4.GIT