/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 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.
#include "gromacs/ewald/pme.h"
#include "gromacs/ewald/pme_load_balancing.h"
#include "gromacs/ewald/pme_pp.h"
-#include "gromacs/gmxlib/network.h"
+#include "gromacs/fileio/checkpoint.h"
#include "gromacs/gmxlib/nrnb.h"
#include "gromacs/listed_forces/listed_forces.h"
#include "gromacs/mdlib/checkpointhandler.h"
#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/nbnxm/nbnxm.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/topology/topology.h"
+#include "gromacs/trajectory/trajectoryframe.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/int64_to_int.h"
-#include "compositesimulatorelement.h"
+#include "andersentemperaturecoupling.h"
#include "computeglobalselement.h"
#include "constraintelement.h"
-#include "energydata.h"
+#include "expandedensembleelement.h"
+#include "firstorderpressurecoupling.h"
#include "forceelement.h"
-#include "freeenergyperturbationdata.h"
+#include "mttk.h"
+#include "nosehooverchains.h"
#include "parrinellorahmanbarostat.h"
-#include "propagator.h"
-#include "signallers.h"
+#include "pullelement.h"
#include "simulatoralgorithm.h"
#include "statepropagatordata.h"
-#include "trajectoryelement.h"
-#include "vrescalethermostat.h"
+#include "velocityscalingtemperaturecoupling.h"
namespace gmx
{
.asParagraph()
.appendText("Using the modular simulator.");
- ModularSimulatorAlgorithmBuilder algorithmBuilder(compat::make_not_null(legacySimulatorData_.get()));
- auto algorithm = algorithmBuilder.build();
+ ModularSimulatorAlgorithmBuilder algorithmBuilder(compat::make_not_null(legacySimulatorData_),
+ std::move(checkpointDataHolder_));
+ addIntegrationElements(&algorithmBuilder);
+ auto algorithm = algorithmBuilder.build();
while (const auto* task = algorithm.getNextTask())
{
}
}
-std::unique_ptr<ISimulatorElement> ModularSimulatorAlgorithmBuilder::buildForces(
- SignallerBuilder<NeighborSearchSignaller>* neighborSearchSignallerBuilder,
- SignallerBuilder<EnergySignaller>* energySignallerBuilder,
- StatePropagatorData* statePropagatorDataPtr,
- EnergyData* energyDataPtr,
- FreeEnergyPerturbationData* freeEnergyPerturbationDataPtr,
- TopologyHolder::Builder* topologyHolderBuilder)
+void ModularSimulator::addIntegrationElements(ModularSimulatorAlgorithmBuilder* builder)
{
- const bool isVerbose = legacySimulatorData_->mdrunOptions.verbose;
- const bool isDynamicBox = inputrecDynamicBox(legacySimulatorData_->inputrec);
-
- auto forceElement = std::make_unique<ForceElement>(
- statePropagatorDataPtr, energyDataPtr, freeEnergyPerturbationDataPtr, isVerbose, isDynamicBox,
- legacySimulatorData_->fplog, legacySimulatorData_->cr, legacySimulatorData_->inputrec,
- legacySimulatorData_->mdAtoms, legacySimulatorData_->nrnb, legacySimulatorData_->fr,
- legacySimulatorData_->wcycle, legacySimulatorData_->runScheduleWork,
- legacySimulatorData_->vsite, legacySimulatorData_->imdSession,
- legacySimulatorData_->pull_work, legacySimulatorData_->constr,
- legacySimulatorData_->top_global, legacySimulatorData_->enforcedRotation);
- topologyHolderBuilder->registerClient(forceElement.get());
- neighborSearchSignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
- energySignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
-
- // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
- return std::move(forceElement);
-}
-
-std::unique_ptr<ISimulatorElement> ModularSimulatorAlgorithmBuilder::buildIntegrator(
- SignallerBuilder<NeighborSearchSignaller>* neighborSearchSignallerBuilder,
- SignallerBuilder<LastStepSignaller>* lastStepSignallerBuilder,
- SignallerBuilder<EnergySignaller>* energySignallerBuilder,
- SignallerBuilder<LoggingSignaller>* loggingSignallerBuilder,
- SignallerBuilder<TrajectorySignaller>* trajectorySignallerBuilder,
- TrajectoryElementBuilder* trajectoryElementBuilder,
- std::vector<ICheckpointHelperClient*>* checkpointClients,
- compat::not_null<StatePropagatorData*> statePropagatorDataPtr,
- compat::not_null<EnergyData*> energyDataPtr,
- FreeEnergyPerturbationData* freeEnergyPerturbationDataPtr,
- bool hasReadEkinState,
- TopologyHolder::Builder* topologyHolderBuilder,
- GlobalCommunicationHelper* globalCommunicationHelper)
-{
- auto forceElement = buildForces(neighborSearchSignallerBuilder, energySignallerBuilder,
- statePropagatorDataPtr, energyDataPtr,
- freeEnergyPerturbationDataPtr, topologyHolderBuilder);
-
- // list of elements owned by the simulator composite object
- std::vector<std::unique_ptr<ISimulatorElement>> elementsOwnershipList;
- // call list of the simulator composite object
- std::vector<compat::not_null<ISimulatorElement*>> elementCallList;
-
- std::function<void()> needToCheckNumberOfBondedInteractions;
- if (legacySimulatorData_->inputrec->eI == eiMD)
+ const bool isTrotter = inputrecNvtTrotter(legacySimulatorData_->inputrec)
+ || inputrecNptTrotter(legacySimulatorData_->inputrec)
+ || inputrecNphTrotter(legacySimulatorData_->inputrec);
+ if (legacySimulatorData_->inputrec->eI == IntegrationAlgorithm::MD)
{
- auto computeGlobalsElement =
- std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>>(
- statePropagatorDataPtr, energyDataPtr, freeEnergyPerturbationDataPtr,
- globalCommunicationHelper->simulationSignals(),
- globalCommunicationHelper->nstglobalcomm(), legacySimulatorData_->fplog,
- legacySimulatorData_->mdlog, legacySimulatorData_->cr,
- legacySimulatorData_->inputrec, legacySimulatorData_->mdAtoms,
- legacySimulatorData_->nrnb, legacySimulatorData_->wcycle,
- legacySimulatorData_->fr, legacySimulatorData_->top_global,
- legacySimulatorData_->constr, hasReadEkinState);
- topologyHolderBuilder->registerClient(computeGlobalsElement.get());
- energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
- trajectorySignallerBuilder->registerSignallerClient(
- compat::make_not_null(computeGlobalsElement.get()));
-
- globalCommunicationHelper->setCheckBondedInteractionsCallback(
- computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback());
-
- auto propagator = std::make_unique<Propagator<IntegrationStep::LeapFrog>>(
- legacySimulatorData_->inputrec->delta_t, statePropagatorDataPtr,
- legacySimulatorData_->mdAtoms, legacySimulatorData_->wcycle);
+ // 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::NoseHoover)
+ {
+ builder->add<VelocityScalingTemperatureCoupling>(Offset(-1),
+ UseFullStepKE::No,
+ ReportPreviousStepConservedEnergy::No,
+ PropagatorTag("LeapFrogPropagator"));
+ }
+ builder->add<Propagator<IntegrationStage::LeapFrog>>(
+ PropagatorTag("LeapFrogPropagator"), TimeStep(legacySimulatorData_->inputrec->delta_t));
+ if (legacySimulatorData_->constr)
+ {
+ builder->add<ConstraintsElement<ConstraintVariable::Positions>>();
+ }
- addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
- auto stateElement = compat::make_not_null(statePropagatorDataPtr->element());
- trajectoryElementBuilder->registerWriterClient(stateElement);
- trajectorySignallerBuilder->registerSignallerClient(stateElement);
- lastStepSignallerBuilder->registerSignallerClient(stateElement);
- checkpointClients->emplace_back(stateElement);
- // we have a full microstate at time t here!
- addToCallList(stateElement, elementCallList);
- if (legacySimulatorData_->inputrec->etc == etcVRESCALE)
+ if (legacySimulatorData_->inputrec->bPull)
{
- // TODO: With increased complexity of the propagator, this will need further development,
- // e.g. using propagators templated for velocity propagation policies and a builder
- propagator->setNumVelocityScalingVariables(legacySimulatorData_->inputrec->opts.ngtc);
- auto thermostat = std::make_unique<VRescaleThermostat>(
- legacySimulatorData_->inputrec->nsttcouple, -1, false,
- legacySimulatorData_->inputrec->ld_seed, legacySimulatorData_->inputrec->opts.ngtc,
- legacySimulatorData_->inputrec->delta_t * legacySimulatorData_->inputrec->nsttcouple,
- legacySimulatorData_->inputrec->opts.ref_t,
- legacySimulatorData_->inputrec->opts.tau_t, legacySimulatorData_->inputrec->opts.nrdf,
- energyDataPtr, propagator->viewOnVelocityScaling(),
- propagator->velocityScalingCallback(), legacySimulatorData_->state_global,
- legacySimulatorData_->cr, legacySimulatorData_->inputrec->bContinuation);
- checkpointClients->emplace_back(thermostat.get());
- energyDataPtr->setVRescaleThermostat(thermostat.get());
- addToCallListAndMove(std::move(thermostat), elementCallList, elementsOwnershipList);
+ builder->add<PullElement>();
}
- std::unique_ptr<ParrinelloRahmanBarostat> prBarostat = nullptr;
- if (legacySimulatorData_->inputrec->epc == epcPARRINELLORAHMAN)
+ builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>>();
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::ParrinelloRahman)
+ {
+ builder->add<ParrinelloRahmanBarostat>(Offset(-1), PropagatorTag("LeapFrogPropagator"));
+ }
+ else if (legacySimulatorData_->inputrec->epc == PressureCoupling::Berendsen
+ || legacySimulatorData_->inputrec->epc == PressureCoupling::CRescale)
{
- // Building the PR barostat here since it needs access to the propagator
- // and we want to be able to move the propagator object
- prBarostat = std::make_unique<ParrinelloRahmanBarostat>(
- legacySimulatorData_->inputrec->nstpcouple, -1,
- legacySimulatorData_->inputrec->delta_t * legacySimulatorData_->inputrec->nstpcouple,
- legacySimulatorData_->inputrec->init_step, propagator->viewOnPRScalingMatrix(),
- propagator->prScalingCallback(), statePropagatorDataPtr, energyDataPtr,
- legacySimulatorData_->fplog, legacySimulatorData_->inputrec,
- legacySimulatorData_->mdAtoms, legacySimulatorData_->state_global,
- legacySimulatorData_->cr, legacySimulatorData_->inputrec->bContinuation);
- energyDataPtr->setParrinelloRahamnBarostat(prBarostat.get());
- checkpointClients->emplace_back(prBarostat.get());
+ builder->add<FirstOrderPressureCoupling>(0, ReportPreviousStepConservedEnergy::No);
}
- addToCallListAndMove(std::move(propagator), elementCallList, elementsOwnershipList);
+ }
+ else if (legacySimulatorData_->inputrec->eI == IntegrationAlgorithm::VV && !isTrotter)
+ {
+ // The velocity verlet integration algorithm
+ builder->add<ForceElement>();
+ builder->add<Propagator<IntegrationStage::VelocitiesOnly>>(
+ PropagatorTag("VelocityHalfStep"), TimeStep(0.5 * legacySimulatorData_->inputrec->delta_t));
if (legacySimulatorData_->constr)
{
- auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Positions>>(
- legacySimulatorData_->constr, statePropagatorDataPtr, energyDataPtr,
- freeEnergyPerturbationDataPtr, MASTER(legacySimulatorData_->cr),
- legacySimulatorData_->fplog, legacySimulatorData_->inputrec,
- legacySimulatorData_->mdAtoms->mdatoms());
- auto constraintElementPtr = compat::make_not_null(constraintElement.get());
- energySignallerBuilder->registerSignallerClient(constraintElementPtr);
- trajectorySignallerBuilder->registerSignallerClient(constraintElementPtr);
- loggingSignallerBuilder->registerSignallerClient(constraintElementPtr);
+ 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>(
+ Offset(0),
+ UseFullStepKE::Yes,
+ ReportPreviousStepConservedEnergy::Yes,
+ PropagatorTag("VelocityHalfAndPositionFullStep"));
+ }
+ 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>>();
+ }
- addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
+ if (legacySimulatorData_->inputrec->bPull)
+ {
+ builder->add<PullElement>();
}
- addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
- auto energyElement = compat::make_not_null(energyDataPtr->element());
- trajectoryElementBuilder->registerWriterClient(energyElement);
- trajectorySignallerBuilder->registerSignallerClient(energyElement);
- energySignallerBuilder->registerSignallerClient(energyElement);
- checkpointClients->emplace_back(energyElement);
- // we have the energies at time t here!
- addToCallList(energyElement, elementCallList);
- if (prBarostat)
+ builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>();
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::ParrinelloRahman)
{
- addToCallListAndMove(std::move(prBarostat), elementCallList, elementsOwnershipList);
+ 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 == eiVV)
+ else if (legacySimulatorData_->inputrec->eI == IntegrationAlgorithm::VV && isTrotter)
{
- auto computeGlobalsElement =
- std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>(
- statePropagatorDataPtr, energyDataPtr, freeEnergyPerturbationDataPtr,
- globalCommunicationHelper->simulationSignals(),
- globalCommunicationHelper->nstglobalcomm(), legacySimulatorData_->fplog,
- legacySimulatorData_->mdlog, legacySimulatorData_->cr,
- legacySimulatorData_->inputrec, legacySimulatorData_->mdAtoms,
- legacySimulatorData_->nrnb, legacySimulatorData_->wcycle,
- legacySimulatorData_->fr, legacySimulatorData_->top_global,
- legacySimulatorData_->constr, hasReadEkinState);
- topologyHolderBuilder->registerClient(computeGlobalsElement.get());
- energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
- trajectorySignallerBuilder->registerSignallerClient(
- compat::make_not_null(computeGlobalsElement.get()));
-
- globalCommunicationHelper->setCheckBondedInteractionsCallback(
- computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback());
-
- auto propagatorVelocities = std::make_unique<Propagator<IntegrationStep::VelocitiesOnly>>(
- legacySimulatorData_->inputrec->delta_t * 0.5, statePropagatorDataPtr,
- legacySimulatorData_->mdAtoms, legacySimulatorData_->wcycle);
- auto propagatorVelocitiesAndPositions =
- std::make_unique<Propagator<IntegrationStep::VelocityVerletPositionsAndVelocities>>(
- legacySimulatorData_->inputrec->delta_t, statePropagatorDataPtr,
- legacySimulatorData_->mdAtoms, legacySimulatorData_->wcycle);
+ // 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>>();
- addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
+ // 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>();
+ }
- std::unique_ptr<ParrinelloRahmanBarostat> prBarostat = nullptr;
- if (legacySimulatorData_->inputrec->epc == epcPARRINELLORAHMAN)
+ // Propagate extended system variables from t to t+dt/2
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::Mttk)
{
- // Building the PR barostat here since it needs access to the propagator
- // and we want to be able to move the propagator object
- prBarostat = std::make_unique<ParrinelloRahmanBarostat>(
- legacySimulatorData_->inputrec->nstpcouple, -1,
- legacySimulatorData_->inputrec->delta_t * legacySimulatorData_->inputrec->nstpcouple,
- legacySimulatorData_->inputrec->init_step,
- propagatorVelocities->viewOnPRScalingMatrix(),
- propagatorVelocities->prScalingCallback(), statePropagatorDataPtr,
- energyDataPtr, legacySimulatorData_->fplog, legacySimulatorData_->inputrec,
- legacySimulatorData_->mdAtoms, legacySimulatorData_->state_global,
- legacySimulatorData_->cr, legacySimulatorData_->inputrec->bContinuation);
- energyDataPtr->setParrinelloRahamnBarostat(prBarostat.get());
- checkpointClients->emplace_back(prBarostat.get());
+ builder->add<NoseHooverChainsElement>(NhcUsage::Barostat,
+ Offset(0),
+ UseFullStepKE::Yes,
+ ScheduleOnInitStep::Yes,
+ mttkPropagatorConnectionDetails);
}
- addToCallListAndMove(std::move(propagatorVelocities), elementCallList, elementsOwnershipList);
- if (legacySimulatorData_->constr)
+ 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)
{
- auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Velocities>>(
- legacySimulatorData_->constr, statePropagatorDataPtr, energyDataPtr,
- freeEnergyPerturbationDataPtr, MASTER(legacySimulatorData_->cr),
- legacySimulatorData_->fplog, legacySimulatorData_->inputrec,
- legacySimulatorData_->mdAtoms->mdatoms());
- energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
- trajectorySignallerBuilder->registerSignallerClient(
- compat::make_not_null(constraintElement.get()));
- loggingSignallerBuilder->registerSignallerClient(
- compat::make_not_null(constraintElement.get()));
+ builder->add<MttkElement>(Offset(0), ScheduleOnInitStep::Yes, mttkPropagatorConnectionDetails);
+ builder->add<Propagator<IntegrationStage::ScaleVelocities>>(
+ PropagatorTag("ScaleMTTKVPre2"));
+ }
- addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
+ // 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"));
}
- addToCallList(compat::make_not_null(computeGlobalsElement.get()), elementCallList);
- auto stateElement = compat::make_not_null(statePropagatorDataPtr->element());
- trajectoryElementBuilder->registerWriterClient(stateElement);
- trajectorySignallerBuilder->registerSignallerClient(stateElement);
- lastStepSignallerBuilder->registerSignallerClient(stateElement);
- checkpointClients->emplace_back(stateElement);
- // we have a full microstate at time t here!
- addToCallList(stateElement, elementCallList);
- if (legacySimulatorData_->inputrec->etc == etcVRESCALE)
+ // 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)
{
- // TODO: With increased complexity of the propagator, this will need further development,
- // e.g. using propagators templated for velocity propagation policies and a builder
- propagatorVelocitiesAndPositions->setNumVelocityScalingVariables(
- legacySimulatorData_->inputrec->opts.ngtc);
- auto thermostat = std::make_unique<VRescaleThermostat>(
- legacySimulatorData_->inputrec->nsttcouple, 0, true,
- legacySimulatorData_->inputrec->ld_seed, legacySimulatorData_->inputrec->opts.ngtc,
- legacySimulatorData_->inputrec->delta_t * legacySimulatorData_->inputrec->nsttcouple,
- legacySimulatorData_->inputrec->opts.ref_t,
- legacySimulatorData_->inputrec->opts.tau_t, legacySimulatorData_->inputrec->opts.nrdf,
- energyDataPtr, propagatorVelocitiesAndPositions->viewOnVelocityScaling(),
- propagatorVelocitiesAndPositions->velocityScalingCallback(),
- legacySimulatorData_->state_global, legacySimulatorData_->cr,
- legacySimulatorData_->inputrec->bContinuation);
- checkpointClients->emplace_back(thermostat.get());
- energyDataPtr->setVRescaleThermostat(thermostat.get());
- addToCallListAndMove(std::move(thermostat), elementCallList, elementsOwnershipList);
+ builder->add<Propagator<IntegrationStage::ScalePositions>>(
+ PropagatorTag("ScaleMTTKXPost"));
}
- addToCallListAndMove(std::move(propagatorVelocitiesAndPositions), elementCallList,
- elementsOwnershipList);
if (legacySimulatorData_->constr)
{
- auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Positions>>(
- legacySimulatorData_->constr, statePropagatorDataPtr, energyDataPtr,
- freeEnergyPerturbationDataPtr, MASTER(legacySimulatorData_->cr),
- legacySimulatorData_->fplog, legacySimulatorData_->inputrec,
- legacySimulatorData_->mdAtoms->mdatoms());
- energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
- trajectorySignallerBuilder->registerSignallerClient(
- compat::make_not_null(constraintElement.get()));
- loggingSignallerBuilder->registerSignallerClient(
- compat::make_not_null(constraintElement.get()));
+ builder->add<ConstraintsElement<ConstraintVariable::Positions>>();
+ }
+
+ if (legacySimulatorData_->inputrec->bPull)
+ {
+ builder->add<PullElement>();
+ }
+
+ builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>();
- addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
+ // Propagate box from t to t+dt
+ if (legacySimulatorData_->inputrec->epc == PressureCoupling::Mttk)
+ {
+ builder->add<MttkBoxScaling>(mttkPropagatorConnectionDetails);
}
- addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
- auto energyElement = compat::make_not_null(energyDataPtr->element());
- trajectoryElementBuilder->registerWriterClient(energyElement);
- trajectorySignallerBuilder->registerSignallerClient(energyElement);
- energySignallerBuilder->registerSignallerClient(energyElement);
- checkpointClients->emplace_back(energyElement);
- // we have the energies at time t here!
- addToCallList(energyElement, elementCallList);
- if (prBarostat)
+ else if (legacySimulatorData_->inputrec->epc == PressureCoupling::CRescale)
{
- addToCallListAndMove(std::move(prBarostat), elementCallList, elementsOwnershipList);
+ // Legacy implementation allows combination of C-Rescale with Trotter Nose-Hoover
+ builder->add<FirstOrderPressureCoupling>(0, ReportPreviousStepConservedEnergy::Yes);
}
}
else
{
gmx_fatal(FARGS, "Integrator not implemented for the modular simulator.");
}
-
- auto integrator = std::make_unique<CompositeSimulatorElement>(std::move(elementCallList),
- std::move(elementsOwnershipList));
- // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
- return std::move(integrator);
+ builder->add<EnergyData::Element>();
}
bool ModularSimulator::isInputCompatible(bool exitOnFailure,
return condition;
};
- bool isInputCompatible = true;
-
// GMX_USE_MODULAR_SIMULATOR allows to use modular simulator also for non-standard uses,
// such as the leap-frog integrator
const auto modularSimulatorExplicitlyTurnedOn = (getenv("GMX_USE_MODULAR_SIMULATOR") != nullptr);
"or unset both to recover default behavior.");
GMX_RELEASE_ASSERT(
- !(modularSimulatorExplicitlyTurnedOff && inputrec->eI == eiVV
- && inputrec->epc == epcPARRINELLORAHMAN),
+ !(modularSimulatorExplicitlyTurnedOff && inputrec->eI == IntegrationAlgorithm::VV
+ && inputrec->epc == PressureCoupling::ParrinelloRahman),
"Cannot use a Parrinello-Rahman barostat with md-vv and "
"GMX_DISABLE_MODULAR_SIMULATOR=ON, "
"as the Parrinello-Rahman barostat is not implemented in the legacy simulator. Unset "
"GMX_DISABLE_MODULAR_SIMULATOR or use a different pressure control algorithm.");
- isInputCompatible =
- isInputCompatible
- && conditionalAssert(
- inputrec->eI == eiMD || inputrec->eI == eiVV,
- "Only integrators md and md-vv are supported by the modular simulator.");
+ bool isInputCompatible = conditionalAssert(
+ inputrec->eI == IntegrationAlgorithm::MD || inputrec->eI == IntegrationAlgorithm::VV,
+ "Only integrators md and md-vv are supported by the modular simulator.");
isInputCompatible = isInputCompatible
- && conditionalAssert(inputrec->eI != eiMD || modularSimulatorExplicitlyTurnedOn,
+ && conditionalAssert(inputrec->eI != IntegrationAlgorithm::MD
+ || modularSimulatorExplicitlyTurnedOn,
"Set GMX_USE_MODULAR_SIMULATOR=ON to use the modular "
"simulator with integrator md.");
- isInputCompatible =
- isInputCompatible
- && conditionalAssert(!doRerun, "Rerun is not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
&& conditionalAssert(
- inputrec->etc == etcNO || inputrec->etc == etcVRESCALE,
- "Only v-rescale thermostat is supported by the modular simulator.");
+ !inputrec->useMts,
+ "Multiple time stepping is not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
- && conditionalAssert(
- inputrec->epc == epcNO || inputrec->epc == epcPARRINELLORAHMAN,
- "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 == efepNO || inputrec->efep == efepYES
- || inputrec->efep == efepSLOWGROWTH,
- "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.");
+ && conditionalAssert(!doRerun, "Rerun is not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
- && conditionalAssert(inputrec->opts.ngacc == 1 && inputrec->opts.acc[0][XX] == 0.0
- && inputrec->opts.acc[0][YY] == 0.0
- && inputrec->opts.acc[0][ZZ] == 0.0 && inputrec->cos_accel == 0.0,
+ && conditionalAssert(!inputrec->useConstantAcceleration && inputrec->cos_accel == 0.0,
"Acceleration is not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
- && conditionalAssert(inputrec->opts.ngfrz == 1 && inputrec->opts.nFreeze[0][XX] == 0
- && inputrec->opts.nFreeze[0][YY] == 0
- && inputrec->opts.nFreeze[0][ZZ] == 0,
+ && conditionalAssert(!inputrecFrozenAtoms(inputrec),
"Freeze groups are not supported by the modular simulator.");
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,
"Essential dynamics is not supported by the modular simulator.");
isInputCompatible = isInputCompatible
- && conditionalAssert(inputrec->eSwapCoords == eswapNO,
+ && conditionalAssert(inputrec->eSwapCoords == SwapType::No,
"Ion / water position swapping is not supported by "
"the modular simulator.");
isInputCompatible =
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.");
-
+ 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;
}
-ModularSimulator::ModularSimulator(std::unique_ptr<LegacySimulatorData> legacySimulatorData) :
- legacySimulatorData_(std::move(legacySimulatorData))
+ModularSimulator::ModularSimulator(std::unique_ptr<LegacySimulatorData> legacySimulatorData,
+ std::unique_ptr<ReadCheckpointDataHolder> checkpointDataHolder) :
+ legacySimulatorData_(std::move(legacySimulatorData)),
+ checkpointDataHolder_(std::move(checkpointDataHolder))
{
checkInputForDisabledFunctionality();
}
void ModularSimulator::checkInputForDisabledFunctionality()
{
- isInputCompatible(true, legacySimulatorData_->inputrec,
- legacySimulatorData_->mdrunOptions.rerun, *legacySimulatorData_->top_global,
- legacySimulatorData_->ms, legacySimulatorData_->replExParams,
- &legacySimulatorData_->fr->listedForces->fcdata(),
+ isInputCompatible(true,
+ legacySimulatorData_->inputrec,
+ legacySimulatorData_->mdrunOptions.rerun,
+ legacySimulatorData_->top_global,
+ legacySimulatorData_->ms,
+ legacySimulatorData_->replExParams,
+ legacySimulatorData_->fr->fcdata.get(),
opt2bSet("-ei", legacySimulatorData_->nfile, legacySimulatorData_->fnm),
legacySimulatorData_->membed != nullptr);
if (legacySimulatorData_->observablesHistory->edsamHistory)
}
}
+void ModularSimulator::readCheckpointToTrxFrame(t_trxframe* fr,
+ ReadCheckpointDataHolder* readCheckpointDataHolder,
+ const CheckpointHeaderContents& checkpointHeaderContents)
+{
+ GMX_RELEASE_ASSERT(checkpointHeaderContents.isModularSimulatorCheckpoint,
+ "ModularSimulator::readCheckpointToTrxFrame can only read checkpoints "
+ "written by modular simulator.");
+ fr->bStep = true;
+ fr->step = int64_to_int(checkpointHeaderContents.step, "conversion of checkpoint to trajectory");
+ fr->bTime = true;
+ fr->time = checkpointHeaderContents.t;
+
+ fr->bAtoms = false;
+
+ StatePropagatorData::readCheckpointToTrxFrame(
+ fr, readCheckpointDataHolder->checkpointData(StatePropagatorData::checkpointID()));
+ if (readCheckpointDataHolder->keyExists(FreeEnergyPerturbationData::checkpointID()))
+ {
+ FreeEnergyPerturbationData::readCheckpointToTrxFrame(
+ fr, readCheckpointDataHolder->checkpointData(FreeEnergyPerturbationData::checkpointID()));
+ }
+ else
+ {
+ FreeEnergyPerturbationData::readCheckpointToTrxFrame(fr, std::nullopt);
+ }
+}
+
} // namespace gmx
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff