*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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 <algorithm>
+#include "gromacs/compat/make_unique.h"
#include "gromacs/domdec/domdec_struct.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/gmxlib/network.h"
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;
}
} /* 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++)
{
}
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 */
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)
{
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);
* 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[],
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;
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.
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;
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,
{
// 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,
else
{
doSDUpdateGeneral<SDUpdate::Combined>
- (upd->sd,
+ (*upd->sd,
start_th, end_th, dt,
inputrec->opts.acc, inputrec->opts.nFreeze,
md->invmass, md->ptype,
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):
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;
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff