#include "gromacs/mdlib/energyoutput.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/resethandler.h"
-#include "gromacs/mdlib/update.h"
#include "gromacs/mdrun/replicaexchange.h"
#include "gromacs/mdrun/shellfc.h"
#include "gromacs/mdrunutility/handlerestart.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/int64_to_int.h"
+#include "andersentemperaturecoupling.h"
#include "computeglobalselement.h"
#include "constraintelement.h"
+#include "expandedensembleelement.h"
+#include "firstorderpressurecoupling.h"
#include "forceelement.h"
+#include "mttk.h"
+#include "nosehooverchains.h"
#include "parrinellorahmanbarostat.h"
+#include "pullelement.h"
#include "simulatoralgorithm.h"
#include "statepropagatordata.h"
#include "velocityscalingtemperaturecoupling.h"
void ModularSimulator::addIntegrationElements(ModularSimulatorAlgorithmBuilder* builder)
{
+ const bool isTrotter = inputrecNvtTrotter(legacySimulatorData_->inputrec)
+ || inputrecNptTrotter(legacySimulatorData_->inputrec)
+ || inputrecNphTrotter(legacySimulatorData_->inputrec);
if (legacySimulatorData_->inputrec->eI == IntegrationAlgorithm::MD)
{
// The leap frog integration algorithm
builder->add<ForceElement>();
builder->add<StatePropagatorData::Element>();
if (legacySimulatorData_->inputrec->etc == TemperatureCoupling::VRescale
- || legacySimulatorData_->inputrec->etc == TemperatureCoupling::Berendsen)
+ || legacySimulatorData_->inputrec->etc == TemperatureCoupling::Berendsen
+ || legacySimulatorData_->inputrec->etc == TemperatureCoupling::NoseHoover)
{
- builder->add<VelocityScalingTemperatureCoupling>(-1,
+ builder->add<VelocityScalingTemperatureCoupling>(Offset(-1),
UseFullStepKE::No,
ReportPreviousStepConservedEnergy::No,
PropagatorTag("LeapFrogPropagator"));
}
- builder->add<Propagator<IntegrationStep::LeapFrog>>(
- PropagatorTag("LeapFrogPropagator"), legacySimulatorData_->inputrec->delta_t);
+ builder->add<Propagator<IntegrationStage::LeapFrog>>(
+ PropagatorTag("LeapFrogPropagator"), TimeStep(legacySimulatorData_->inputrec->delta_t));
if (legacySimulatorData_->constr)
{
builder->add<ConstraintsElement<ConstraintVariable::Positions>>();
}
+
+ if (legacySimulatorData_->inputrec->bPull)
+ {
+ builder->add<PullElement>();
+ }
+
builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>>();
if (legacySimulatorData_->inputrec->epc == PressureCoupling::ParrinelloRahman)
{
- builder->add<ParrinelloRahmanBarostat>(-1, PropagatorTag("LeapFrogPropagator"));
+ builder->add<ParrinelloRahmanBarostat>(Offset(-1), PropagatorTag("LeapFrogPropagator"));
+ }
+ else if (legacySimulatorData_->inputrec->epc == PressureCoupling::Berendsen
+ || legacySimulatorData_->inputrec->epc == PressureCoupling::CRescale)
+ {
+ builder->add<FirstOrderPressureCoupling>(0, ReportPreviousStepConservedEnergy::No);
}
}
- else if (legacySimulatorData_->inputrec->eI == IntegrationAlgorithm::VV)
+ else if (legacySimulatorData_->inputrec->eI == IntegrationAlgorithm::VV && !isTrotter)
{
// The velocity verlet integration algorithm
builder->add<ForceElement>();
- builder->add<Propagator<IntegrationStep::VelocitiesOnly>>(
- PropagatorTag("VelocityHalfStep"), 0.5 * legacySimulatorData_->inputrec->delta_t);
+ builder->add<Propagator<IntegrationStage::VelocitiesOnly>>(
+ PropagatorTag("VelocityHalfStep"), TimeStep(0.5 * legacySimulatorData_->inputrec->delta_t));
if (legacySimulatorData_->constr)
{
builder->add<ConstraintsElement<ConstraintVariable::Velocities>>();
}
builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>();
+ // Here, we have x / v / f at the full time step
builder->add<StatePropagatorData::Element>();
+ if (legacySimulatorData_->inputrec->bExpanded)
+ {
+ builder->add<ExpandedEnsembleElement>();
+ }
if (legacySimulatorData_->inputrec->etc == TemperatureCoupling::VRescale
|| legacySimulatorData_->inputrec->etc == TemperatureCoupling::Berendsen)
{
builder->add<VelocityScalingTemperatureCoupling>(
- 0,
+ Offset(0),
UseFullStepKE::Yes,
ReportPreviousStepConservedEnergy::Yes,
PropagatorTag("VelocityHalfAndPositionFullStep"));
}
- builder->add<Propagator<IntegrationStep::VelocityVerletPositionsAndVelocities>>(
- PropagatorTag("VelocityHalfAndPositionFullStep"), legacySimulatorData_->inputrec->delta_t);
+ else if (ETC_ANDERSEN(legacySimulatorData_->inputrec->etc))
+ {
+ builder->add<AndersenTemperatureCoupling>();
+ }
+ builder->add<Propagator<IntegrationStage::VelocityVerletPositionsAndVelocities>>(
+ PropagatorTag("VelocityHalfAndPositionFullStep"),
+ TimeStep(legacySimulatorData_->inputrec->delta_t));
if (legacySimulatorData_->constr)
{
builder->add<ConstraintsElement<ConstraintVariable::Positions>>();
}
+
+ if (legacySimulatorData_->inputrec->bPull)
+ {
+ builder->add<PullElement>();
+ }
+
builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>();
if (legacySimulatorData_->inputrec->epc == PressureCoupling::ParrinelloRahman)
{
- builder->add<ParrinelloRahmanBarostat>(-1, PropagatorTag("VelocityHalfStep"));
+ builder->add<ParrinelloRahmanBarostat>(Offset(-1), PropagatorTag("VelocityHalfStep"));
+ }
+ else if (legacySimulatorData_->inputrec->epc == PressureCoupling::Berendsen
+ || legacySimulatorData_->inputrec->epc == PressureCoupling::CRescale)
+ {
+ builder->add<FirstOrderPressureCoupling>(0, ReportPreviousStepConservedEnergy::Yes);
+ }
+ }
+ else if (legacySimulatorData_->inputrec->eI == IntegrationAlgorithm::VV && isTrotter)
+ {
+ // For a new simulation, avoid the first Trotter half step
+ const auto scheduleTrotterFirstHalfOnInitStep =
+ ((legacySimulatorData_->startingBehavior == StartingBehavior::NewSimulation)
+ ? ScheduleOnInitStep::No
+ : ScheduleOnInitStep::Yes);
+ // Define the tags and offsets for MTTK pressure scaling
+ const MttkPropagatorConnectionDetails mttkPropagatorConnectionDetails = {
+ PropagatorTag("ScaleMTTKXPre"), PropagatorTag("ScaleMTTKXPost"), Offset(0),
+ PropagatorTag("ScaleMTTKVPre1"), PropagatorTag("ScaleMTTKVPost1"), Offset(1),
+ PropagatorTag("ScaleMTTKVPre2"), PropagatorTag("ScaleMTTKVPost2"), Offset(0)
+ };
+
+ builder->add<ForceElement>();
+ // Propagate velocities from t-dt/2 to t
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::Mttk)
+ {
+ builder->add<Propagator<IntegrationStage::ScaleVelocities>>(
+ PropagatorTag("ScaleMTTKVPre1"));
+ }
+ builder->add<Propagator<IntegrationStage::VelocitiesOnly>>(
+ PropagatorTag("VelocityHalfStep1"),
+ TimeStep(0.5 * legacySimulatorData_->inputrec->delta_t));
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::Mttk)
+ {
+ builder->add<Propagator<IntegrationStage::ScaleVelocities>>(
+ PropagatorTag("ScaleMTTKVPost1"));
+ }
+ if (legacySimulatorData_->constr)
+ {
+ builder->add<ConstraintsElement<ConstraintVariable::Velocities>>();
+ }
+ builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>();
+
+ // Propagate extended system variables from t-dt/2 to t
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::Mttk)
+ {
+ builder->add<MttkElement>(
+ Offset(-1), scheduleTrotterFirstHalfOnInitStep, mttkPropagatorConnectionDetails);
+ }
+ if (legacySimulatorData_->inputrec->etc == TemperatureCoupling::NoseHoover)
+ {
+ builder->add<NoseHooverChainsElement>(NhcUsage::System,
+ Offset(-1),
+ UseFullStepKE::Yes,
+ scheduleTrotterFirstHalfOnInitStep,
+ PropagatorTag("ScaleNHC"));
+ builder->add<Propagator<IntegrationStage::ScaleVelocities>>(PropagatorTag("ScaleNHC"));
+ }
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::Mttk)
+ {
+ builder->add<NoseHooverChainsElement>(NhcUsage::Barostat,
+ Offset(-1),
+ UseFullStepKE::Yes,
+ scheduleTrotterFirstHalfOnInitStep,
+ mttkPropagatorConnectionDetails);
+ }
+ // We have a full state at time t here
+ builder->add<StatePropagatorData::Element>();
+ if (legacySimulatorData_->inputrec->bExpanded)
+ {
+ builder->add<ExpandedEnsembleElement>();
+ }
+
+ // Propagate extended system variables from t to t+dt/2
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::Mttk)
+ {
+ builder->add<NoseHooverChainsElement>(NhcUsage::Barostat,
+ Offset(0),
+ UseFullStepKE::Yes,
+ ScheduleOnInitStep::Yes,
+ mttkPropagatorConnectionDetails);
+ }
+ if (legacySimulatorData_->inputrec->etc == TemperatureCoupling::NoseHoover)
+ {
+ builder->add<NoseHooverChainsElement>(NhcUsage::System,
+ Offset(0),
+ UseFullStepKE::Yes,
+ ScheduleOnInitStep::Yes,
+ PropagatorTag("VelocityHalfStep2"));
+ }
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::Mttk)
+ {
+ builder->add<MttkElement>(Offset(0), ScheduleOnInitStep::Yes, mttkPropagatorConnectionDetails);
+ builder->add<Propagator<IntegrationStage::ScaleVelocities>>(
+ PropagatorTag("ScaleMTTKVPre2"));
+ }
+
+ // Propagate velocities from t to t+dt/2
+ builder->add<Propagator<IntegrationStage::VelocitiesOnly>>(
+ PropagatorTag("VelocityHalfStep2"),
+ TimeStep(0.5 * legacySimulatorData_->inputrec->delta_t));
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::Mttk)
+ {
+ builder->add<Propagator<IntegrationStage::ScaleVelocities>>(
+ PropagatorTag("ScaleMTTKVPost2"));
+ builder->add<Propagator<IntegrationStage::ScalePositions>>(
+ PropagatorTag("ScaleMTTKXPre"));
+ }
+ // Propagate positions from t to t+dt
+ builder->add<Propagator<IntegrationStage::PositionsOnly>>(
+ PropagatorTag("PositionFullStep"), TimeStep(legacySimulatorData_->inputrec->delta_t));
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::Mttk)
+ {
+ builder->add<Propagator<IntegrationStage::ScalePositions>>(
+ PropagatorTag("ScaleMTTKXPost"));
+ }
+ if (legacySimulatorData_->constr)
+ {
+ builder->add<ConstraintsElement<ConstraintVariable::Positions>>();
+ }
+
+ if (legacySimulatorData_->inputrec->bPull)
+ {
+ builder->add<PullElement>();
+ }
+
+ builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>();
+
+ // Propagate box from t to t+dt
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::Mttk)
+ {
+ builder->add<MttkBoxScaling>(mttkPropagatorConnectionDetails);
+ }
+ else if (legacySimulatorData_->inputrec->epc == PressureCoupling::CRescale)
+ {
+ // Legacy implementation allows combination of C-Rescale with Trotter Nose-Hoover
+ builder->add<FirstOrderPressureCoupling>(0, ReportPreviousStepConservedEnergy::Yes);
}
}
else
isInputCompatible =
isInputCompatible
&& conditionalAssert(
- !inputrec->useMts,
- "Multiple time stepping is not supported by the modular simulator.");
+ !inputrec->useMts,
+ "Multiple time stepping is not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
&& conditionalAssert(!doRerun, "Rerun is not supported by the modular simulator.");
- isInputCompatible = isInputCompatible
- && conditionalAssert(inputrec->etc == TemperatureCoupling::No
- || inputrec->etc == TemperatureCoupling::VRescale
- || inputrec->etc == TemperatureCoupling::Berendsen,
- "Only v-rescale and Berendsen thermostat are "
- "supported by the modular simulator.");
isInputCompatible =
isInputCompatible
- && conditionalAssert(
- inputrec->epc == PressureCoupling::No
- || inputrec->epc == PressureCoupling::ParrinelloRahman,
- "Only Parrinello-Rahman barostat is supported by the modular simulator.");
- isInputCompatible =
- isInputCompatible
- && conditionalAssert(
- !(inputrecNptTrotter(inputrec) || inputrecNphTrotter(inputrec)
- || inputrecNvtTrotter(inputrec)),
- "Legacy Trotter decomposition is not supported by the modular simulator.");
- isInputCompatible =
- isInputCompatible
- && conditionalAssert(inputrec->efep == FreeEnergyPerturbationType::No
- || inputrec->efep == FreeEnergyPerturbationType::Yes
- || inputrec->efep == FreeEnergyPerturbationType::SlowGrowth,
- "Expanded ensemble free energy calculation is not "
- "supported by the modular simulator.");
- isInputCompatible = isInputCompatible
- && conditionalAssert(!inputrec->bPull,
- "Pulling 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
isInputCompatible =
isInputCompatible
&& conditionalAssert(
- inputrec->deform[XX][XX] == 0.0 && inputrec->deform[XX][YY] == 0.0
- && inputrec->deform[XX][ZZ] == 0.0 && inputrec->deform[YY][XX] == 0.0
- && inputrec->deform[YY][YY] == 0.0 && inputrec->deform[YY][ZZ] == 0.0
- && inputrec->deform[ZZ][XX] == 0.0 && inputrec->deform[ZZ][YY] == 0.0
- && inputrec->deform[ZZ][ZZ] == 0.0,
- "Deformation is not supported by the modular simulator.");
+ inputrec->deform[XX][XX] == 0.0 && inputrec->deform[XX][YY] == 0.0
+ && inputrec->deform[XX][ZZ] == 0.0 && inputrec->deform[YY][XX] == 0.0
+ && inputrec->deform[YY][YY] == 0.0 && inputrec->deform[YY][ZZ] == 0.0
+ && inputrec->deform[ZZ][XX] == 0.0 && inputrec->deform[ZZ][YY] == 0.0
+ && inputrec->deform[ZZ][ZZ] == 0.0,
+ "Deformation is not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
&& conditionalAssert(gmx_mtop_interaction_count(globalTopology, IF_VSITE) == 0,
isInputCompatible =
isInputCompatible
&& conditionalAssert(
- gmx_mtop_ftype_count(globalTopology, F_ORIRES) == 0,
- "Orientation restraints are not supported by the modular simulator.");
+ gmx_mtop_ftype_count(globalTopology, F_ORIRES) == 0,
+ "Orientation restraints are not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
&& conditionalAssert(ms == nullptr,
}
else
{
- auto distantRestraintEnsembleEnvVar = getenv("GMX_DISRE_ENSEMBLE_SIZE");
+ auto* distantRestraintEnsembleEnvVar = getenv("GMX_DISRE_ENSEMBLE_SIZE");
numEnsembleRestraintSystems =
(ms != nullptr && distantRestraintEnsembleEnvVar != nullptr)
? static_cast<int>(strtol(distantRestraintEnsembleEnvVar, nullptr, 10))
isInputCompatible
&& conditionalAssert(!inputrec->bSimTemp,
"Simulated tempering is not supported by the modular simulator.");
- isInputCompatible = isInputCompatible
- && conditionalAssert(!inputrec->bExpanded,
- "Expanded ensemble simulations are not supported by "
- "the modular simulator.");
isInputCompatible =
isInputCompatible
&& conditionalAssert(!doEssentialDynamics,
isInputCompatible =
isInputCompatible
&& conditionalAssert(
- getenv("GMX_FORCE_UPDATE_DEFAULT_GPU") == nullptr,
- "Integration on the GPU is not supported by the modular simulator.");
+ getenv("GMX_FORCE_UPDATE_DEFAULT_GPU") == nullptr,
+ "Integration on the GPU is not supported by the modular simulator.");
// Modular simulator is centered around NS updates
// TODO: think how to handle nstlist == 0
isInputCompatible = isInputCompatible
isInputCompatible = isInputCompatible
&& conditionalAssert(!GMX_FAHCORE,
"GMX_FAHCORE not supported by the modular simulator.");
- GMX_RELEASE_ASSERT(
- isInputCompatible
- || !(inputrec->eI == IntegrationAlgorithm::VV
- && inputrec->epc == PressureCoupling::ParrinelloRahman),
- "Requested Parrinello-Rahman barostat with md-vv, but other options are not compatible "
- "with the modular simulator. The Parrinello-Rahman barostat is not implemented for "
- "md-vv in the legacy simulator. Use a different pressure control algorithm.");
-
+ if (!isInputCompatible
+ && (inputrec->eI == IntegrationAlgorithm::VV && inputrec->epc == PressureCoupling::ParrinelloRahman))
+ {
+ gmx_fatal(FARGS,
+ "Requested Parrinello-Rahman barostat with md-vv. This combination is only "
+ "available in the modular simulator. Some other selected options are, however, "
+ "only available in the legacy simulator. Use a different pressure control "
+ "algorithm.");
+ }
return isInputCompatible;
}