From: Berk Hess Date: Fri, 29 Oct 2021 14:50:28 +0000 (+0000) Subject: Reimplement constant acceleration groups X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?p=alexxy%2Fgromacs.git;a=commitdiff_plain;h=f3f4224ea504eef7ba021f0f37fa9d2fc532505f Reimplement constant acceleration groups --- diff --git a/docs/reference-manual/algorithms/group-concept.rst b/docs/reference-manual/algorithms/group-concept.rst index 1c97ea602c..f62eae61eb 100644 --- a/docs/reference-manual/algorithms/group-concept.rst +++ b/docs/reference-manual/algorithms/group-concept.rst @@ -40,6 +40,14 @@ freeze group 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 diff --git a/docs/release-notes/2022/major/bugs-fixed.rst b/docs/release-notes/2022/major/bugs-fixed.rst index c41151c576..9be30ced09 100644 --- a/docs/release-notes/2022/major/bugs-fixed.rst +++ b/docs/release-notes/2022/major/bugs-fixed.rst @@ -26,11 +26,18 @@ both dhdl.xvg and the energy file, which is used by e.g. gmx bar, was correct. 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` - diff --git a/docs/release-notes/2022/major/removed-functionality.rst b/docs/release-notes/2022/major/removed-functionality.rst index 440bb6413e..f21e37f0c8 100644 --- a/docs/release-notes/2022/major/removed-functionality.rst +++ b/docs/release-notes/2022/major/removed-functionality.rst @@ -1,13 +1,6 @@ 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. diff --git a/docs/user-guide/mdp-options.rst b/docs/user-guide/mdp-options.rst index d0e990632a..a1b96fa84a 100644 --- a/docs/user-guide/mdp-options.rst +++ b/docs/user-guide/mdp-options.rst @@ -3023,6 +3023,23 @@ Expanded Ensemble calculations 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 diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index c71bc38f2d..d190d663ed 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -139,6 +139,7 @@ enum tpxv 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 */ }; @@ -1579,10 +1580,14 @@ static void do_inputrec(gmx::ISerializer* serializer, t_inputrec* ir, int file_v { 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); @@ -1597,6 +1602,7 @@ static void do_inputrec(gmx::ISerializer* serializer, t_inputrec* ir, int file_v 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) @@ -1609,18 +1615,20 @@ static void do_inputrec(gmx::ISerializer* serializer, t_inputrec* ir, int file_v { serializer->doIvecArray(ir->opts.nFreeze, ir->opts.ngfrz); } - if (serializer->reading() && file_version < tpxv_RemovedConstantAcceleration && removedOptsNgacc > 0) + if (ir->opts.ngacc > 0) { - std::vector dummy; - dummy.resize(removedOptsNgacc); - serializer->doRvecArray(reinterpret_cast(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); @@ -2632,11 +2640,6 @@ static void do_mtop(gmx::ISerializer* serializer, gmx_mtop_t* mtop, int file_ver } 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; diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp index ff47ab5a35..557a361c8f 100644 --- a/src/gromacs/gmxpreprocess/readir.cpp +++ b/src/gromacs/gmxpreprocess/readir.cpp @@ -107,10 +107,11 @@ using gmx::BasicVector; 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 fep_lambda; char lambda_weights[STRLEN]; std::vector pullGroupNames; @@ -1771,6 +1772,33 @@ static void convertReals(warninp_t wi, gmx::ArrayRef inputs, } } +static void convertRvecs(warninp_t wi, gmx::ArrayRef inputs, const char* name, rvec* outputs) +{ + int i = 0, d = 0; + for (const auto& input : inputs) + { + try + { + outputs[i][d] = gmx::fromString(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; @@ -2340,6 +2368,8 @@ void get_ir(const char* mdparin, /* 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); @@ -3907,6 +3937,31 @@ void do_index(const char* mdparin, *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()) @@ -4370,6 +4425,7 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_ char err_buf[STRLEN]; int i, m, c, nmol; bool bCharge; + real * mgrp, mt; gmx_mtop_atomloop_block_t aloopb; char warn_buf[STRLEN]; @@ -4571,6 +4627,57 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_ "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)) { diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsDefineParametersWithValuesIncludingAssignment.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsDefineParametersWithValuesIncludingAssignment.xml index a7cc35f439..0cf43f926c 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsDefineParametersWithValuesIncludingAssignment.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsDefineParametersWithValuesIncludingAssignment.xml @@ -264,6 +264,8 @@ dh_hist_size = 0 dh_hist_spacing = 0.1 ; Non-equilibrium MD stuff +acc-grps = +accelerate = freezegrps = freezedim = cos-acceleration = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricField.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricField.xml index cbf5084ddc..184ec9c259 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricField.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricField.xml @@ -264,6 +264,8 @@ dh_hist_size = 0 dh_hist_spacing = 0.1 ; Non-equilibrium MD stuff +acc-grps = +accelerate = freezegrps = freezedim = cos-acceleration = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldOscillating.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldOscillating.xml index 7e439fa8b5..53a4f45999 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldOscillating.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldOscillating.xml @@ -264,6 +264,8 @@ dh_hist_size = 0 dh_hist_spacing = 0.1 ; Non-equilibrium MD stuff +acc-grps = +accelerate = freezegrps = freezedim = cos-acceleration = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldPulsed.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldPulsed.xml index 71826f94aa..9b6639cc64 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldPulsed.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldPulsed.xml @@ -264,6 +264,8 @@ dh_hist_size = 0 dh_hist_spacing = 0.1 ; Non-equilibrium MD stuff +acc-grps = +accelerate = freezegrps = freezedim = cos-acceleration = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsEmptyLines.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsEmptyLines.xml index d6fc03c7cb..2298676e16 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsEmptyLines.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsEmptyLines.xml @@ -264,6 +264,8 @@ dh_hist_size = 0 dh_hist_spacing = 0.1 ; Non-equilibrium MD stuff +acc-grps = +accelerate = freezegrps = freezedim = cos-acceleration = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsImplicitSolventNo.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsImplicitSolventNo.xml index d6fc03c7cb..2298676e16 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsImplicitSolventNo.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsImplicitSolventNo.xml @@ -264,6 +264,8 @@ dh_hist_size = 0 dh_hist_spacing = 0.1 ; Non-equilibrium MD stuff +acc-grps = +accelerate = freezegrps = freezedim = cos-acceleration = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsKeyWithoutValue.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsKeyWithoutValue.xml index d6fc03c7cb..2298676e16 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsKeyWithoutValue.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsKeyWithoutValue.xml @@ -264,6 +264,8 @@ dh_hist_size = 0 dh_hist_spacing = 0.1 ; Non-equilibrium MD stuff +acc-grps = +accelerate = freezegrps = freezedim = cos-acceleration = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsMimic.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsMimic.xml index 391a10118c..38dfb294cb 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsMimic.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsMimic.xml @@ -264,6 +264,8 @@ dh_hist_size = 0 dh_hist_spacing = 0.1 ; Non-equilibrium MD stuff +acc-grps = +accelerate = freezegrps = freezedim = cos-acceleration = 0 diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesDifferentKindsOfMdpLines.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesDifferentKindsOfMdpLines.xml index 202eac0de7..f1a0281dc9 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesDifferentKindsOfMdpLines.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesDifferentKindsOfMdpLines.xml @@ -264,6 +264,8 @@ dh_hist_size = 0 dh_hist_spacing = 0.1 ; Non-equilibrium MD stuff +acc-grps = +accelerate = freezegrps = freezedim = cos-acceleration = 0 diff --git a/src/gromacs/mdlib/mdatoms.cpp b/src/gromacs/mdlib/mdatoms.cpp index 3749019a9d..e4ee7b7b57 100644 --- a/src/gromacs/mdlib/mdatoms.cpp +++ b/src/gromacs/mdlib/mdatoms.cpp @@ -87,6 +87,7 @@ MDAtoms::~MDAtoms() sfree(mdatoms_->ptype); sfree(mdatoms_->cTC); sfree(mdatoms_->cENER); + sfree(mdatoms_->cACC); sfree(mdatoms_->cFREEZE); sfree(mdatoms_->cVCM); sfree(mdatoms_->cORF); @@ -288,6 +289,10 @@ void atoms2md(const gmx_mtop_t& mtop, /* 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); @@ -474,6 +479,10 @@ void atoms2md(const gmx_mtop_t& mtop, 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]; diff --git a/src/gromacs/mdlib/tests/leapfrogtestdata.cpp b/src/gromacs/mdlib/tests/leapfrogtestdata.cpp index 3f66855450..cf80cda412 100644 --- a/src/gromacs/mdlib/tests/leapfrogtestdata.cpp +++ b/src/gromacs/mdlib/tests/leapfrogtestdata.cpp @@ -175,7 +175,8 @@ LeapFrogTestData::LeapFrogTestData(int numAtoms, update_ = std::make_unique(inputRecord_, nullptr); update_->updateAfterPartition(numAtoms, gmx::ArrayRef(), - gmx::arrayRefFromArray(mdAtoms_.cTC, mdAtoms_.nr)); + gmx::arrayRefFromArray(mdAtoms_.cTC, mdAtoms_.nr), + gmx::ArrayRef()); doPressureCouple_ = (nstpcouple != 0); diff --git a/src/gromacs/mdlib/update.cpp b/src/gromacs/mdlib/update.cpp index 800d68aa1d..3a7c4b0961 100644 --- a/src/gromacs/mdlib/update.cpp +++ b/src/gromacs/mdlib/update.cpp @@ -174,6 +174,8 @@ public: gmx::ArrayRef cFREEZE_; //! Group index for temperature coupling gmx::ArrayRef cTC_; + //! Group index for accleration groups + gmx::ArrayRef cAcceleration_; private: //! stochastic dynamics struct @@ -517,6 +519,7 @@ updateMDLeapfrogSimpleSimd(int start, enum class AccelerationType { none, + group, cosine }; @@ -530,6 +533,8 @@ enum class AccelerationType * \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. @@ -548,6 +553,8 @@ static void updateMDLeapfrogGeneral(int start, real dt, real dtPressureCouple, gmx::ArrayRef cTC, + gmx::ArrayRef cAcceleration, + const rvec* gmx_restrict acceleration, gmx::ArrayRef invMassPerDim, const gmx_ekindata_t* ekind, const matrix box, @@ -568,6 +575,7 @@ static void updateMDLeapfrogGeneral(int start, gmx::ArrayRef 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(M_PI) / box[ZZ][ZZ]; @@ -585,6 +593,14 @@ static void updateMDLeapfrogGeneral(int start, 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; @@ -613,6 +629,10 @@ static void updateMDLeapfrogGeneral(int start, 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) { @@ -641,6 +661,9 @@ static void do_update_md(int start, const int nsttcouple, const int nstpcouple, gmx::ArrayRef cTC, + const bool useConstantAcceleration, + gmx::ArrayRef cAcceleration, + const rvec* acceleration, gmx::ArrayRef gmx_unused invmass, gmx::ArrayRef invMassPerDim, const gmx_ekindata_t* ekind, @@ -662,8 +685,8 @@ static void do_update_md(int start, 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) { @@ -686,6 +709,8 @@ static void do_update_md(int start, dt, dtPressureCouple, cTC, + cAcceleration, + acceleration, invMassPerDim, ekind, box, @@ -697,6 +722,27 @@ static void do_update_md(int start, nsttcouple, stepM); } + else if (useConstantAcceleration) + { + updateMDLeapfrogGeneral(start, + nrend, + doNoseHoover, + dt, + dtPressureCouple, + cTC, + cAcceleration, + acceleration, + invMassPerDim, + ekind, + box, + x, + xprime, + v, + f, + nh_vxi, + nsttcouple, + stepM); + } else { updateMDLeapfrogGeneral(start, @@ -705,6 +751,8 @@ static void do_update_md(int start, dt, dtPressureCouple, cTC, + cAcceleration, + acceleration, invMassPerDim, ekind, box, @@ -820,6 +868,8 @@ static void do_update_vv_vel(int start, int nrend, real dt, gmx::ArrayRef nFreeze, + gmx::ArrayRef cAcceleration, + const rvec* acceleration, gmx::ArrayRef invmass, gmx::ArrayRef ptype, gmx::ArrayRef cFREEZE, @@ -829,7 +879,7 @@ static void do_update_vv_vel(int start, real veta, real alpha) { - int gf = 0; + int gf = 0, ga = 0; int n, d; real g, mv1, mv2; @@ -851,12 +901,17 @@ static void do_update_vv_vel(int start, { 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 { @@ -1009,12 +1064,14 @@ Update::Impl::Impl(const t_inputrec& inputRecord, BoxDeformation* boxDeformation void Update::updateAfterPartition(int numAtoms, gmx::ArrayRef cFREEZE, - gmx::ArrayRef cTC) + gmx::ArrayRef cTC, + gmx::ArrayRef 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 */ @@ -1048,6 +1105,8 @@ static void doSDUpdateGeneral(const gmx_stochd_t& sd, gmx::ArrayRef ptype, gmx::ArrayRef cFREEZE, gmx::ArrayRef cTC, + gmx::ArrayRef cAcceleration, + const rvec* acceleration, const rvec x[], rvec xprime[], rvec v[], @@ -1056,7 +1115,7 @@ static void doSDUpdateGeneral(const gmx_stochd_t& sd, 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) @@ -1067,6 +1126,8 @@ static void doSDUpdateGeneral(const gmx_stochd_t& sd, 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) { @@ -1086,8 +1147,9 @@ static void doSDUpdateGeneral(const gmx_stochd_t& sd, 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++) { @@ -1095,7 +1157,7 @@ static void doSDUpdateGeneral(const gmx_stochd_t& sd, { 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; @@ -1112,7 +1174,7 @@ static void doSDUpdateGeneral(const gmx_stochd_t& sd, } 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 @@ -1149,6 +1211,8 @@ static void do_update_sd(int start, gmx::ArrayRef ptype, gmx::ArrayRef cFREEZE, gmx::ArrayRef cTC, + gmx::ArrayRef cAcceleration, + const rvec* acceleration, int seed, const t_commrec* cr, const gmx_stochd_t& sd, @@ -1166,6 +1230,8 @@ static void do_update_sd(int start, ptype, cFREEZE, gmx::ArrayRef(), + cAcceleration, + acceleration, x, xprime, v, @@ -1186,6 +1252,8 @@ static void do_update_sd(int start, ptype, cFREEZE, cTC, + cAcceleration, + acceleration, x, xprime, v, @@ -1426,6 +1494,8 @@ void Update::Impl::update_sd_second_half(const t_inputrec& input ptype, cFREEZE_, cTC_, + cAcceleration_, + inputRecord.opts.acceleration, state->x.rvec_array(), xp_.rvec_array(), state->v.rvec_array(), @@ -1582,6 +1652,9 @@ void Update::Impl::update_coords(const t_inputrec& inputRecord, inputRecord.nsttcouple, inputRecord.nstpcouple, cTC_, + inputRecord.useConstantAcceleration, + cAcceleration_, + inputRecord.opts.acceleration, invMass, invMassPerDim, ekind, @@ -1604,6 +1677,8 @@ void Update::Impl::update_coords(const t_inputrec& inputRecord, ptype, cFREEZE_, cTC_, + cAcceleration_, + inputRecord.opts.acceleration, inputRecord.ld_seed, cr, sd_, @@ -1646,6 +1721,8 @@ void Update::Impl::update_coords(const t_inputrec& inputRecord, dt, gmx::arrayRefFromArray(inputRecord.opts.nFreeze, inputRecord.opts.ngfrz), + cAcceleration_, + inputRecord.opts.acceleration, invMass, ptype, cFREEZE_, diff --git a/src/gromacs/mdlib/update.h b/src/gromacs/mdlib/update.h index a3c4844e2e..e8f95bc24c 100644 --- a/src/gromacs/mdlib/update.h +++ b/src/gromacs/mdlib/update.h @@ -97,10 +97,12 @@ public: * \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 cFREEZE, - gmx::ArrayRef cTC); + gmx::ArrayRef cTC, + gmx::ArrayRef cAcceleration); /*! \brief Perform numerical integration step. * diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 3898e07d3a..735b5daa06 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -395,7 +395,9 @@ void gmx::LegacySimulator::do_md() md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr) : gmx::ArrayRef(), md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr) - : gmx::ArrayRef()); + : gmx::ArrayRef(), + md->cACC ? gmx::arrayRefFromArray(md->cACC, md->nr) + : gmx::ArrayRef()); fr->longRangeNonbondeds->updateAfterPartition(*md); } else @@ -409,7 +411,9 @@ void gmx::LegacySimulator::do_md() md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr) : gmx::ArrayRef(), md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr) - : gmx::ArrayRef()); + : gmx::ArrayRef(), + md->cACC ? gmx::arrayRefFromArray(md->cACC, md->nr) + : gmx::ArrayRef()); fr->longRangeNonbondeds->updateAfterPartition(*md); } @@ -1004,7 +1008,9 @@ void gmx::LegacySimulator::do_md() md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr) : gmx::ArrayRef(), md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr) - : gmx::ArrayRef()); + : gmx::ArrayRef(), + md->cACC ? gmx::arrayRefFromArray(md->cACC, md->nr) + : gmx::ArrayRef()); fr->longRangeNonbondeds->updateAfterPartition(*md); } } @@ -1968,7 +1974,9 @@ void gmx::LegacySimulator::do_md() md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr) : gmx::ArrayRef(), md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr) - : gmx::ArrayRef()); + : gmx::ArrayRef(), + md->cACC ? gmx::arrayRefFromArray(md->cACC, md->nr) + : gmx::ArrayRef()); fr->longRangeNonbondeds->updateAfterPartition(*md); } diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp index 43551e0fdf..313b4ebb34 100644 --- a/src/gromacs/mdrun/runner.cpp +++ b/src/gromacs/mdrun/runner.cpp @@ -970,13 +970,6 @@ int Mdrunner::mdrunner() 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 diff --git a/src/gromacs/mdrun/shellfc.cpp b/src/gromacs/mdrun/shellfc.cpp index 2f4c2d6e93..7657596888 100644 --- a/src/gromacs/mdrun/shellfc.cpp +++ b/src/gromacs/mdrun/shellfc.cpp @@ -852,7 +852,7 @@ static void init_adir(gmx_shellfc_t* shfc, 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; diff --git a/src/gromacs/mdtypes/inputrec.cpp b/src/gromacs/mdtypes/inputrec.cpp index 9b51dd8ced..d93e399c4f 100644 --- a/src/gromacs/mdtypes/inputrec.cpp +++ b/src/gromacs/mdtypes/inputrec.cpp @@ -330,6 +330,7 @@ void done_inputrec(t_inputrec* ir) 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); @@ -405,6 +406,17 @@ static void pr_grp_opts(FILE* out, int indent, const char* title, const t_grpopt } } + 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++) @@ -1077,6 +1089,7 @@ static void cmp_grpopts(FILE* fp, const t_grpopts* opt1, const t_grpopts* opt2, 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++) @@ -1108,6 +1121,10 @@ static void cmp_grpopts(FILE* fp, const t_grpopts* opt1, const t_grpopts* opt2, } } } + 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]); diff --git a/src/gromacs/mdtypes/inputrec.h b/src/gromacs/mdtypes/inputrec.h index a7e8e9c4e4..4c01f0b93b 100644 --- a/src/gromacs/mdtypes/inputrec.h +++ b/src/gromacs/mdtypes/inputrec.h @@ -69,6 +69,8 @@ struct t_grpopts 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 @@ -87,6 +89,8 @@ struct t_grpopts 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 @@ -566,7 +570,7 @@ struct t_inputrec // NOLINT (clang-analyzer-optin.performance.Padding) 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. diff --git a/src/gromacs/mdtypes/mdatom.h b/src/gromacs/mdtypes/mdatom.h index 91ee9f00d7..c80ffdeaed 100644 --- a/src/gromacs/mdtypes/mdatom.h +++ b/src/gromacs/mdtypes/mdatom.h @@ -117,6 +117,8 @@ typedef struct t_mdatoms 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 diff --git a/src/gromacs/modularsimulator/modularsimulator.cpp b/src/gromacs/modularsimulator/modularsimulator.cpp index fac5ea0d86..03e64cc6e3 100644 --- a/src/gromacs/modularsimulator/modularsimulator.cpp +++ b/src/gromacs/modularsimulator/modularsimulator.cpp @@ -401,7 +401,7 @@ bool ModularSimulator::isInputCompatible(bool exitOn && 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 diff --git a/src/gromacs/taskassignment/decidegpuusage.cpp b/src/gromacs/taskassignment/decidegpuusage.cpp index 5bac5adcef..1c9b27b729 100644 --- a/src/gromacs/taskassignment/decidegpuusage.cpp +++ b/src/gromacs/taskassignment/decidegpuusage.cpp @@ -637,6 +637,10 @@ bool decideWhetherToUseGpuForUpdate(const bool isDomainDecom "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 diff --git a/src/gromacs/tools/dump.cpp b/src/gromacs/tools/dump.cpp index 541d6da1bf..d6506af7f9 100644 --- a/src/gromacs/tools/dump.cpp +++ b/src/gromacs/tools/dump.cpp @@ -168,10 +168,6 @@ void list_tpr(const char* fn, { for (auto group : keysOf(gcount)) { - if (group == SimulationAtomGroupType::AccelerationUnused) - { - continue; - } gcount[group][getGroupType(groups, group, i)]++; } } diff --git a/src/gromacs/topology/topology.h b/src/gromacs/topology/topology.h index ed809ab418..9278a2a864 100644 --- a/src/gromacs/topology/topology.h +++ b/src/gromacs/topology/topology.h @@ -56,7 +56,7 @@ enum class SimulationAtomGroupType : int { TemperatureCoupling, EnergyOutput, - AccelerationUnused, + Acceleration, Freeze, User1, User2, diff --git a/src/programs/mdrun/tests/CMakeLists.txt b/src/programs/mdrun/tests/CMakeLists.txt index 74c4b9c356..086dc308ca 100644 --- a/src/programs/mdrun/tests/CMakeLists.txt +++ b/src/programs/mdrun/tests/CMakeLists.txt @@ -111,6 +111,7 @@ gmx_add_gtest_executable(${exename} swapcoords.cpp tabulated_bonded_interactions.cpp freezegroups.cpp + constantacceleration.cpp # pseudo-library for code for mdrun $ ) diff --git a/src/programs/mdrun/tests/constantacceleration.cpp b/src/programs/mdrun/tests/constantacceleration.cpp new file mode 100644 index 0000000000..08e1b21c61 --- /dev/null +++ b/src/programs/mdrun/tests/constantacceleration.cpp @@ -0,0 +1,160 @@ +/* + * 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 + * \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; +class AccelerationGroupTest : + public MdrunTestFixture, + public ::testing::WithParamInterface +{ +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 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 referenceVelocities(c_groupSize, referenceVelocity); + + SCOPED_TRACE("Checking velocities of non-accelerated atoms"); + ArrayRef nonAcceleratedVelocities = frame.v().subArray(0, c_groupSize); + EXPECT_THAT(nonAcceleratedVelocities, Pointwise(RVecEq(tolerance), zeroVelocities)); + + SCOPED_TRACE("Checking velocities of accelerated atoms"); + ArrayRef 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 diff --git a/src/testutils/testmatchers.cpp b/src/testutils/testmatchers.cpp index 782d0bd162..1e7f5bd9fa 100644 --- a/src/testutils/testmatchers.cpp +++ b/src/testutils/testmatchers.cpp @@ -1,7 +1,7 @@ /* * 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. @@ -80,10 +80,14 @@ public: { 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.