Gromacs issue #4571

New multi-sim protocol for MD in trajectory space (path sampling)

Papers:

Davide's Slides

Replica exchange, multisim, others, excluded from modularsimulator at this time. Ref: ModularSimulator::isInputCompatible It's probably not too exciting just to enable it and see what breaks :-) Mark can do this properly in January 2023, but not right now

Assumptions

  • one domain per replica (for simplicity)

Preconditions we must construct during setup

  • pair list searching needs a larger buffer so that the same pairlist can be used for all three force calculations, even though atoms will have extra displacements in the two extra force evaluations. For now, just add a big constant like 0.1 nm and later try to be more subtle, in order to improve performance.
  • modular simulator needs an object that has tasks to do the communication, so needs communicators that talk to the neighbouring replicas

The pipeline that we want to implement in modular simulator looks like

  1. Post receives for coordinates \(R^{\tau+1}\) and \(R^{\tau-1}\)
  2. Post receives for forces \(F^{\tau-1}\)
  3. Send \(R^{\tau}\) to left and right neighbours
  4. Compute \(F^{\tau}\) (communication runs in the background)
  5. Wait for \(R^{\tau+1}\) (send initiated in step 3.)
  6. Send \(F^{\tau}\) to right neighbour
  7. compute \(\eta^\tau\)
  8. compute \(F(R^{\tau}+\epsilon\eta^{\tau})\) (communication runs in the background)
  9. compute \(F(R^{\tau}-\epsilon\eta^{\tau})\)
  10. Wait for \(F^{\tau-1}\) (send initiated in step 6.) and \(R^{\tau-1}\) (send initiated in step 3.)
  11. compute \(G^{\tau}\)
  12. add external forces (e.g. AWH, pulling, or PLUMED)
  13. update and constraints and other normal things like checkpoint signalling, end-of-run signalling
  14. go to 1.

A note on program flow during simulation launch:

All of the ModularSimulator set-up starts after https://gitlab.com/gromacs/gromacs/-/blob/main/src/gromacs/mdrun/runner.cpp#L2200

The actual computation that will be performed for force calculation (by the existing ForceElement) during the MD step is extended by the contents of the MDModules container (not to be confused with the Modular Simulator framework!) which must be final before https://gitlab.com/gromacs/gromacs/-/blob/main/src/gromacs/mdrun/runner.cpp#L2218

See also https://manual.gromacs.org/current/doxygen/html-full/page_mdmodules.xhtml, but Mark is proposing that you try to defer integrating with these frameworks right away, and try to use existing facilities and/or localized tightly-scoped changes for handling input parameters, new communication calls, and new blocks of computation.

A multi-sim implementation for modular simulator exists: https://gitlab.com/gromacs/gromacs/-/commit/84a5f710734aed6725ee15cae6565a8e5e2110ed. Mark will put that on main branch to use as a starting point for this work.

Pre-conditions have been met, the integration of path MD algorithm in Gromacs is underway.

We still need a way to output the effective potential energy of the system (e.g. in the .trr?).

Hints for implementing multi-domain support

MPI tests e.g. https://gitlab.com/gromacs/gromacs/-/blob/main/src/gromacs/domdec/tests/CMakeLists.txt and https://gitlab.com/gromacs/gromacs/-/blob/main/src/gromacs/domdec/tests/haloexchange_mpi.cpp

repartitioning happens in https://gitlab.com/gromacs/gromacs/-/blob/main/src/gromacs/modularsimulator/domdechelper.cpp, perhaps you want to register a IDomDecHelperClient to get called back each time repartitioning happens

Possible workflow for that callback:

nearest neighbour means neighbours plus self, for brevity

Having previously over-allocated the receive buffers for global atom indices during construction,

  • post MPI_Irecv for vectors of global atom indices A' for all nearest neighbours from the adjacent replicas
  • look up the global atom indices for all particles in my domain (e.g. from ForceElement->cr_->dd->globalAtomIndices) -> make a vector A of them
  • send A to each nearest neighbour in each adjacent replica
  • MPI_Waitall on the receives of A', passing an array of matching MPI_Status objects
  • Use the status.COUNT field (I think that's the name) to find out the logical size of the vector from each neighbour

Now this domain knows which nearest neighbour has each particle that it cares about.

  • Post MPI_Irecv for the incoming vectors B' for particles it will have to communicate each step
  • Construct vectors B to return to each neighbour describing global indices of atoms that it cares about in future
  • MPI_Waitall on the receives for B', again using MPI_Status to find out the lengths of the messages

Now this domain knows which particles are expected by each of its nearest neighbours.

  • turn that back into vectors of local particle indices C using stuff in src/gromacs/domdec/ga2la.h

On each MD step, after the normal force calculation

  • post MPI_Irecv for vectors of coordinates (NB 3 floats per particle) D'
  • Fill coordinate vector D for each nearest neighbour using C
  • send D to each neighbour in adjacent replicas
  • MPI_Waitall on the receives for D'
  • Compute eta using D'
  • Update local coordinates using eta
  • Compute new force
  • Compute local coordinates again
  • Compute new force

Advanced technique (only do this when everything is working!) post the MPI_Irecv for the next step immediately after the MPI_Waitall. This helps the MPI library be efficient when some domains run faster than others (which they always will). But the book keeping is a bit trickier.

Select a repo