asyncmd.trajectory.InPartsTrajectoryPropagator#
Useful for making efficient use of backfilling and/or running simulations that are longer than the time-limit, i.e. when using slurm or another queuing system. Use it together with e.g. the SlurmGmxEngine and not locally like here done for demonstration purposes.
Imports and some basic checks that everything is available#
%%bash
# if using the module system to make gromacs and friends available:
# check that they are loaded!
#module list
%%bash
# unix only, check that gmx is available
which gmx
/usr/local/gromacs-2022.4/bin/gmx
%matplotlib inline
import os
import asyncio
import asyncmd
from asyncmd import gromacs as asyncgmx
from asyncmd import trajectory as asynctraj
Could not initialize SLURM cluster handling. If you are sure SLURM (sinfo/sacct/etc) is available try calling `asyncmd.config.set_slurm_settings()` with the appropriate arguments.
Setup working directory#
We will write the trajectory output to it.
workdir = "."
Load two different configurations as asyncmd.Trajectory#
We will use them as starting configurations.
# create an asyncmd.Trajectory of the initial configuration from the `GmxEngine.ipynb` notebook
conf_in_alphaR = asyncmd.Trajectory(trajectory_files="../resources/gromacs/capped_alanine_dipeptide/conf_in_alphaR.trr",
structure_file="../resources/gromacs/capped_alanine_dipeptide/conf.gro",
)
# create a second asyncmd.Trajectory of another configuration (in another state)
conf_in_C7eq = asyncmd.Trajectory(trajectory_files="../resources/gromacs/capped_alanine_dipeptide/conf_in_C7eq.trr",
structure_file="../resources/gromacs/capped_alanine_dipeptide/conf.gro",
)
Load and potentially modify the parameter file (mdp file) for the molecular dynamics simulations#
# Pcoupl = C-rescale needs gromacs version >= 2021
mdp = asyncgmx.MDP("../resources/gromacs/capped_alanine_dipeptide/md.mdp")
print(mdp["Pcoupl"])
# set nstxout-compressed, such that the engines will produce XTC trajectories
mdp["nstxout-compressed"] = 20
# and deactivate trr trajectory output
mdp["nstxout"] = mdp["nstvout"] = 0
C-rescale
# if your gmx version is >= 2021 you should comment the next line since C-rescale give the correct ensemble (and Berendsen doesnt!)
#mdp["Pcoupl"] = "Berendsen"
Initialize the InPartsTrajectoryPropagators#
The InPartsTrajectoryPropagator produces a Trajectory of a given total length (n_steps) in parts of a given walltime (walltime_per_part). This is useful to make full use of SLURMs backfilling and also to generate Trajectories of a given total length that exceeds the queues timelimit.
# The walltime per trajectory part determines how long each of the parts of the trajectory will be
walltime = 10 / (60 * 60) # walltime is measured in hours, so this will be 180 s per part!
propas = [asynctraj.InPartsTrajectoryPropagator(
n_steps=2e4,
engine_cls=asyncgmx.GmxEngine,
engine_kwargs={"mdconfig": mdp,
"gro_file": "../resources/gromacs/capped_alanine_dipeptide/conf.gro",
"top_file": "../resources/gromacs/capped_alanine_dipeptide/topol_amber99sbildn.top",
"mdrun_extra_args": "-nt 2", # limit mdrun to 2 threads
},
walltime_per_part=walltime,
)
for _ in range(2)]
Use the propagate_and_concatenate() method which directly concatenates the generated trajectory parts into one trajectory#
The propagate method returns the list of trajectory parts and cut_and_concatenate can make the list into one continuous trajectory. The propagate_and_concatenate method just calls both of them in order for convenience.
# the `propagate_and_concatenate` method returns the concatenated trajectory of the requested (total) length
# Using asyncio.gather as usual to do both MD runs in parallel
wdir_alphaR = os.path.join(workdir, "from_alphaR")
os.mkdir(wdir_alphaR)
wdir_C7eq = os.path.join(workdir, "from_C7eq")
os.mkdir(wdir_C7eq)
traj_from_alphaR, traj_from_C7eq = await asyncio.gather(propas[0].propagate_and_concatenate(
starting_configuration=conf_in_alphaR,
workdir=wdir_alphaR,
deffnm="from_alphaR",
tra_out=os.path.join(wdir_alphaR, "traj_from_alphaR.xtc")
),
propas[1].propagate_and_concatenate(
starting_configuration=conf_in_C7eq,
workdir=wdir_C7eq,
deffnm="from_C7_eq",
tra_out=os.path.join(wdir_C7eq, "traj_from_C7_eq.xtc")
)
)
print(f"The trajectory from alphaR has {len(traj_from_alphaR)} frames, the one from C7_eq has {len(traj_from_C7eq)} frames.")
The trajectory from alphaR has 1001 frames, the one from C7_eq has 1001 frames.
You can easily extend your simulations#
Here we will just reset the number of steps for the existing propagator objects, but it would work the same if we would have initialized two new ones (using engine_cls and engine_kwargs compatible with our previous run, you can change the walltime).
# double the number of integration steps we want to do
propas[0].n_steps *= 2
propas[1].n_steps *= 2
# and run again, this time passing continuation=True and overwrite=True such that we overwrite the old concatenated trajectory
traj_from_alphaR, traj_from_C7eq = await asyncio.gather(propas[0].propagate_and_concatenate(
starting_configuration=conf_in_alphaR,
workdir=wdir_alphaR,
deffnm="from_alphaR",
tra_out=os.path.join(wdir_alphaR, "traj_from_alphaR.xtc"),
continuation=True,
overwrite=True,
),
propas[1].propagate_and_concatenate(
starting_configuration=conf_in_C7eq,
workdir=wdir_C7eq,
deffnm="from_C7_eq",
tra_out=os.path.join(wdir_C7eq, "traj_from_C7_eq.xtc"),
continuation=True,
overwrite=True,
)
)
print(f"The trajectory from alphaR has {len(traj_from_alphaR)} frames, the one from C7_eq has {len(traj_from_C7eq)} frames.")
The trajectory from alphaR has 2001 frames, the one from C7_eq has 2001 frames.