simulation, and afterward use position restraints in conjunction
with constant pressure.
+accelerate group
+
+ On each atom in an “accelerate group” an acceleration
+ :math:`\mathbf{a}^g` is imposed. This is equivalent to
+ a mass-weighted external force. This feature makes it possible to
+ drive the system into a non-equilibrium state to compute,
+ for example, transport properties.
+
energy-monitor group
Mutual interactions between all energy-monitor groups are compiled
Also, please use the syntax :issue:`number` to reference issues on GitLab, without the
a space between the colon and number!
+Removed velocity output for acceleration groups
+"""""""""""""""""""""""""""""""""""""""""""""""
+
+The reported velocity in the energy file for acceleration groups was always
+zero. Now their velocity is no longer reported in the energy file.
+
+:issue:`1354`
+
Use correct c0 parameter in Me2PO4 in OPLSAA
""""""""""""""""""""""""""""""""""""""""""""
-OPLSAA torsions must sum to 0, but the paramters for Me2PO4 did not do so. Changed the c0
+OPLSAA torsions must sum to 0, but the parameters for Me2PO4 did not do so. Changed the c0
parameter to the correct value.
:issue:`4075`
-
Removed functionality
^^^^^^^^^^^^^^^^^^^^^
-Removed constant-acceleration groups support
-""""""""""""""""""""""""""""""""""""""""""""
-This code has been broken since before GROMACS 4.6, so it has been
-removed.
-
-:issue:`1354`
-
.. Note to developers!
Please use """"""" to underline the individual entries for fixed issues in the subfolders,
otherwise the formatting on the webpage is messed up.
Non-equilibrium MD
^^^^^^^^^^^^^^^^^^
+.. mdp:: acc-grps
+
+ groups for constant acceleration (*e.g.* ``Protein Sol``) all atoms
+ in groups Protein and Sol will experience constant acceleration as
+ specified in the :mdp:`accelerate` line. Note that the kinetic energy
+ of the center of mass of accelarated groups contributes to the kinetic
+ energy and temperature of the system. If this is not desired, make
+ each accelerate group also a separate temperature coupling group.
+
+.. mdp:: accelerate
+
+ (0) [nm ps\ :sup:`-2`]
+ acceleration for :mdp:`acc-grps`; x, y and z for each group
+ (*e.g.* ``0.1 0.0 0.0 -0.1 0.0 0.0`` means that first group has
+ constant acceleration of 0.1 nm ps\ :sup:`-2` in X direction, second group
+ the opposite).
+
.. mdp:: freezegrps
Groups that are to be frozen (*i.e.* their X, Y, and/or Z position
tpxv_RemovedConstantAcceleration, /**< Removed support for constant acceleration NEMD. */
tpxv_TransformationPullCoord, /**< Support for transformation pull coordinates */
tpxv_SoftcoreGapsys, /**< Added gapsys softcore function */
+ tpxv_ReaddedConstantAcceleration, /**< Re-added support for constant acceleration NEMD. */
tpxv_Count /**< the total number of tpxv versions */
};
{
ir->opts.nhchainlength = 1;
}
- int removedOptsNgacc = 0;
- if (serializer->reading() && file_version < tpxv_RemovedConstantAcceleration)
+ if (serializer->reading() && file_version >= tpxv_RemovedConstantAcceleration
+ && file_version < tpxv_ReaddedConstantAcceleration)
{
- serializer->doInt(&removedOptsNgacc);
+ ir->opts.ngacc = 0;
+ }
+ else
+ {
+ serializer->doInt(&ir->opts.ngacc);
}
serializer->doInt(&ir->opts.ngfrz);
serializer->doInt(&ir->opts.ngener);
snew(ir->opts.anneal_temp, ir->opts.ngtc);
snew(ir->opts.tau_t, ir->opts.ngtc);
snew(ir->opts.nFreeze, ir->opts.ngfrz);
+ snew(ir->opts.acceleration, ir->opts.ngacc);
snew(ir->opts.egp_flags, ir->opts.ngener * ir->opts.ngener);
}
if (ir->opts.ngtc > 0)
{
serializer->doIvecArray(ir->opts.nFreeze, ir->opts.ngfrz);
}
- if (serializer->reading() && file_version < tpxv_RemovedConstantAcceleration && removedOptsNgacc > 0)
+ if (ir->opts.ngacc > 0)
{
- std::vector<gmx::RVec> dummy;
- dummy.resize(removedOptsNgacc);
- serializer->doRvecArray(reinterpret_cast<rvec*>(dummy.data()), removedOptsNgacc);
- ir->useConstantAcceleration = std::any_of(dummy.begin(), dummy.end(), [](const gmx::RVec& vec) {
- return vec[XX] != 0.0 || vec[YY] != 0.0 || vec[ZZ] != 0.0;
- });
+ serializer->doRvecArray(ir->opts.acceleration, ir->opts.ngacc);
}
- else
+ if (serializer->reading())
{
ir->useConstantAcceleration = false;
+ for (int g = 0; g < ir->opts.ngacc; g++)
+ {
+ if (norm2(ir->opts.acceleration[g]) != 0)
+ {
+ ir->useConstantAcceleration = true;
+ }
+ }
}
serializer->doIntArray(ir->opts.egp_flags, ir->opts.ngener * ir->opts.ngener);
}
do_groups(serializer, &mtop->groups, &(mtop->symtab));
- if (file_version < tpxv_RemovedConstantAcceleration)
- {
- mtop->groups.groups[SimulationAtomGroupType::AccelerationUnused].clear();
- mtop->groups.groupNumbers[SimulationAtomGroupType::AccelerationUnused].clear();
- }
mtop->haveMoleculeIndices = true;
struct gmx_inputrec_strings
{
- char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], freeze[STRLEN], frdim[STRLEN],
- energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
- couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
- wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
+ char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], accelerationGroups[STRLEN],
+ acceleration[STRLEN], freeze[STRLEN], frdim[STRLEN], energy[STRLEN], user1[STRLEN],
+ user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN], couple_moltype[STRLEN],
+ orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN], wall_atomtype[STRLEN],
+ wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, std::string> fep_lambda;
char lambda_weights[STRLEN];
std::vector<std::string> pullGroupNames;
}
}
+static void convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, rvec* outputs)
+{
+ int i = 0, d = 0;
+ for (const auto& input : inputs)
+ {
+ try
+ {
+ outputs[i][d] = gmx::fromString<real>(input);
+ }
+ catch (gmx::GromacsException&)
+ {
+ auto message = gmx::formatString(
+ "Invalid value for mdp option %s. %s should only consist of real numbers "
+ "separated by spaces.",
+ name,
+ name);
+ warning_error(wi, message);
+ }
+ ++d;
+ if (d == DIM)
+ {
+ d = 0;
+ ++i;
+ }
+ }
+}
+
static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
{
opts->wall_atomtype[0] = nullptr;
/* Non-equilibrium MD stuff */
printStringNewline(&inp, "Non-equilibrium MD stuff");
+ setStringEntry(&inp, "acc-grps", inputrecStrings->accelerationGroups, nullptr);
+ setStringEntry(&inp, "accelerate", inputrecStrings->acceleration, nullptr);
setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
*defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
mdModulesNotifiers.preProcessingNotifier_.notify(defaultIndexGroupsAndNames);
+ auto accelerations = gmx::splitString(inputrecStrings->acceleration);
+ auto accelerationGroupNames = gmx::splitString(inputrecStrings->accelerationGroups);
+ if (accelerationGroupNames.size() * DIM != accelerations.size())
+ {
+ gmx_fatal(FARGS,
+ "Invalid Acceleration input: %zu groups and %zu acc. values",
+ accelerationGroupNames.size(),
+ accelerations.size());
+ }
+ do_numbering(natoms,
+ groups,
+ accelerationGroupNames,
+ defaultIndexGroups,
+ gnames,
+ SimulationAtomGroupType::Acceleration,
+ restnm,
+ GroupCoverage::AllGenerateRest,
+ bVerbose,
+ wi);
+ nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
+ snew(ir->opts.acceleration, nr);
+ ir->opts.ngacc = nr;
+
+ convertRvecs(wi, accelerations, "accelerations", ir->opts.acceleration);
+
auto freezeDims = gmx::splitString(inputrecStrings->frdim);
auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
if (freezeDims.size() != DIM * freezeGroupNames.size())
char err_buf[STRLEN];
int i, m, c, nmol;
bool bCharge;
+ real * mgrp, mt;
gmx_mtop_atomloop_block_t aloopb;
char warn_buf[STRLEN];
"constant by hand.");
}
+ ir->useConstantAcceleration = false;
+ for (int i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
+ {
+ if (norm2(ir->opts.acceleration[i]) != 0)
+ {
+ ir->useConstantAcceleration = true;
+ }
+ }
+ if (ir->useConstantAcceleration)
+ {
+ gmx::RVec acceleration = { 0.0_real, 0.0_real, 0.0_real };
+ snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
+ for (const AtomProxy atomP : AtomRange(*sys))
+ {
+ const t_atom& local = atomP.atom();
+ int i = atomP.globalAtomNumber();
+ mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
+ }
+ mt = 0.0;
+ for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
+ {
+ for (m = 0; (m < DIM); m++)
+ {
+ acceleration[m] += ir->opts.acceleration[i][m] * mgrp[i];
+ }
+ mt += mgrp[i];
+ }
+ for (m = 0; (m < DIM); m++)
+ {
+ if (fabs(acceleration[m]) > 1e-6)
+ {
+ const char* dim[DIM] = { "X", "Y", "Z" };
+ fprintf(stderr,
+ "Net Acceleration in %s direction, will %s be corrected\n",
+ dim[m],
+ ir->nstcomm != 0 ? "" : "not");
+ if (ir->nstcomm != 0 && m < ndof_com(ir))
+ {
+ acceleration[m] /= mt;
+ for (i = 0;
+ (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration]));
+ i++)
+ {
+ ir->opts.acceleration[i][m] -= acceleration[m];
+ }
+ }
+ }
+ }
+ sfree(mgrp);
+ }
+
if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->sc_alpha != 0
&& !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
{
dh_hist_spacing = 0.1
; Non-equilibrium MD stuff
+acc-grps =
+accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
dh_hist_spacing = 0.1
; Non-equilibrium MD stuff
+acc-grps =
+accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
dh_hist_spacing = 0.1
; Non-equilibrium MD stuff
+acc-grps =
+accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
dh_hist_spacing = 0.1
; Non-equilibrium MD stuff
+acc-grps =
+accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
dh_hist_spacing = 0.1
; Non-equilibrium MD stuff
+acc-grps =
+accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
dh_hist_spacing = 0.1
; Non-equilibrium MD stuff
+acc-grps =
+accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
dh_hist_spacing = 0.1
; Non-equilibrium MD stuff
+acc-grps =
+accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
dh_hist_spacing = 0.1
; Non-equilibrium MD stuff
+acc-grps =
+accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
dh_hist_spacing = 0.1
; Non-equilibrium MD stuff
+acc-grps =
+accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
sfree(mdatoms_->ptype);
sfree(mdatoms_->cTC);
sfree(mdatoms_->cENER);
+ sfree(mdatoms_->cACC);
sfree(mdatoms_->cFREEZE);
sfree(mdatoms_->cVCM);
sfree(mdatoms_->cORF);
/* We always copy cTC with domain decomposition */
}
srenew(md->cENER, md->nalloc);
+ if (inputrec.useConstantAcceleration)
+ {
+ srenew(md->cACC, md->nalloc);
+ }
if (inputrecFrozenAtoms(&inputrec))
{
srenew(md->cFREEZE, md->nalloc);
md->cTC[i] = groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ag];
}
md->cENER[i] = getGroupType(groups, SimulationAtomGroupType::EnergyOutput, ag);
+ if (md->cACC)
+ {
+ md->cACC[i] = groups.groupNumbers[SimulationAtomGroupType::Acceleration][ag];
+ }
if (md->cVCM)
{
md->cVCM[i] = groups.groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][ag];
update_ = std::make_unique<Update>(inputRecord_, nullptr);
update_->updateAfterPartition(numAtoms,
gmx::ArrayRef<const unsigned short>(),
- gmx::arrayRefFromArray(mdAtoms_.cTC, mdAtoms_.nr));
+ gmx::arrayRefFromArray(mdAtoms_.cTC, mdAtoms_.nr),
+ gmx::ArrayRef<const unsigned short>());
doPressureCouple_ = (nstpcouple != 0);
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
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,
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)
{
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,
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,
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
{
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
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,
ptype,
cFREEZE,
cTC,
+ cAcceleration,
+ acceleration,
x,
xprime,
v,
ptype,
cFREEZE_,
cTC_,
+ cAcceleration_,
+ inputRecord.opts.acceleration,
state->x.rvec_array(),
xp_.rvec_array(),
state->v.rvec_array(),
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_,
dt,
gmx::arrayRefFromArray(inputRecord.opts.nFreeze,
inputRecord.opts.ngfrz),
+ cAcceleration_,
+ inputRecord.opts.acceleration,
invMass,
ptype,
cFREEZE_,
* \param[in] numAtoms Updated number of atoms.
* \param[in] cFREEZE Group index for freezing
* \param[in] cTC Group index for center of mass motion removal
+ * \param[in] cAcceleration Group index for constant acceleration groups
*/
void 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);
/*! \brief Perform numerical integration step.
*
md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr)
: gmx::ArrayRef<const unsigned short>(),
md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
- : gmx::ArrayRef<const unsigned short>());
+ : gmx::ArrayRef<const unsigned short>(),
+ md->cACC ? gmx::arrayRefFromArray(md->cACC, md->nr)
+ : gmx::ArrayRef<const unsigned short>());
fr->longRangeNonbondeds->updateAfterPartition(*md);
}
else
md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr)
: gmx::ArrayRef<const unsigned short>(),
md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
- : gmx::ArrayRef<const unsigned short>());
+ : gmx::ArrayRef<const unsigned short>(),
+ md->cACC ? gmx::arrayRefFromArray(md->cACC, md->nr)
+ : gmx::ArrayRef<const unsigned short>());
fr->longRangeNonbondeds->updateAfterPartition(*md);
}
md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr)
: gmx::ArrayRef<const unsigned short>(),
md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
- : gmx::ArrayRef<const unsigned short>());
+ : gmx::ArrayRef<const unsigned short>(),
+ md->cACC ? gmx::arrayRefFromArray(md->cACC, md->nr)
+ : gmx::ArrayRef<const unsigned short>());
fr->longRangeNonbondeds->updateAfterPartition(*md);
}
}
md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr)
: gmx::ArrayRef<const unsigned short>(),
md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
- : gmx::ArrayRef<const unsigned short>());
+ : gmx::ArrayRef<const unsigned short>(),
+ md->cACC ? gmx::arrayRefFromArray(md->cACC, md->nr)
+ : gmx::ArrayRef<const unsigned short>());
fr->longRangeNonbondeds->updateAfterPartition(*md);
}
GMX_RELEASE_ASSERT(inputrec != nullptr, "All ranks should have a valid inputrec now");
partialDeserializedTpr.reset(nullptr);
- GMX_RELEASE_ASSERT(
- !inputrec->useConstantAcceleration,
- "Linear acceleration has been removed in GROMACS 2022, and was broken for many years "
- "before that. Use GROMACS 4.5 or earlier if you need this feature.");
-
- // Now we decide whether to use the domain decomposition machinery.
- // Note that this does not necessarily imply actually using multiple domains.
// Now the number of ranks is known to all ranks, and each knows
// the inputrec read by the master rank. The ranks can now all run
// the task-deciding functions and will agree on the result
auto invmass = gmx::arrayRefFromArray(md.invmass, md.nr);
dt = ir->delta_t;
- /* Does NOT work with freeze groups (yet) */
+ /* Does NOT work with freeze or acceleration groups (yet) */
for (n = 0; n < end; n++)
{
w_dt = invmass[n] * dt;
sfree(ir->opts.anneal_time);
sfree(ir->opts.anneal_temp);
sfree(ir->opts.tau_t);
+ sfree(ir->opts.acceleration);
sfree(ir->opts.nFreeze);
sfree(ir->opts.egp_flags);
}
}
+ pr_indent(out, indent);
+ fprintf(out, "acc:\t");
+ for (i = 0; (i < opts->ngacc); i++)
+ {
+ for (m = 0; (m < DIM); m++)
+ {
+ fprintf(out, " %10g", opts->acceleration[i][m]);
+ }
+ }
+ fprintf(out, "\n");
+
pr_indent(out, indent);
fprintf(out, "nfreeze:");
for (i = 0; (i < opts->ngfrz); i++)
char buf1[256], buf2[256];
cmp_int(fp, "inputrec->grpopts.ngtc", -1, opt1->ngtc, opt2->ngtc);
+ cmp_int(fp, "inputrec->grpopts.ngacc", -1, opt1->ngacc, opt2->ngacc);
cmp_int(fp, "inputrec->grpopts.ngfrz", -1, opt1->ngfrz, opt2->ngfrz);
cmp_int(fp, "inputrec->grpopts.ngener", -1, opt1->ngener, opt2->ngener);
for (i = 0; (i < std::min(opt1->ngtc, opt2->ngtc)); i++)
}
}
}
+ for (i = 0; (i < std::min(opt1->ngacc, opt2->ngacc)); i++)
+ {
+ cmp_rvec(fp, "inputrec->grpopts.acceleration", i, opt1->acceleration[i], opt2->acceleration[i], ftol, abstol);
+ }
for (i = 0; (i < std::min(opt1->ngfrz, opt2->ngfrz)); i++)
{
cmp_ivec(fp, "inputrec->grpopts.nFreeze", i, opt1->nFreeze[i], opt2->nFreeze[i]);
int ngtc = 0;
//! Number of of Nose-Hoover chains per group
int nhchainlength = 0;
+ //! Number of Accelerate groups
+ int ngacc;
//! Number of Freeze groups
int ngfrz = 0;
//! Number of Energy groups
real** anneal_temp = nullptr;
//! Tau coupling time
real* tau_t = nullptr;
+ //! Acceleration per group
+ rvec* acceleration = nullptr;
//! Whether the group will be frozen in each direction
ivec* nFreeze = nullptr;
//! Exclusions/tables of energy group pairs
bool bAdress = false;
//! Whether twin-range scheme is active - always false if a valid .tpr was read
bool useTwinRange = false;
- //! Whether we have constant acceleration - removed in GROMACS 2022
+ //! Whether we have constant acceleration
bool useConstantAcceleration = false;
//! KVT object that contains input parameters converted to the new style.
unsigned short* cTC;
//! Group index for energy matrix
unsigned short* cENER;
+ //! Group index for acceleration
+ unsigned short* cACC;
//! Group index for freezing
unsigned short* cFREEZE;
//! Group index for center of mass motion removal
&& conditionalAssert(!doRerun, "Rerun is not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
- && conditionalAssert(inputrec->cos_accel == 0.0,
+ && conditionalAssert(!inputrec->useConstantAcceleration && inputrec->cos_accel == 0.0,
"Acceleration is not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
"Only Parrinello-Rahman, Berendsen, and C-rescale pressure coupling are "
"supported.\n";
}
+ if (inputrec.cos_accel != 0 || inputrec.useConstantAcceleration)
+ {
+ errorMessage += "Acceleration is not supported.\n";
+ }
if (EEL_PME_EWALD(inputrec.coulombtype) && inputrec.epsilon_surface != 0)
{
// The graph is needed, but not supported
{
for (auto group : keysOf(gcount))
{
- if (group == SimulationAtomGroupType::AccelerationUnused)
- {
- continue;
- }
gcount[group][getGroupType(groups, group, i)]++;
}
}
{
TemperatureCoupling,
EnergyOutput,
- AccelerationUnused,
+ Acceleration,
Freeze,
User1,
User2,
swapcoords.cpp
tabulated_bonded_interactions.cpp
freezegroups.cpp
+ constantacceleration.cpp
# pseudo-library for code for mdrun
$<TARGET_OBJECTS:mdrun_objlib>
)
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ * \brief End-to-end tests checking sanity of results of simulations
+ * containing constant acceleration groups
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_mdrun_integration_tests
+ */
+#include "gmxpre.h"
+
+#include "config.h"
+
+#include "gromacs/topology/ifunc.h"
+#include "gromacs/utility/stringutil.h"
+
+#include "testutils/mpitest.h"
+#include "testutils/simulationdatabase.h"
+#include "testutils/testmatchers.h"
+
+#include "moduletest.h"
+#include "simulatorcomparison.h"
+#include "trajectoryreader.h"
+
+namespace gmx::test
+{
+namespace
+{
+/*! \brief Test fixture checking the velocities of aomts
+ *
+ * This tests that velocities of non-accelerated atoms are zero
+ * and that velocities of accelarated atoms match acceleration*time.
+ * This code assumes the first half the atoms are non-accelerated
+ * and the second half are accelerated
+ */
+using AccelerationGroupTestParams = std::tuple<std::string, std::string>;
+class AccelerationGroupTest :
+ public MdrunTestFixture,
+ public ::testing::WithParamInterface<AccelerationGroupTestParams>
+{
+public:
+ //! Check that the velocities are zero or accelerated
+ static void checkVelocities(const std::string& trajectoryName,
+ const RVec acceleration,
+ const FloatingPointTolerance& tolerance)
+ {
+ const size_t c_groupSize = 3;
+
+ const std::vector<RVec> zeroVelocities(c_groupSize, RVec{ 0, 0, 0 });
+
+ TrajectoryFrameReader trajectoryFrameReader(trajectoryName);
+ while (trajectoryFrameReader.readNextFrame())
+ {
+ const auto frame = trajectoryFrameReader.frame();
+ GMX_RELEASE_ASSERT(frame.v().size() == 2 * c_groupSize,
+ "Expect velocities for both atom groups");
+
+ const RVec referenceVelocity = real(frame.time()) * acceleration;
+ const std::vector<RVec> referenceVelocities(c_groupSize, referenceVelocity);
+
+ SCOPED_TRACE("Checking velocities of non-accelerated atoms");
+ ArrayRef<const RVec> nonAcceleratedVelocities = frame.v().subArray(0, c_groupSize);
+ EXPECT_THAT(nonAcceleratedVelocities, Pointwise(RVecEq(tolerance), zeroVelocities));
+
+ SCOPED_TRACE("Checking velocities of accelerated atoms");
+ ArrayRef<const RVec> acceleratedVelocities = frame.v().subArray(c_groupSize, c_groupSize);
+ EXPECT_THAT(acceleratedVelocities, Pointwise(RVecEq(tolerance), referenceVelocities));
+ }
+ }
+};
+
+TEST_P(AccelerationGroupTest, WithinTolerances)
+{
+ const auto& params = GetParam();
+ const auto& integrator = std::get<0>(params);
+ const auto& tcoupling = std::get<1>(params);
+ const auto& simulationName = "spc2";
+
+ // Prepare mdp input
+ auto mdpFieldValues = prepareMdpFieldValues(simulationName, integrator, tcoupling, "no");
+ mdpFieldValues["nsteps"] = "8";
+ mdpFieldValues["dt"] = "0.002";
+ mdpFieldValues["nstxout"] = "0";
+ mdpFieldValues["nstvout"] = "8";
+ mdpFieldValues["nstfout"] = "0";
+ mdpFieldValues["comm-mode"] = "none";
+ // The two groups will not see each other when the cut-off is 0.9 nm
+ mdpFieldValues["coulombtype"] = "reaction-field";
+ mdpFieldValues["rcoulomb"] = "0.8";
+ mdpFieldValues["rvdw"] = "0.8";
+ mdpFieldValues["verlet-buffer-tolerance"] = "-1";
+ mdpFieldValues["rlist"] = "0.9";
+ // Couple the (non-)accelerated groups separately, so their velocties are independent
+ mdpFieldValues["tc-grps"] = "firstWaterMolecule secondWaterMolecule";
+ mdpFieldValues["ref-t"] = "0.001 0.001";
+ // Use weak coupling so we can check vecolities of atoms with a tight tolerance
+ mdpFieldValues["tau-t"] = "10.0 10.0";
+ const RVec c_acceleration = { 2.0, 3.0, 1.5 };
+ mdpFieldValues["acc-grps"] = "secondWaterMolecule";
+ mdpFieldValues["accelerate"] = "2.0 3.0 1.5";
+ // Set initial velocities to zero
+ mdpFieldValues["gen-vel"] = "yes";
+ mdpFieldValues["gen-temp"] = "0";
+
+ // Run grompp
+ runner_.useTopGroAndNdxFromDatabase(simulationName);
+ runner_.useStringAsMdpFile(prepareMdpFileContents(mdpFieldValues));
+ runGrompp(&runner_);
+ // Run mdrun
+ runMdrun(&runner_);
+
+ // T-coupling causes changes in the velocities up to 1e-4
+ auto tolerance = absoluteTolerance((GMX_DOUBLE && tcoupling == "no") ? 1e-10 : 1e-4);
+
+ // Check the velocities of the non-accelerated and accelerated groups
+ checkVelocities(runner_.fullPrecisionTrajectoryFileName_, c_acceleration, tolerance);
+}
+
+// The v-rescale case check that Ekin computation and temperature coupling
+// can be performed independently for atoms groups, so the accelerations
+// are not affected. This can be useful in practice.
+INSTANTIATE_TEST_SUITE_P(AccelerationWorks,
+ AccelerationGroupTest,
+ ::testing::Combine(::testing::Values("md", "md-vv"),
+ ::testing::Values("no", "v-rescale")));
+} // namespace
+} // namespace gmx::test
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, 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.
{
return true;
}
- *listener->stream() << " Actual value: " << value2 << std::endl
- << "Expected value: " << value1 << std::endl
- << " Difference: " << diff.toString() << std::endl
- << " Tolerance: " << tolerance_.toString(diff);
+ if (listener->stream())
+ {
+ *listener->stream() << " Actual value: " << value2 << std::endl
+ << "Expected value: " << value1 << std::endl
+ << " Difference: " << diff.toString() << std::endl
+ << " Tolerance: " << tolerance_.toString(diff);
+ }
+
return false;
}
//! Describe to a human what matching means.
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff