AMBER Constant pH MD Setup

What this workflow does

This workflow sets up and runs a constant pH molecular dynamics simulation with AmberTools. The example system is Bovine Pancreatic Trypsin Inhibitor (BPTI, PDB 6PTI).

A standard MD run fixes the protonation state of every residue for the whole trajectory. A constant pH run does not. It lets titratable residues change their protonation state during the simulation, at a fixed solvent pH. The result is a predicted pKa per residue.

The workflow prepares the structure with pdb4amber. LEaP builds the topology with the ff14SB force field and the constph residue templates. LEaP then solvates the system in a truncated octahedron box of TIP3P water and adds ions.

cpinutil generates the constant pH input file. It marks the titratable residues: ASP, GLU, CYS, LYS and TYR. sander then runs the minimization with backbone restraints, the heating to 300 K, the NVT equilibration and the NPT equilibration.

The production stage runs the constant pH MD at solvent pH 7.0. cphstats reads the protonation record and predicts a pKa for each titratable residue. It also reports the population of each protonation state.

The compute problem

Constant pH MD costs more than plain MD. sander performs a Monte Carlo protonation-state attempt every ntcnstph steps. Each attempt evaluates a Generalized Born energy. The production stage is therefore the long pole, and it gets longer with the number of titratable residues.

BPTI is small. A realistic target is not. When you scale to a larger protein, or when you scan several solvent pH values, the production stage needs a cluster. The setup stages still take seconds.

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

How Horus solves it

Horus runs the AMBER stages in the quay.io/biocontainers/biobb_amber image and the biobb_io fetch stage in the native conda environment. The split is declarative. Each stage names its executor with one field.

A pH scan is the case that pays. Copy the sander_cph_free stage once per pH value and change solvph in each config file. Point every copy at a remote executor. Horus runs them in parallel and collects every cpout file. The topology, the cpin file and the equilibrated coordinates are built once and shared.

The equilibration stages are a different profile again. They are serial and they must finish in order before any production stage starts. Keep them on one host. Horus enforces the order from the data dependencies.

The alternative is a shell script per pH value on the cluster, plus manual bookkeeping of which cpin file matches which trajectory. The cphstats analysis then needs the matching pair. Horus keeps the pairing in the graph.

Pipeline

fetch_pdb          Download 6PTI from the PDB (biobb_io)
   │
pdb4amber          Prepare structure for constant-pH AMBER
   │
leap_gen_top       Build protein topology (ff14SB + constph, leap)
   │
leap_solvate       Solvate in truncated octahedron box (TIP3P, leap)
   │
leap_add_ions      Neutralize and add ions (leap)
   │
parmed_cpinutil    Generate constant-pH input file (cpinutil, titratable residues)
   │
sander_min ──► process_min       Energy minimize with backbone restraints
   │
sander_heat ──► process_heat      Heat 0 → 300 K
   │
sander_nvt ──► process_nvt         NVT equilibration
   │
sander_npt ──► process_npt          NPT equilibration
   │
sander_cph_free    Constant pH production MD (sander, solvph 7.0)
   │
cphstats           Analyse protonation states, predict pKa values

Inputs and outputs

Inputs

This workflow has no file inputs. It fetches PDB 6PTI from the RCSB at run time. Set the code in configs/fetch_pdb.yaml.

Outputs land in workflow_results/results/:

• 6PTI.pdb, structure.pdb4amber.pdb: the fetched and prepared structures.
• structure.leap.top, structure.leap.crd: the ff14SB constph topology.
• structure.ions.parmtop: the solvated, ion-neutralized system.
• structure.cpin, structure.cpH.parmtop: the constant pH input and topology.
• sander.pH.netcdf, sander.pH.cpout: the production trajectory and the protonation record.
• cphstats.pH.dat: the predicted pKa per titratable residue.
• cphstats.pH.pop.dat: the protonation-state populations.

configs/parmed_cpinutil.yaml sets the Generalized Born model and the titratable residue list. configs/sander_*.yaml set solvph, icnstph and ntcnstph. configs/cphstats.yaml sets the running-average window.

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
• AMBER Constant pH MD tutorial