Also adds a test for pulling with MTS.
(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.
.. 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
#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"
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.
*
{
std::string opt;
+ checkMtsConsistency(*ir, wi);
+
opt = "awh-nstout";
if (awhParams->nstOut <= 0)
{
{
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");
+ }
+ }
}
}
* \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
const t_mdatoms* mdatoms,
gmx::ArrayRef<const real> lambda,
const StepWorkload& stepWork,
- gmx::ForceWithVirial* forceWithVirial,
+ gmx::ForceWithVirial* forceWithVirialMtsLevel0,
+ gmx::ForceWithVirial* forceWithVirialMtsLevel1,
gmx_enerdata_t* enerd,
gmx_edsam* ed,
bool didNeighborSearch)
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);
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<double> 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<double> 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)
}
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");
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<int>(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)
{
#include <bitset>
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/enumerationhelpers.h"
struct t_inputrec;
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<MtsForceGroups, std::string> 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
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<const MtsLevel> 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<int>(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.
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"
"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;
::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
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff