ABC DNA MD Setup

What this workflow does

This workflow implements the Ascona B-DNA Consortium (ABC) structure-preparation protocol with AmberTools. The ABC protocol gives every consortium member one reproducible recipe for DNA MD setup. The same recipe makes the trajectories comparable across laboratories.

The example system is the Drew-Dickerson Dodecamer, sequence CGCGAATTCGCG.

LEaP builds the DNA topology with the parmbsc1 force field. LEaP then solvates the system in a truncated octahedron box of SPC/E water. It neutralizes the charge with potassium and adds KCl to 150 mM.

cpptraj randomizes the ion positions. This removes the placement bias that LEaP introduces. parmed then applies hydrogen-mass repartitioning. This lets the production run use a 4 fs timestep instead of 2 fs.

sander runs the ten-stage Roe and Brooks equilibration protocol. The first stages minimize and heat under strong restraints. The restraints relax over the later stages. The last stages run NPT without restraints. A free production MD run follows at the 4 fs timestep.

The compute problem

The ABC protocol is a long serial chain. Ten equilibration stages run in order. Each one reads the restart file of the one before it. No stage overlaps another.

The chain is also brittle in a specific way. If stage seven fails, you do not want to repeat stages one to six. Those stages already produced correct restart files. A shell script that runs all ten commands does repeat them.

The cost profile is skewed. The minimization stages finish in minutes. The NPT stages and the production run need a GPU. The production run at 4 fs is the stage you extend to microseconds.

The biobb_amber package has no osx-arm64 conda build. Every LEaP stage, every sander stage and the parmed stage need a Linux container on Apple Silicon.

How Horus solves it

Horus tracks each equilibration stage as its own node with its own outputs. A rerun starts at the first stage whose inputs changed. The earlier restart files stay valid. You edit step7.in and Horus repeats stages seven to ten only.

The executor routing follows the cost. Run the LEaP stages, the cpptraj ion randomization and the parmed repartitioning locally. They are cheap and they fix the system definition. Point sander_eq6 through sander_md at a GPU host with one executor: field each.

Horus moves the repartitioned topology and the current restart file to that host. It brings each trajectory back. The ten sander command strings never change. The same workflow file therefore runs end to end on a workstation for a test, and across two machines for the real run.

This matters for reproducibility, which is the point of the ABC protocol. The workflow file is the protocol. A consortium member reads one YAML file and sees every parameter. A member on different hardware changes only the executor fields.

Pipeline

leap_gen_top              Build DNA topology (parmbsc1 force field, leap)
   │
leap_solvate               Solvate in truncated octahedron box (SPC/E water, leap)
   │
leap_add_ions               Neutralize + 100 mM ion concentration (leap)
   │
cpptraj_randomize_ions        Randomize ion positions (cpptraj)
   │
parmed_hmassrepartition         Hydrogen-mass repartition topology for 4fs timestep (parmed)
   │
sander_eq1  ──► process_eq1       Eq step 1: minimization with restraints
   │
sander_eq2  ──► process_eq2       Eq step 2: heating with restraints
   │
sander_eq3  ──► process_eq3       Eq step 3: minimization with restraints
   │
sander_eq4  ──► process_eq4       Eq step 4: minimization with restraints
   │
sander_eq5  ──► process_eq5       Eq step 5: minimization without restraints
   │
sander_eq6  ──► process_eq6       Eq step 6: NPT with restraints
   │
sander_eq7  ──► process_eq7       Eq step 7: NPT with restraints
   │
sander_eq8  ──► process_eq8       Eq step 8: NPT with backbone restraints
   │
sander_eq9  ──► process_eq9       Eq step 9: NPT without restraints
   │
sander_eq10 ──► process_eq10      Eq step 10: NPT production equilibration
   │
sander_md                  Free production MD (4 fs timestep)

Inputs and outputs

Inputs

This workflow uses real input data checked into the directory. It fetches nothing at run time.

• CGCGAATTCGCG.pdb: the starting Drew-Dickerson Dodecamer structure.
• ABCix_config_files/step1.in … step10.in: the ten sander MDIN files of the ABC equilibration protocol.
• ABCix_config_files/md.in: the production MD sander input file.

Outputs land in workflow_results/results/:

• structure.leap.top, structure.leap.crd: the parmbsc1 DNA topology.
• structure.solv.parmtop, structure.ions.parmtop: the solvated and the neutralized system.
• structure.randIons.crd: the system with randomized ion placement.
• structure.leap.4fs.top: the hydrogen-mass-repartitioned topology.
• sander.eq1.nc … sander.eq10.nc: the ten equilibration trajectories, each with its log, mdinfo, restart and energy files.
• sander.md.nc: the free production MD trajectory at 4 fs.

The configs/leap_*.yaml files hold the force field, the water model, the box geometry and the ion concentration. The configs/sander_eq*.yaml files select which MDIN file each stage runs.

Run the workflow

Install the horus-runtime one time:

uv sync

If you do not have uv, install it first:

curl -LsSf https://astral.sh/uv/install.sh | sh

You can also install the packages with pip:

pip install horus-runtime horus-environments

Then run the workflow:

uv run horus run workflow.yaml

References

• BioExcel Building Blocks (biobb)
• biobb_amber
• AmberTools
• Ascona B-DNA Consortium / NAFlex ABC
• The static and dynamic structural heterogeneities of B-DNA