gmx::ArrayRef<const rvec> invMassPerDim,
t_state* state,
const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
- const t_fcdata& fcdata,
+ t_fcdata* fcdata,
const gmx_ekindata_t* ekind,
const matrix M,
int UpdatePart,
gmx::ArrayRef<const unsigned short> cFREEZE_;
//! Group index for temperature coupling
gmx::ArrayRef<const unsigned short> cTC_;
+ //! Group index for accleration groups
+ gmx::ArrayRef<const unsigned short> cAcceleration_;
private:
//! stochastic dynamics struct
gmx::ArrayRef<const rvec> invMassPerDim,
t_state* state,
const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
- const t_fcdata& fcdata,
+ t_fcdata* fcdata,
const gmx_ekindata_t* ekind,
const matrix M,
int updatePart,
gmx::ArrayRef<const t_grp_tcstat> tcstat,
gmx::ArrayRef<const unsigned short> cTC,
const rvec pRVScaleMatrixDiagonal,
- const rvec* gmx_restrict x,
- rvec* gmx_restrict xprime,
- VelocityType* gmx_restrict v,
- const rvec* gmx_restrict f)
+ const rvec* gmx_restrict x,
+ rvec* gmx_restrict xprime,
+ VelocityType* gmx_restrict v,
+ const rvec* gmx_restrict f)
{
real lambdaGroup;
{
vNew -= dtPressureCouple * pRVScaleMatrixDiagonal[d] * v[a][d];
}
- // TODO: Remove NOLINTs once clang-tidy is updated to v11, it should be able to handle constexpr.
- if constexpr (storeUpdatedVelocities == StoreUpdatedVelocities::yes) // NOLINT // NOLINTNEXTLINE
+ if constexpr (storeUpdatedVelocities == StoreUpdatedVelocities::yes)
{
v[a][d] = vNew;
}
- xprime[a][d] = x[a][d] + vNew * dt; // NOLINT(readability-misleading-indentation)
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ xprime[a][d] = x[a][d] + vNew * dt;
}
}
}
real dt,
gmx::ArrayRef<const real> invMass,
gmx::ArrayRef<const t_grp_tcstat> tcstat,
- const rvec* gmx_restrict x,
- rvec* gmx_restrict xprime,
- VelocityType* gmx_restrict v,
- const rvec* gmx_restrict f)
+ const rvec* gmx_restrict x,
+ rvec* gmx_restrict xprime,
+ VelocityType* gmx_restrict v,
+ const rvec* gmx_restrict f)
{
SimdReal timestep(dt);
SimdReal lambdaSystem(tcstat[0].lambda);
v1 = fma(f1 * invMass1, timestep, lambdaSystem * v1);
v2 = fma(f2 * invMass2, timestep, lambdaSystem * v2);
- // TODO: Remove NOLINTs once clang-tidy is updated to v11, it should be able to handle constexpr.
- if constexpr (storeUpdatedVelocities == StoreUpdatedVelocities::yes) // NOLINT // NOLINTNEXTLINE
+ if constexpr (storeUpdatedVelocities == StoreUpdatedVelocities::yes)
{
simdStoreRvecs(v, a, v0, v1, v2);
}
- SimdReal x0, x1, x2; // NOLINT(readability-misleading-indentation)
+ // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
+ SimdReal x0, x1, x2;
simdLoadRvecs(x, a, &x0, &x1, &x2);
SimdReal xprime0 = fma(v0, timestep, x0);
enum class AccelerationType
{
none,
+ group,
cosine
};
* \param[in] dtPressureCouple Time step for pressure coupling, is 0 when no pressure
* coupling should be applied at this step.
* \param[in] cTC Temperature coupling group indices
+ * \param[in] cAcceleration Acceleration group indices
+ * \param[in] acceleration Acceleration per group.
* \param[in] invMassPerDim Inverse mass per dimension
* \param[in] ekind Kinetic energy data.
* \param[in] box The box dimensions.
real dt,
real dtPressureCouple,
gmx::ArrayRef<const unsigned short> cTC,
+ gmx::ArrayRef<const unsigned short> cAcceleration,
+ const rvec* gmx_restrict acceleration,
gmx::ArrayRef<const rvec> invMassPerDim,
const gmx_ekindata_t* ekind,
const matrix box,
- const rvec* gmx_restrict x,
- rvec* gmx_restrict xprime,
- rvec* gmx_restrict v,
- const rvec* gmx_restrict f,
- const double* gmx_restrict nh_vxi,
- const int nsttcouple,
- const matrix M)
+ const rvec* gmx_restrict x,
+ rvec* gmx_restrict xprime,
+ rvec* gmx_restrict v,
+ const rvec* gmx_restrict f,
+ const double* gmx_restrict nh_vxi,
+ const int nsttcouple,
+ const matrix M)
{
/* This is a version of the leap-frog integrator that supports
* all combinations of T-coupling, P-coupling and NEMD.
gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind->tcstat;
/* Initialize group values, changed later when multiple groups are used */
+ int ga = 0;
int gt = 0;
real omega_Z = 2 * static_cast<real>(M_PI) / box[ZZ][ZZ];
switch (accelerationType)
{
case AccelerationType::none: copy_rvec(v[n], vRel); break;
+ case AccelerationType::group:
+ if (!cAcceleration.empty())
+ {
+ ga = cAcceleration[n];
+ }
+ /* With constant acceleration we do scale the velocity of the accelerated groups */
+ copy_rvec(v[n], vRel);
+ break;
case AccelerationType::cosine:
cosineZ = std::cos(x[n][ZZ] * omega_Z);
vCosine = cosineZ * ekind->cosacc.vcos;
switch (accelerationType)
{
case AccelerationType::none: break;
+ case AccelerationType::group:
+ /* Apply the constant acceleration */
+ vNew += acceleration[ga][d] * dt;
+ break;
case AccelerationType::cosine:
if (d == XX)
{
}
/*! \brief Handles the Leap-frog MD x and v integration */
-static void do_update_md(int start,
- int nrend,
- real dt,
- int64_t step,
- const rvec* gmx_restrict x,
- rvec* gmx_restrict xprime,
- rvec* gmx_restrict v,
- const rvec* gmx_restrict f,
- const TemperatureCoupling etc,
- const PressureCoupling epc,
- const int nsttcouple,
- const int nstpcouple,
- gmx::ArrayRef<const unsigned short> cTC,
+static void do_update_md(int start,
+ int nrend,
+ real dt,
+ int64_t step,
+ const rvec* gmx_restrict x,
+ rvec* gmx_restrict xprime,
+ rvec* gmx_restrict v,
+ const rvec* gmx_restrict f,
+ const TemperatureCoupling etc,
+ const PressureCoupling epc,
+ const int nsttcouple,
+ const int nstpcouple,
+ gmx::ArrayRef<const unsigned short> cTC,
+ const bool useConstantAcceleration,
+ gmx::ArrayRef<const unsigned short> cAcceleration,
+ const rvec* acceleration,
gmx::ArrayRef<const real> gmx_unused invmass,
gmx::ArrayRef<const rvec> invMassPerDim,
const gmx_ekindata_t* ekind,
const matrix box,
- const double* gmx_restrict nh_vxi,
- const matrix M,
- bool gmx_unused havePartiallyFrozenAtoms)
+ const double* gmx_restrict nh_vxi,
+ const matrix M,
+ bool gmx_unused havePartiallyFrozenAtoms)
{
GMX_ASSERT(nrend == start || xprime != x,
"For SIMD optimization certain compilers need to have xprime != x");
real dtPressureCouple = (doParrinelloRahman ? nstpcouple * dt : 0);
- /* NEMD cosine acceleration is applied in updateMDLeapFrogGeneral */
- bool doAcceleration = (ekind->cosacc.cos_accel != 0);
+ /* NEMD (also cosine) acceleration is applied in updateMDLeapFrogGeneral */
+ const bool doAcceleration = (useConstantAcceleration || ekind->cosacc.cos_accel != 0);
if (doNoseHoover || doPROffDiagonal || doAcceleration)
{
dt,
dtPressureCouple,
cTC,
+ cAcceleration,
+ acceleration,
invMassPerDim,
ekind,
box,
nsttcouple,
stepM);
}
+ else if (useConstantAcceleration)
+ {
+ updateMDLeapfrogGeneral<AccelerationType::group>(start,
+ nrend,
+ doNoseHoover,
+ dt,
+ dtPressureCouple,
+ cTC,
+ cAcceleration,
+ acceleration,
+ invMassPerDim,
+ ekind,
+ box,
+ x,
+ xprime,
+ v,
+ f,
+ nh_vxi,
+ nsttcouple,
+ stepM);
+ }
else
{
updateMDLeapfrogGeneral<AccelerationType::cosine>(start,
dt,
dtPressureCouple,
cTC,
+ cAcceleration,
+ acceleration,
invMassPerDim,
ekind,
box,
}
}
/*! \brief Handles the Leap-frog MD x and v integration */
-static void doUpdateMDDoNotUpdateVelocities(int start,
- int nrend,
- real dt,
+static void doUpdateMDDoNotUpdateVelocities(int start,
+ int nrend,
+ real dt,
const rvec* gmx_restrict x,
- rvec* gmx_restrict xprime,
+ rvec* gmx_restrict xprime,
const rvec* gmx_restrict v,
const rvec* gmx_restrict f,
- bool gmx_unused havePartiallyFrozenAtoms,
+ bool gmx_unused havePartiallyFrozenAtoms,
gmx::ArrayRef<const real> gmx_unused invmass,
gmx::ArrayRef<const rvec> invMassPerDim,
const gmx_ekindata_t& ekind)
int nrend,
real dt,
gmx::ArrayRef<const ivec> nFreeze,
+ gmx::ArrayRef<const unsigned short> cAcceleration,
+ const rvec* acceleration,
gmx::ArrayRef<const real> invmass,
gmx::ArrayRef<const ParticleType> ptype,
gmx::ArrayRef<const unsigned short> cFREEZE,
real veta,
real alpha)
{
- int gf = 0;
+ int gf = 0, ga = 0;
int n, d;
real g, mv1, mv2;
{
gf = cFREEZE[n];
}
+ if (!cAcceleration.empty())
+ {
+ ga = cAcceleration[n];
+ }
for (d = 0; d < DIM; d++)
{
if ((ptype[n] != ParticleType::Shell) && !nFreeze[gf][d])
{
- v[n][d] = mv1 * (mv1 * v[n][d] + 0.5 * (w_dt * mv2 * f[n][d]));
+ v[n][d] = mv1 * (mv1 * v[n][d] + 0.5 * (w_dt * mv2 * f[n][d]))
+ + 0.5 * acceleration[ga][d] * dt;
}
else
{
}
Update::Impl::Impl(const t_inputrec& inputRecord, BoxDeformation* boxDeformation) :
- sd_(inputRecord),
- deform_(boxDeformation)
+ sd_(inputRecord), deform_(boxDeformation)
{
update_temperature_constants(inputRecord);
xp_.resizeWithPadding(0);
void Update::updateAfterPartition(int numAtoms,
gmx::ArrayRef<const unsigned short> cFREEZE,
- gmx::ArrayRef<const unsigned short> cTC)
+ gmx::ArrayRef<const unsigned short> cTC,
+ gmx::ArrayRef<const unsigned short> cAcceleration)
{
impl_->xp()->resizeWithPadding(numAtoms);
- impl_->cFREEZE_ = cFREEZE;
- impl_->cTC_ = cTC;
+ impl_->cFREEZE_ = cFREEZE;
+ impl_->cTC_ = cTC;
+ impl_->cAcceleration_ = cAcceleration;
}
/*! \brief Sets the SD update type */
gmx::ArrayRef<const ParticleType> ptype,
gmx::ArrayRef<const unsigned short> cFREEZE,
gmx::ArrayRef<const unsigned short> cTC,
+ gmx::ArrayRef<const unsigned short> cAcceleration,
+ const rvec* acceleration,
const rvec x[],
rvec xprime[],
rvec v[],
int seed,
const int* gatindex)
{
- // cTC and cFREEZE can be nullptr any time, but various
+ // cTC, cACC and cFREEZE can be nullptr any time, but various
// instantiations do not make sense with particular pointer
// values.
if (updateType == SDUpdate::ForcesOnly)
if (updateType == SDUpdate::FrictionAndNoiseOnly)
{
GMX_ASSERT(f == nullptr, "SD update with only noise cannot handle forces");
+ GMX_ASSERT(cAcceleration.empty(),
+ "SD update with only noise cannot handle acceleration groups");
}
if (updateType == SDUpdate::Combined)
{
real inverseMass = invmass[n];
real invsqrtMass = std::sqrt(inverseMass);
- int freezeGroup = !cFREEZE.empty() ? cFREEZE[n] : 0;
- int temperatureGroup = !cTC.empty() ? cTC[n] : 0;
+ int freezeGroup = !cFREEZE.empty() ? cFREEZE[n] : 0;
+ int accelerationGroup = !cAcceleration.empty() ? cAcceleration[n] : 0;
+ int temperatureGroup = !cTC.empty() ? cTC[n] : 0;
for (int d = 0; d < DIM; d++)
{
{
if (updateType == SDUpdate::ForcesOnly)
{
- real vn = v[n][d] + inverseMass * f[n][d] * dt;
+ real vn = v[n][d] + (inverseMass * f[n][d] + acceleration[accelerationGroup][d]) * dt;
v[n][d] = vn;
// Simple position update.
xprime[n][d] = x[n][d] + v[n][d] * dt;
}
else
{
- real vn = v[n][d] + inverseMass * f[n][d] * dt;
+ real vn = v[n][d] + (inverseMass * f[n][d] + acceleration[accelerationGroup][d]) * dt;
v[n][d] = (vn * sd.sdc[temperatureGroup].em
+ invsqrtMass * sd.sdsig[temperatureGroup].V * dist(rng));
// Here we include half of the friction+noise
}
}
-static void do_update_sd(int start,
- int nrend,
- real dt,
- int64_t step,
- const rvec* gmx_restrict x,
- rvec* gmx_restrict xprime,
- rvec* gmx_restrict v,
+static void do_update_sd(int start,
+ int nrend,
+ real dt,
+ int64_t step,
+ const rvec* gmx_restrict x,
+ rvec* gmx_restrict xprime,
+ rvec* gmx_restrict v,
const rvec* gmx_restrict f,
gmx::ArrayRef<const ivec> nFreeze,
gmx::ArrayRef<const real> invmass,
gmx::ArrayRef<const ParticleType> ptype,
gmx::ArrayRef<const unsigned short> cFREEZE,
gmx::ArrayRef<const unsigned short> cTC,
+ gmx::ArrayRef<const unsigned short> cAcceleration,
+ const rvec* acceleration,
int seed,
const t_commrec* cr,
const gmx_stochd_t& sd,
ptype,
cFREEZE,
gmx::ArrayRef<const unsigned short>(),
+ cAcceleration,
+ acceleration,
x,
xprime,
v,
}
else
{
- doSDUpdateGeneral<SDUpdate::Combined>(sd,
- start,
- nrend,
- dt,
- nFreeze,
- invmass,
- ptype,
- cFREEZE,
- cTC,
- x,
- xprime,
- v,
- f,
- step,
- seed,
- DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
+ doSDUpdateGeneral<SDUpdate::Combined>(
+ sd,
+ start,
+ nrend,
+ dt,
+ nFreeze,
+ invmass,
+ ptype,
+ cFREEZE,
+ cTC,
+ cAcceleration,
+ acceleration,
+ x,
+ xprime,
+ v,
+ f,
+ step,
+ seed,
+ haveDDAtomOrdering(*cr) ? cr->dd->globalAtomIndices.data() : nullptr);
}
}
-static void do_update_bd(int start,
- int nrend,
- real dt,
- int64_t step,
- const rvec* gmx_restrict x,
- rvec* gmx_restrict xprime,
- rvec* gmx_restrict v,
+static void do_update_bd(int start,
+ int nrend,
+ real dt,
+ int64_t step,
+ const rvec* gmx_restrict x,
+ rvec* gmx_restrict xprime,
+ rvec* gmx_restrict v,
const rvec* gmx_restrict f,
gmx::ArrayRef<const ivec> nFreeze,
gmx::ArrayRef<const real> invmass,
ptype,
cFREEZE_,
cTC_,
+ cAcceleration_,
+ inputRecord.opts.acceleration,
state->x.rvec_array(),
xp_.rvec_array(),
state->v.rvec_array(),
nullptr,
step,
inputRecord.ld_seed,
- DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
+ haveDDAtomOrdering(*cr) ? cr->dd->globalAtomIndices.data() : nullptr);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
gmx::ArrayRef<const rvec> invMassPerDim,
t_state* state,
const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
- const t_fcdata& fcdata,
+ t_fcdata* fcdata,
const gmx_ekindata_t* ekind,
const matrix M,
int updatePart,
/* We need to update the NMR restraint history when time averaging is used */
if (state->flags & enumValueToBitMask(StateEntry::DisreRm3Tav))
{
- update_disres_history(*fcdata.disres, &state->hist);
+ update_disres_history(*fcdata->disres, &state->hist);
}
if (state->flags & enumValueToBitMask(StateEntry::OrireDtav))
{
- update_orires_history(*fcdata.orires, &state->hist);
+ GMX_ASSERT(fcdata, "Need valid fcdata");
+ fcdata->orires->updateHistory();
}
/* ############# START The update of velocities and positions ######### */
inputRecord.nsttcouple,
inputRecord.nstpcouple,
cTC_,
+ inputRecord.useConstantAcceleration,
+ cAcceleration_,
+ inputRecord.opts.acceleration,
invMass,
invMassPerDim,
ekind,
ptype,
cFREEZE_,
cTC_,
+ cAcceleration_,
+ inputRecord.opts.acceleration,
inputRecord.ld_seed,
cr,
sd_,
inputRecord.bd_fric,
sd_.bd_rf.data(),
inputRecord.ld_seed,
- DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
+ haveDDAtomOrdering(*cr) ? cr->dd->globalAtomIndices.data() : nullptr);
break;
case (IntegrationAlgorithm::VV):
case (IntegrationAlgorithm::VVAK):
dt,
gmx::arrayRefFromArray(inputRecord.opts.nFreeze,
inputRecord.opts.ngfrz),
+ cAcceleration_,
+ inputRecord.opts.acceleration,
invMass,
ptype,
cFREEZE_,
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff