From 457e3d820e7e0e388d50d04f81a6267eb236e3aa Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Fri, 9 Oct 2020 15:19:10 +0200 Subject: [PATCH] Add MTS support for pull and AWH Also adds a test for pulling with MTS. --- docs/user-guide/mdp-options.rst | 10 ++-- .../applied_forces/awh/read_params.cpp | 55 +++++++++++++++++++ src/gromacs/gmxpreprocess/readir.cpp | 13 +++++ src/gromacs/mdlib/sim_util.cpp | 50 +++++++++++------ src/gromacs/mdtypes/multipletimestepping.cpp | 7 +-- src/gromacs/mdtypes/multipletimestepping.h | 21 ++++++- .../mdrun/tests/multiple_time_stepping.cpp | 25 ++++++++- 7 files changed, 150 insertions(+), 31 deletions(-) diff --git a/docs/user-guide/mdp-options.rst b/docs/user-guide/mdp-options.rst index d6f9d3c63f..251172ae96 100644 --- a/docs/user-guide/mdp-options.rst +++ b/docs/user-guide/mdp-options.rst @@ -254,9 +254,9 @@ Run control (longrange-nonbonded nonbonded pair dihedral) A list of force groups that will be evaluated only every :mdp:`mts-level2-factor` steps. Supported entries are: - ``longrange-nonbonded``, ``nonbonded``, ``pair``, ``dihedral`` and - ``angle``. With ``pair`` the listed pair forces (such as 1-4) are - selected. With ``dihedral`` all dihedrals are selected, including cmap. + ``longrange-nonbonded``, ``nonbonded``, ``pair``, ``dihedral``, ``angle``, + ``pull`` and ``awh``. With ``pair`` the listed pair forces (such as 1-4) + are selected. With ``dihedral`` all dihedrals are selected, including cmap. All other forces, including all restraints, are evaluated and integrated every step. When PME or Ewald is used for electrostatics and/or LJ interactions, ``longrange-nonbonded`` has to be entered here. @@ -2092,13 +2092,15 @@ AWH adaptive biasing .. mdp-value:: pull The pull module is providing the reaction coordinate for this dimension. + With multiple time-stepping, AWH and pull should be in the same MTS level. .. mdp-value:: fep-lambda The free energy lambda state is the reaction coordinate for this dimension. The lambda states to use are specified by :mdp:`fep-lambdas`, :mdp:`vdw-lambdas`, :mdp:`coul-lambdas` etc. This is not compatible with delta-lambda. It also requires - calc-lambda-neighbors to be -1. + calc-lambda-neighbors to be -1. With multiple time-stepping, AWH should + be in the slow level. .. mdp:: awh1-dim1-coord-index diff --git a/src/gromacs/applied_forces/awh/read_params.cpp b/src/gromacs/applied_forces/awh/read_params.cpp index d301bd83d9..e5945f17d6 100644 --- a/src/gromacs/applied_forces/awh/read_params.cpp +++ b/src/gromacs/applied_forces/awh/read_params.cpp @@ -48,6 +48,7 @@ #include "gromacs/mdtypes/awh_params.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/md_enums.h" +#include "gromacs/mdtypes/multipletimestepping.h" #include "gromacs/mdtypes/pull_params.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pulling/pull.h" @@ -75,6 +76,58 @@ const char* eawhcoordprovider_names[eawhcoordproviderNR + 1] = { "pull", "fep-la namespace { + +/*! \brief + * Check multiple time-stepping consistency between AWH and pull and/or FEP + * + * \param[in] inputrec Inputput parameter struct. + * \param[in,out] wi Struct for bookeeping warnings. + */ +void checkMtsConsistency(const t_inputrec& inputrec, warninp_t wi) +{ + if (!inputrec.useMts) + { + return; + } + + GMX_RELEASE_ASSERT(inputrec.mtsLevels.size() == 2, "Only 2 MTS levels supported here"); + + bool usesPull = false; + bool usesFep = false; + for (int b = 0; b < inputrec.awhParams->numBias; b++) + { + const auto& biasParams = inputrec.awhParams->awhBiasParams[b]; + for (int d = 0; d < biasParams.ndim; d++) + { + switch (biasParams.dimParams[d].eCoordProvider) + { + case eawhcoordproviderPULL: usesPull = true; break; + case eawhcoordproviderFREE_ENERGY_LAMBDA: usesFep = true; break; + default: GMX_RELEASE_ASSERT(false, "Unsupported coord provider"); + } + } + } + const int awhMtsLevel = forceGroupMtsLevel(inputrec.mtsLevels, MtsForceGroups::Awh); + if (usesPull && forceGroupMtsLevel(inputrec.mtsLevels, MtsForceGroups::Pull) != awhMtsLevel) + { + warning_error(wi, + "When AWH is applied to pull coordinates, pull and AWH should be computed at " + "the same MTS level"); + } + if (usesFep && awhMtsLevel != ssize(inputrec.mtsLevels) - 1) + { + warning_error(wi, + "When AWH is applied to the free-energy lambda with MTS, AWH should be " + "computed at the slow MTS level"); + } + + if (inputrec.awhParams->nstSampleCoord % inputrec.mtsLevels[awhMtsLevel].stepFactor != 0) + { + warning_error(wi, + "With MTS applied to AWH, awh-nstsample should be a multiple of mts-factor"); + } +} + /*! \brief * Check the parameters of an AWH bias pull dimension. * @@ -751,6 +804,8 @@ void checkAwhParams(const AwhParams* awhParams, const t_inputrec* ir, warninp_t { std::string opt; + checkMtsConsistency(*ir, wi); + opt = "awh-nstout"; if (awhParams->nstOut <= 0) { diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp index 270617bbab..b2d504ca23 100644 --- a/src/gromacs/gmxpreprocess/readir.cpp +++ b/src/gromacs/gmxpreprocess/readir.cpp @@ -311,6 +311,19 @@ static void setupMtsLevels(gmx::ArrayRef mtsLevels, { checkMtsRequirement(ir, "nstdhdl", ir.fepvals->nstdhdl, wi); } + + if (ir.bPull) + { + const int pullMtsLevel = gmx::forceGroupMtsLevel(ir.mtsLevels, gmx::MtsForceGroups::Pull); + if (ir.pull->nstxout % ir.mtsLevels[pullMtsLevel].stepFactor != 0) + { + warning_error(wi, "pull-nstxout should be a multiple of mts-factor"); + } + if (ir.pull->nstfout % ir.mtsLevels[pullMtsLevel].stepFactor != 0) + { + warning_error(wi, "pull-nstfout should be a multiple of mts-factor"); + } + } } } diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index ef27166978..17a9e44bfa 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -580,7 +580,8 @@ static bool haveSpecialForces(const t_inputrec& inputrec, * \param[in] mdatoms Per atom properties * \param[in] lambda Array of free-energy lambda values * \param[in] stepWork Step schedule flags - * \param[in,out] forceWithVirial Force and virial buffers + * \param[in,out] forceWithVirialMtsLevel0 Force and virial for MTS level0 forces + * \param[in,out] forceWithVirialMtsLevel1 Force and virial for MTS level1 forces, can be nullptr * \param[in,out] enerd Energy buffer * \param[in,out] ed Essential dynamics pointer * \param[in] didNeighborSearch Tells if we did neighbor searching this step, used for ED sampling @@ -604,7 +605,8 @@ static void computeSpecialForces(FILE* fplog, const t_mdatoms* mdatoms, gmx::ArrayRef lambda, const StepWorkload& stepWork, - gmx::ForceWithVirial* forceWithVirial, + gmx::ForceWithVirial* forceWithVirialMtsLevel0, + gmx::ForceWithVirial* forceWithVirialMtsLevel1, gmx_enerdata_t* enerd, gmx_edsam* ed, bool didNeighborSearch) @@ -615,7 +617,7 @@ static void computeSpecialForces(FILE* fplog, if (stepWork.computeForces) { gmx::ForceProviderInput forceProviderInput(x, *mdatoms, t, box, *cr); - gmx::ForceProviderOutput forceProviderOutput(forceWithVirial, enerd); + gmx::ForceProviderOutput forceProviderOutput(forceWithVirialMtsLevel0, enerd); /* Collect forces from modules */ forceProviders->calculateForces(forceProviderInput, &forceProviderOutput); @@ -623,26 +625,36 @@ static void computeSpecialForces(FILE* fplog, if (inputrec->bPull && pull_have_potential(pull_work)) { - pull_potential_wrapper(cr, inputrec, box, x, forceWithVirial, mdatoms, enerd, pull_work, - lambda.data(), t, wcycle); + const int mtsLevel = forceGroupMtsLevel(inputrec->mtsLevels, gmx::MtsForceGroups::Pull); + if (mtsLevel == 0 || stepWork.computeSlowForces) + { + auto& forceWithVirial = (mtsLevel == 0) ? forceWithVirialMtsLevel0 : forceWithVirialMtsLevel1; + pull_potential_wrapper(cr, inputrec, box, x, forceWithVirial, mdatoms, enerd, pull_work, + lambda.data(), t, wcycle); + } } if (awh) { - const bool needForeignEnergyDifferences = awh->needForeignEnergyDifferences(step); - std::vector foreignLambdaDeltaH, foreignLambdaDhDl; - if (needForeignEnergyDifferences) + const int mtsLevel = forceGroupMtsLevel(inputrec->mtsLevels, gmx::MtsForceGroups::Pull); + if (mtsLevel == 0 || stepWork.computeSlowForces) { - enerd->foreignLambdaTerms.finalizePotentialContributions(enerd->dvdl_lin, lambda, - *inputrec->fepvals); - std::tie(foreignLambdaDeltaH, foreignLambdaDhDl) = enerd->foreignLambdaTerms.getTerms(cr); - } + const bool needForeignEnergyDifferences = awh->needForeignEnergyDifferences(step); + std::vector foreignLambdaDeltaH, foreignLambdaDhDl; + if (needForeignEnergyDifferences) + { + enerd->foreignLambdaTerms.finalizePotentialContributions(enerd->dvdl_lin, lambda, + *inputrec->fepvals); + std::tie(foreignLambdaDeltaH, foreignLambdaDhDl) = enerd->foreignLambdaTerms.getTerms(cr); + } - enerd->term[F_COM_PULL] += awh->applyBiasForcesAndUpdateBias( - inputrec->pbcType, mdatoms->massT, foreignLambdaDeltaH, foreignLambdaDhDl, box, - forceWithVirial, t, step, wcycle, fplog); + auto& forceWithVirial = (mtsLevel == 0) ? forceWithVirialMtsLevel0 : forceWithVirialMtsLevel1; + enerd->term[F_COM_PULL] += awh->applyBiasForcesAndUpdateBias( + inputrec->pbcType, mdatoms->massT, foreignLambdaDeltaH, foreignLambdaDhDl, box, + forceWithVirial, t, step, wcycle, fplog); + } } - rvec* f = as_rvec_array(forceWithVirial->force_.data()); + rvec* f = as_rvec_array(forceWithVirialMtsLevel0->force_.data()); /* Add the forces from enforced rotation potentials (if any) */ if (inputrec->bRot) @@ -1781,8 +1793,10 @@ void do_force(FILE* fplog, } computeSpecialForces(fplog, cr, inputrec, awh, enforcedRotation, imdSession, pull_work, step, t, - wcycle, fr->forceProviders, box, x.unpaddedArrayRef(), mdatoms, lambda, stepWork, - &forceOutMtsLevel0.forceWithVirial(), enerd, ed, stepWork.doNeighborSearch); + wcycle, fr->forceProviders, box, x.unpaddedArrayRef(), mdatoms, lambda, + stepWork, &forceOutMtsLevel0.forceWithVirial(), + forceOutMtsLevel1 ? &forceOutMtsLevel1->forceWithVirial() : nullptr, enerd, + ed, stepWork.doNeighborSearch); GMX_ASSERT(!(nonbondedAtMtsLevel1 && stepWork.useGpuFBufferOps), "The schedule below does not allow for nonbonded MTS with GPU buffer ops"); diff --git a/src/gromacs/mdtypes/multipletimestepping.cpp b/src/gromacs/mdtypes/multipletimestepping.cpp index be56600040..02bcf41be3 100644 --- a/src/gromacs/mdtypes/multipletimestepping.cpp +++ b/src/gromacs/mdtypes/multipletimestepping.cpp @@ -67,10 +67,9 @@ void assertMtsRequirements(const t_inputrec& ir) GMX_RELEASE_ASSERT(ir.mtsLevels[0].stepFactor == 1, "Base MTS step should be 1"); - GMX_RELEASE_ASSERT( - (!EEL_FULL(ir.coulombtype) && !EVDW_PME(ir.vdwtype)) - || ir.mtsLevels.back().forceGroups[static_cast(MtsForceGroups::LongrangeNonbonded)], - "Long-range nonbondeds should be in the highest MTS level"); + GMX_RELEASE_ASSERT((!EEL_FULL(ir.coulombtype) && !EVDW_PME(ir.vdwtype)) + || forceGroupMtsLevel(ir.mtsLevels, MtsForceGroups::LongrangeNonbonded) > 0, + "Long-range nonbondeds should be in the highest MTS level"); for (const auto& mtsLevel : ir.mtsLevels) { diff --git a/src/gromacs/mdtypes/multipletimestepping.h b/src/gromacs/mdtypes/multipletimestepping.h index a1a953a3a7..5796114005 100644 --- a/src/gromacs/mdtypes/multipletimestepping.h +++ b/src/gromacs/mdtypes/multipletimestepping.h @@ -37,6 +37,7 @@ #include +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/enumerationhelpers.h" struct t_inputrec; @@ -51,12 +52,14 @@ enum class MtsForceGroups : int Nonbonded, //!< Non-bonded pair interactions Pair, //!< Bonded pair interactions Dihedral, //!< Dihedrals, including cmap (not restraints) - Angle, //! Bonded angle potentials (not restraints) - Count //! The number of groups above + Angle, //!< Bonded angle potentials (not restraints) + Pull, //!< COM pulling + Awh, //!< Accelerated weight histogram method + Count //!< The number of groups above }; static const gmx::EnumerationArray mtsForceGroupNames = { - "longrange-nonbonded", "nonbonded", "pair", "dihedral", "angle" + "longrange-nonbonded", "nonbonded", "pair", "dihedral", "angle", "pull", "awh" }; //! Setting for a single level for multiple time step integration @@ -68,6 +71,18 @@ struct MtsLevel int stepFactor; }; +/*! \brief Returns the MTS level at which a force group is to be computed + * + * \param[in] mtsLevels List of force groups for each MTS level, can be empty without MTS + * \param[in] mtsForceGroup The force group to query the MTS level for + */ +static inline int forceGroupMtsLevel(ArrayRef mtsLevels, const MtsForceGroups mtsForceGroup) +{ + GMX_ASSERT(mtsLevels.empty() || mtsLevels.size() == 2, "Only 0 or 2 MTS levels are supported"); + + return (mtsLevels.empty() || mtsLevels[0].forceGroups[static_cast(mtsForceGroup)]) ? 0 : 1; +}; + /*! \brief Returns the interval in steps at which the non-bonded pair forces are calculated * * Note: returns 1 when multiple time-stepping is not activated. diff --git a/src/programs/mdrun/tests/multiple_time_stepping.cpp b/src/programs/mdrun/tests/multiple_time_stepping.cpp index 44d6ffd990..a3022f3d43 100644 --- a/src/programs/mdrun/tests/multiple_time_stepping.cpp +++ b/src/programs/mdrun/tests/multiple_time_stepping.cpp @@ -90,6 +90,8 @@ TEST_P(MtsComparisonTest, WithinTolerances) SCOPED_TRACE(formatString("Comparing for '%s' no MTS with MTS scheme '%s'", simulationName.c_str(), mtsScheme.c_str())); + const bool isPullTest = (mtsScheme.find("pull") != std::string::npos); + const int numSteps = 4; auto sharedMdpOptions = gmx::formatString( "integrator = md\n" @@ -97,12 +99,27 @@ TEST_P(MtsComparisonTest, WithinTolerances) "nsteps = %d\n" "verlet-buffer-tolerance = -1\n" "rlist = 1.0\n" - "coulomb-type = PME\n" + "coulomb-type = %s\n" "vdw-type = cut-off\n" "rcoulomb = 0.9\n" "rvdw = 0.9\n" "constraints = h-bonds\n", - numSteps); + numSteps, isPullTest ? "reaction-field" : "PME"); + + if (isPullTest) + { + sharedMdpOptions += + "pull = yes\n" + "pull-ngroups = 2\n" + "pull-group1-name = FirstWaterMolecule\n" + "pull-group2-name = SecondWaterMolecule\n" + "pull-ncoords = 1\n" + "pull-coord1-type = umbrella\n" + "pull-coord1-geometry = distance\n" + "pull-coord1-groups = 1 2\n" + "pull-coord1-init = 1\n" + "pull-coord1-k = 10000\n"; + } // set nstfout to > numSteps so we only write forces at step 0 const int nstfout = 2 * numSteps; @@ -189,6 +206,10 @@ INSTANTIATE_TEST_CASE_P( ::testing::Values("longrange-nonbonded", "longrange-nonbonded nonbonded pair dihedral"))); +INSTANTIATE_TEST_CASE_P(MultipleTimeSteppingIsNearSingleTimeSteppingPull, + MtsComparisonTest, + ::testing::Combine(::testing::Values("spc2"), ::testing::Values("pull"))); + } // namespace } // namespace test } // namespace gmx -- 2.22.0