GROMACS will automatically use the modular simulator for the velocity
verlet integrator (`integrator = md-vv`), if the functionality chosen
in the other input parameters is implemented in the new framework.
-Currently, this includes NVE simulations, NVT simulations (
-`tcoupl = v-rescale` only), NPH simulation (`pcoupl = parrinello-rahman`
-only), and NPT simulations (`tcoupl = v-rescale` and
-`pcoupl = parrinello-rahman` only), with or without free energy perturbation.
+Currently, this includes NVE, NVT, NPH, and NPT simulations,
+with or without free energy perturbation, using thermodynamic
+boundary conditions
+
+* `tcoupl`: `no`, `v-rescale`, or `berendsen`
+* `pcoupl`: `no` or `parrinello-rahman`
To disable the modular simulator for cases defaulting to the new framework,
the environment variable `GMX_DISABLE_MODULAR_SIMULATOR=ON` can be set. To
order. Over time, with data dependencies clearly defined, this
approach can be modified to have independent tasks run in parallel.
-The approach is most easily displayed using some pseudo code:
-
- class ModularSimulator : public ISimulator
- {
- public:
- //! Run the simulator
- void run() override;
- private:
- std::vector<ISignaller*> signallers_;
- std::vector<ISimulatorElement*> elements_;
- std::queue<SimulatorRunFunction*> taskQueue_;
- }
-
- void ModularSimulator::run()
- {
- constructElementsAndSignallers();
- setupAllElements();
- while (not lastStep)
- {
- // Fill the task queue with new tasks (can be precomputed for many steps)
- populateTaskQueue();
- // Now simply loop through the queue and run one task after the next
- for (auto task : taskQueue)
- {
- (*task)(); // run task
- }
- }
- }
-
-This allows for an important division of tasks.
-
-* `constructElementsAndSignallers()` is responsible to **store the
- elements in the right order**. This includes the different order of
- element in different algorithms (e.g. leap-frog vs. velocity
- verlet), but also logical dependencies (energy output after compute
- globals).
-* `populateTaskQueue()` is responsible to **decide if elements need to
- run at a specific time step**. The elements get called in order, and
- decide whether they need to run at a specific step. This can be
- pre-computed for multiple steps. In the current implementation, the
- tasks are pre-computed for the entire life-time of the neighbor
- list.
-* **Running the actual simulation tasks** is done after the task queue
- was filled. This is achieved by simply looping over the task list,
- no conditionals or branching needed.
-
### Simulator elements
The task scheduler holds a list of *simulator elements*, defined by
called later during the actual (partial) simulation run. From the
point of view of the builder of the task scheduler, it is important to
note that the order of the elements determines the order in which
-computation is performed. The task scheduler periodically loops over
+computation is performed.
+
+ class ISimulatorElement
+ {
+ public:
+ /*! \\brief Query whether element wants to run at step / time
+ *
+ * Element can register one or more functions to be run at that step through
+ * the registration pointer.
+ */
+ virtual void scheduleTask(Step, Time, const RegisterRunFunction&) = 0;
+ //! Method guaranteed to be called after construction, before simulator run
+ virtual void elementSetup() = 0;
+ //! Method guaranteed to be called after simulator run, before deconstruction
+ virtual void elementTeardown() = 0;
+ //! Standard virtual destructor
+ virtual ~ISimulatorElement() = default;
+ };
+
+
+The task scheduler periodically loops over
its list of elements, builds a queue of function pointers to run, and
returns this list of tasks. As an example, a possible application
would be to build a new queue after each domain-decomposition (DD) /
register an energy calculation during that step, being ready to write
to trajectory when requested.
+ class ISignaller
+ {
+ public:
+ //! Function run before every step of scheduling
+ virtual void signal(Step, Time) = 0;
+ //! Method guaranteed to be called after construction, before simulator run
+ virtual void setup() = 0;
+ };
+
+### The modular simulator
+
+The approach is most easily displayed using some simplified (pseudo) code.
+
+The simulator itself is responsible to **store the elements in the
+right order** (in `addIntegrationElements`) This includes the different
+order of elements in different algorithms (e.g. leap-frog vs. velocity
+verlet), but also logical dependencies (energy output after compute
+globals). Once the algorithm has been built, the simulator simply
+executes one task after the next, until the end of the simulation is
+reached.
+
+ class ModularSimulator : public ISimulator
+ {
+ public:
+ //! Run the simulator
+ void run() override;
+ }
+
+ void ModularSimulator::run()
+ {
+
+ ModularSimulatorAlgorithmBuilder algorithmBuilder();
+ addIntegrationElements(&algorithmBuilder);
+ auto algorithm = algorithmBuilder.build();
+
+ while (const auto* task = algorithm.getNextTask())
+ {
+ // execute task
+ (*task)();
+ }
+ }
+
+The following snippet illustrates building a leap-frog integration
+algorithm. The algorithm builder allows for a concise description of
+the simulator algorithm.
+
+ void ModularSimulator::addIntegrationElements(ModularSimulatorAlgorithmBuilder* builder)
+ {
+ if (legacySimulatorData_->inputrec->eI == eiMD)
+ {
+ // The leap frog integration algorithm
+ builder->add<ForceElement>();
+ // We have a full state here (positions(t), velocities(t-dt/2), forces(t)
+ builder->add<StatePropagatorData::Element>();
+ if (legacySimulatorData_->inputrec->etc == etcVRESCALE)
+ {
+ builder->add<VRescaleThermostat>(-1, VRescaleThermostatUseFullStepKE::No);
+ }
+ builder->add<Propagator<IntegrationStep::LeapFrog>>(legacySimulatorData_->inputrec->delta_t,
+ RegisterWithThermostat::True,
+ RegisterWithBarostat::True);
+ if (legacySimulatorData_->constr)
+ {
+ builder->add<ConstraintsElement<ConstraintVariable::Positions>>();
+ }
+ builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>>();
+ // We have the energies at time t here
+ builder->add<EnergyData::Element>();
+ if (legacySimulatorData_->inputrec->epc == epcPARRINELLORAHMAN)
+ {
+ builder->add<ParrinelloRahmanBarostat>(-1);
+ }
+ }
+ }
+
+### The simulator algorithm
+
+The simulator algorithm is responsible to **decide if elements need to
+run at a specific time step**. The elements get called in order, and
+decide whether they need to run at a specific step. This can be
+pre-computed for multiple steps. In the current implementation, the
+tasks are pre-computed for the entire life-time of the neighbor
+list.
+
+The simulator algorithm offers functionality to get the next task
+from the queue. It owns all elements involved in the simulation
+and is hence controlling their lifetime. This ensures that pointers and
+callbacks exchanged between elements remain valid throughout the duration
+of the simulation run. It also maintains the list of tasks,
+and updates it when needed.
+
+ class ModularSimulatorAlgorithm
+ {
+ public:
+ //! Get next task in queue
+ [[nodiscard]] const SimulatorRunFunction* getNextTask();
+ private:
+ //! List of signalers
+ std::vector<std::unique_ptr<ISignaller>> signallerList_;
+ //! List of elements
+ std::vector<std::unique_ptr<ISimulatorElement>> elementsList_;
+
+ //! The run queue
+ std::vector<SimulatorRunFunction> taskQueue_;
+ //! The task iterator
+ std::vector<SimulatorRunFunction>::const_iterator taskIterator_;
+
+ //! Update task queue
+ void updateTaskQueue();
+ }
+
+The `getNextTask()` function is returning the next task in the task
+queue. It rebuilds the task list when needed.
+
+ const SimulatorRunFunction* ModularSimulatorAlgorithm::getNextTask()
+ {
+ if (!taskQueue_.empty())
+ {
+ taskIterator_++;
+ }
+ if (taskIterator_ == taskQueue_.end())
+ {
+ if (runFinished_)
+ {
+ return nullptr;
+ }
+ updateTaskQueue();
+ taskIterator_ = taskQueue_.begin();
+ }
+ return &*taskIterator_;
+ }
+
+Updating the task queue involves calling all signallers and
+elements for every step of the scheduling period. This refills
+the task queue. It is important to keep in mind that the *scheduling step* is not
+necessarily identical to the *current step* of the simulation. Most of
+the time, the scheduling step is ahead, as we are pre-scheduling steps.
+
+ void ModularSimulatorAlgorithm::updateTaskQueue()
+ {
+ for (Step schedulingStep = currentStep;
+ schedulingStep < currentStep + schedulingPeriod;
+ schedulingStep++)
+ {
+ Time time = getTime(schedulingStep);
+ // Have signallers signal any special treatment of scheduling step
+ for (const auto& signaller : signallerList)
+ {
+ signaller.signal(schedulingStep, time);
+ }
+ // Query all elements whether they need to run at scheduling step
+ for (const auto& element : signallerList)
+ {
+ element.schedule(schedulingStep, time, registerRunFunction_);
+ }
+ }
+ }
+
### Sequence diagrams
#### Pre-loop
\msc
hscale="2";
-ModularSimulator,
-Signallers [label="ModularSimulator::\nSignallers"],
-Elements [label="ModularSimulator::\nElements"],
-TaskQueue [label="ModularSimulator::\nTaskQueue"];
+ModularSimulatorBuilder [label="ModularSimulatorAlgorithmBuilder"],
+ModularSimulator [label="ModularSimulatorAlgorithm"],
+Signallers [label="ModularSimulatorAlgorithm::\nSignallers"],
+Elements [label="ModularSimulatorAlgorithm::\nElements"],
+TaskQueue [label="ModularSimulatorAlgorithm::\nTaskQueue"];
--- [ label = "constructElementsAndSignallers()" ];
- ModularSimulator => Signallers [ label = "Create signallers\nand order them" ];
- ModularSimulator => Elements [ label = "Create elements\nand order them" ];
+ ModularSimulatorBuilder => Signallers [ label = "Create signallers\nand order them" ];
+ ModularSimulatorBuilder => Elements [ label = "Create elements\nand order them" ];
--- [ label = "constructElementsAndSignallers()" ];
|||;
|||;
\msc
hscale="2";
-ModularSimulator,
-Signallers [label="ModularSimulator::\nSignallers"],
-Elements [label="ModularSimulator::\nElements"],
-TaskQueue [label="ModularSimulator::\nTaskQueue"];
+ModularSimulator [label="ModularSimulatorAlgorithm"],
+Signallers [label="ModularSimulatorAlgorithm::\nSignallers"],
+Elements [label="ModularSimulatorAlgorithm::\nElements"],
+TaskQueue [label="ModularSimulatorAlgorithm::\nTaskQueue"];
ModularSimulator box TaskQueue [ label = "loop: while(not lastStep)" ];
ModularSimulator note TaskQueue [ label = "The task queue is empty. The simulation state is at step N.", textbgcolor="yellow" ];
\msc
hscale="2";
-ModularSimulator,
-Signallers [label="ModularSimulator::\nSignallers"],
-Elements [label="ModularSimulator::\nElements"],
-TaskQueue [label="ModularSimulator::\nTaskQueue"];
+ModularSimulator [label="ModularSimulatorAlgorithm"],
+Signallers [label="ModularSimulatorAlgorithm::\nSignallers"],
+Elements [label="ModularSimulatorAlgorithm::\nElements"],
+TaskQueue [label="ModularSimulatorAlgorithm::\nTaskQueue"];
--- [ label = "populateTaskQueue()" ];
ModularSimulator box ModularSimulator [ label = "doDomainDecomposition()\ndoPmeLoadBalancing()" ];
## Acceptance tests and further plans
-In January 2019, we defined acceptance tests which need to be
+Acceptance tests which need to be
fulfilled to make the modular simulator the default code path:
* End-to-end tests pass on both `do_md` and the new loop in
- Jenkins pre- and post-submit matrices
+ Gitlab pre- and post-submit pipelines
* Physical validation cases pass on the new loop
* Performance on different sized benchmark cases, x86 CPU-only
and NVIDIA GPU are at most 1% slower -
* simultaneous execution of independent tasks
We will probably not prioritize support for (and might consider
-deprecating from do_md for GROMACS 2020)
+deprecating from do_md in a future GROMACS version)
* Simulated annealing
* REMD
* Simulated tempering
* Time-averaged restraints
* Freeze, deform, cos-acceleration
-## Signallers and elements
+## Signaller and element details
The current implementation of the modular simulator consists of
the following signallers and elements:
* `EnergySignaller`: Informs its clients about energy related
special steps, namely energy calculation steps, virial
calculation steps, and free energy calculation steps.
-* `TrajectoryElement`: Informs its clients if writing to
+* `TrajectorySignaller`: Informs its clients if writing to
trajectory (state [x/v/f] and/or energy) is planned for the
- current step. Note that the `TrajectoryElement` is not a
- pure signaller, but also implements the `ISimulatorElement`
- interface (see section "Simulator Elements" below).
+ current step.
### Simulator Elements
#### `TrajectoryElement`
-The `TrajectoryElement` is a special element, as it
-is both implementing the `ISimulatorElement` and the `ISignaller`
-interfaces. During the signaller phase, it is signalling its
-_signaller clients_ that the trajectory will be written at the
-end of the current step. During the simulator run phase, it is
-calling its _trajectory clients_ (which do not necessarily need
-to be identical with the signaller clients), passing them a valid
+The `TrajectoryElement` is calling its trajectory clients, passing them a valid
output pointer and letting them write to trajectory. Unlike the
legacy implementation, the trajectory element itself knows nothing
about the data that is written to file - it is only responsible
to inform clients about trajectory steps, and providing a valid
file pointer to the objects that need to write to trajectory.
-#### `StatePropagatorData::Element`
-The `StatePropagatorData::Element` takes part in the simulator run, as it might
-have to save a valid state at the right moment during the
-integration. Placing the StatePropagatorData correctly is for now the
-duty of the simulator builder - this might be automated later
-if we have enough meta-data of the variables (i.e., if
-`StatePropagatorData` knows at which time the variables currently are,
-and can decide when a valid state (full-time step of all
-variables) is reached. The `StatePropagatorData::Element` is also a client of
-both the trajectory signaller and writer - it will save a
-state for later writeout during the simulator step if it
-knows that trajectory writing will occur later in the step,
-and it knows how to write to file given a file pointer by
-the `TrajectoryElement`.
-
-#### `EnergyData::Element`
-The `EnergyData::Element` takes part in the simulator run, as it
-does either add data (at energy calculation steps), or
-record a non-calculation step (all other steps). It is the
-responsibility of the simulator builder to ensure that the
-`EnergyData::Element` is called at a point of the simulator run
-at which it has access to a valid energy state.
-
#### `ComputeGlobalsElement`
The `ComputeGlobalsElement` encapsulates the legacy calls to
`compute_globals`. While a new approach to the global reduction
calls to `do_force` and `do_shellfc`, respectively. It is the
responsibility of the simulator builder to place them at the right
place of the integration algorithm. Moving forward, a version of these
-elements which would allow to only calculate forces of subsets of
-degrees of freedom would be desirable to pave the way towards multiple
-time step integrators, allowing to integrate slower degrees of freedom
+elements which would allow calling of do_force with subsets of the topology
+would be desirable to pave the way towards multiple time step integrators
+within modular simulator, allowing to integrate slower degrees of freedom
at a different frequency than faster degrees of freedom.
#### `ConstraintElement`
callback to the propagator to update the velocity scaling factor, and
scales the box and the positions of the system.
+#### `StatePropagatorData::Element`
+The `StatePropagatorData::Element` takes part in the simulator run, as it might
+have to save a valid state at the right moment during the
+integration. Placing the StatePropagatorData correctly is for now the
+duty of the simulator builder - this might be automated later
+if we have enough meta-data of the variables (i.e., if
+`StatePropagatorData` knows at which time the variables currently are,
+and can decide when a valid state (full-time step of all
+variables) is reached. The `StatePropagatorData::Element` is also a client of
+both the trajectory signaller and writer - it will save a
+state for later writeout during the simulator step if it
+knows that trajectory writing will occur later in the step,
+and it knows how to write to file given a file pointer by
+the `TrajectoryElement`.
+
+#### `EnergyData::Element`
+The `EnergyData::Element` takes part in the simulator run,
+either adding data (at energy calculation steps), or
+recording a non-calculation step (all other steps). It is the
+responsibility of the simulator builder to ensure that the
+`EnergyData::Element` is called at a point of the simulator run
+at which it has access to a valid energy state.
+
+It subscribes to the trajectory signaller, the energy signaller,
+and the logging signaller to know when an energy calculation is
+needed and when a non-recording step is enough. The EnergyData
+element is also a subscriber to the trajectory writer element,
+as it is responsible to write energy data to trajectory.
+
#### `FreeEnergyPerturbationData::Element`
The `FreeEnergyPerturbationData::Element` is a member class of
`FreeEnergyPerturbationData` that updates the lambda
Note that the `StatePropagatorData` can be converted to and from the
legacy `t_state` object. This is useful when dealing with
functionality which has not yet been adapted to use the new
-data approach - of the elements currently implemented, only
+data approach. Of the elements currently implemented, only
domain decomposition, PME load balancing, and the initial
constraining are using this.
### `EnergyData`
The EnergyData owns the EnergyObject, and is hence responsible
-for saving energy data and writing it to trajectory. It also owns
-the tensors for the different virials and the pressure as well as
-the total dipole vector.
-
-It subscribes to the trajectory signaller, the energy signaller,
-and the logging signaller to know when an energy calculation is
-needed and when a non-recording step is enough. The simulator
-builder is responsible to place the element in a location at
-which a valid energy state is available. The EnergyData is
-also a subscriber to the trajectory writer element, as it is
-responsible to write energy data to trajectory.
-
+for saving energy data and writing it to trajectory.
The EnergyData offers an interface to add virial contributions,
but also allows access to the raw pointers to tensor data, the
dipole vector, and the legacy energy data structures.
current FEP state, offering access to its values via getter
functions.
-### `TopologyHolder`
-The topology object owns the local topology and holds a constant reference
-to the global topology owned by the ISimulator.
-
-The local topology is only infrequently changed if domain decomposition is
-on, and never otherwise. The topology holder therefore offers elements to register
-as ITopologyHolderClients. If they do so, they get a handle to the updated local
-topology whenever it is changed, and can rely that their handle is valid
-until the next update. The domain decomposition element is defined as friend
-class to be able to update the local topology when needed.
+## Simulator algorithm builder
+Elements that define the integration algorithm (i.e. which are
+added using the templated `ModularSimulatorAlgorithmBuilder::add`
+method) need to implement a `getElementPointerImpl` factory function.
+This gives them access to the data structures and some other
+infrastructure, but also allows elements to accept additional
+arguments (e.g frequency, offset, ...).
-The topology holder is not a `ISimulatorElement`, i.e. it does not take part in the
-simulator loop.
+ template<typename Element, typename... Args>
+ void ModularSimulatorAlgorithmBuilder::add(Args&&... args)
+ {
+ // Get element from factory method
+ auto* element = static_cast<Element*>(getElementPointer<Element>(
+ legacySimulatorData_, &elementAdditionHelper_, statePropagatorData_.get(),
+ energyData_.get(), freeEnergyPerturbationData_.get(), &globalCommunicationHelper_,
+ std::forward<Args>(args)...));
+
+ // Make sure returned element pointer is owned by *this
+ // Ensuring this makes sure we can control the life time
+ if (!elementExists(element))
+ {
+ throw ElementNotFoundError("Tried to append non-existing element to call list.");
+ }
+
+ // Register element with infrastructure
+ }
+
+Note that `getElementPointer<Element>` will call `Element::getElementPointerImpl`,
+which needs to be implemented by the different elements.
## Infrastructure
### `DomDecHelper` and `PmeLoadBalanceHelper`
steps. This allows elements to receive the new topology / state before
deciding what functionality they need to run.
-### `Checkpointing`
+### Checkpointing
The `CheckpointHelper` is responsible to write checkpoints, and to offer
its clients access to the data read from checkpoint.
that checkpointing is about to occur. At the next NS step, the checkpoint
is written. On the last step, checkpointing happens immediately before the
step (no signalling). To be able to react to last step being signalled,
-the CheckpointHelper does also implement the `ISimulatorElement` interface,
-but does only register a function if the last step has been called.
+the CheckpointHelper also implements the `ISimulatorElement` interface,
+but only registers a function if the last step has been called.
Checkpointing happens at the top of a simulation step, which gives a
straightforward re-entry point at the top of the simulator loop.
edit the respective `do_cpt_XXX` function in `checkpoint.cpp` which interacts
with the XDR library. The alternative would be to write an entirely new data
structure, changing the function signature of all checkpoint-related functions,
-and write a corresponding low-level routing interacting with the XDR library.
+and write a corresponding low-level routine interacting with the XDR library.
**The MdModule approach:** To allow for modules to write checkpoints, the legacy
checkpoint was extended by a KVTree. When writing to checkpoint, this tree gets
the KVTree gets deserialized, and then distributed to the modules which can read
back the data previously stored.
-**The modular simulator approach in GROMACS 2020:** To allow for checkpointing
-before a truly modularized checkpoint infrastructure existed, modular simulator
-as implemented in GROMACS 2020 passed a t_state object to the checkpoint clients,
-which they could fill with their data. This is obviously not a permanent solution.
-
##### Modular simulator design
The MdModule checks off almost all requirements to a modularized checkpointing format.
* Hide the implementation details of the data structure holding the data (currently,
a KV-Tree) from the clients. This allows to change the implementation details of
reading / writing checkpoints without touching client code.
-* Offer a unified way to read and write to data, allowing clients write one
+* Offer a unified way to read and write to data, allowing clients to write one
(templated) function to read to and write from checkpoint. This allows to
- eliminate code duplication and the danger of having read and write functions
- getting out of sync.
+ eliminate code duplication and the danger of having read and write functions
+ getting out of sync.
The modular simulator checkpointing does not currently change the way that the
legacy simulator is checkpointing. Some data structures involved in the legacy
the help of InMemorySerializer, a KV-Tree could likely be emulated (with sub-objects
that serialize to memory, and then a parent object that serializes this memory to
file), but that doesn't present a clear advantage anymore.
+
+### `TopologyHolder`
+The topology object owns the local topology and holds a constant reference
+to the global topology owned by the ISimulator.
+
+The local topology is only infrequently changed if domain decomposition is
+on, and never otherwise. The topology holder therefore offers elements to register
+as ITopologyHolderClients. If they do so, they get a handle to the updated local
+topology whenever it is changed, and can rely that their handle is valid
+until the next update. The domain decomposition element is defined as friend
+class to be able to update the local topology when needed.
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff