Run an MD simulation with ACEMD#

You will learn: how to take a built, solvated, ionised system and drive it through a complete MD protocol - equilibration then production - using ACEMD.

Prerequisites:

HTMD installed.
ACEMD installed (pip install acemd matched to your CUDA build).
The Build a protein tutorial walks through the build steps reused below.

The flow#

ACEMD’s Python API exposes two convenience setups - one per protocol stage - plus the runner itself:

acemd.protocols.setup_equilibration() - writes the equilibration input.yaml and copies the structure / parameter / coordinate files into a fresh equilibration directory.
The acemd CLI - run from the shell against the equilibration directory, reads the input file and produces trajectory + log + restart files alongside.
acemd.protocols.setup_production() - reads the equilibrated state from step 2 and writes a production input.yaml into a new production directory.
The acemd CLI again - run against the production directory.

Equilibration is short (a few ns of NPT with restraints relaxing); production is the long unrestrained NVT/NPT run that generates the trajectory you’ll analyse.

Note

This tutorial executes the setup_equilibration step so you can see the output it generates. The two acemd invocations and the production setup are shown but not executed - a real run takes hours to days on a GPU, well outside a tutorial. Copy the code blocks into a script and launch it on a workstation or queue when you’re ready to simulate for real.

Setup#

from moleculekit.molecule import Molecule
from acemd.protocols import setup_equilibration, setup_production

Step 1 - Get a built, solvated, ionised system#

For this tutorial we reuse a Trp-cage (1L2Y) system that the canonical protein build already produced - load → segment → prepare → solvate → AMBER build with ionisation. ACEMD doesn’t care which builder produced the files - an AMBER structure.prmtop + structure.pdb pair, a CHARMM PSF + coords, or an OpenMM XML force field all work as inputs. We solvate and ionise here because for explicit-solvent runs you’ll usually want both - implicit-solvent runs (gbsa=True on the build) are also supported if you want to skip the water box.

mol = Molecule("1L2Y")
mol = autoSegment(mol, fields=("segid", "chain"))
prepared, _ = systemPrepare(mol, pH=7.4)
solvated = solvate(prepared, pad=12)
amber.build(solvated, outdir="./build", ionize=True, saltconc=0.15)

To keep the docs build fast we skip the actual build here and pick up from the resulting ./build/ directory. Load it back into a Molecule to verify:

mol = Molecule("./build/structure.prmtop")
mol.read("./build/structure.pdb")

Copyright by Acellera Ltd. By executing you are accepting the License. In order to register, run htmd_register on your terminal.
The registration information must be valid so that it might be verified.

Step 2 - Set up the equilibration#

setup_equilibration("./build", "./equil", run="4ns")

2026-06-06 13:11:47,348 - acemd - INFO - Copying ./build/structure.pdb to ./equil/structure.pdb
2026-06-06 13:11:47,349 - acemd - INFO - Copying ./build/structure.prmtop to ./equil/structure.prmtop

setup_equilibration() reads the structure + parameters from ./build/, copies them into ./equil/, and writes an input.yaml that configures:

A short steepest-descent minimisation (minimize: 500 steps).
An NPT ensemble: thermostat: true at thermostattemperature: 300 K, barostat: true at the default 1 atm.
velocities: 300 K - Maxwell-Boltzmann velocities seeded at this temperature when no checkpoint exists; on a resumed run the velocity field is loaded from the checkpoint instead.
restart: true - on re-invocation of acemd --input ./equil, the run resumes from ./equil/restart.chk rather than starting over.
Default positional restraints in extforces: four entries by default - protein Cα with 1@0, protein heavy atoms (non-Cα) with 0.1@0, nucleic backbone with 1@0, nucleic non-backbone heavy with 0.1@0. Each decays linearly to zero by step 500000, gradually releasing the system into the equilibrated water box. Pass defaultrestraints=False to drop them, or add your own via extforces=[...].
The integration length you passed in (run: 4ns), and the periodic box size taken from the build (boxsize: [...]).

Inspect the produced files:

sorted(os.listdir("./equil"))

['input.yaml', 'run.sh', 'structure.pdb', 'structure.prmtop']

input.yaml is the file ACEMD reads; the other files are the structure / parameters / starting coordinates it references. The YAML itself is human-readable - inspect it to see (and override) any of the integration parameters:

print(open("./equil/input.yaml").read())

barostat: true
barostatconstratio: false
boxsize:
- 47.91499710083008
- 45.03700256347656
- 38.792999267578125
coordinates: structure.pdb
extforces:
- sel: protein and name CA
  setpoints:
  - 1@0
  - 0@500000
  type: positionalRestraint
- sel: protein and noh and not name CA
  setpoints:
  - 0.1@0
  - 0@500000
  type: positionalRestraint
minimize: 500
restart: true
run: 4ns
structure: structure.prmtop
thermostat: true
thermostattemperature: 300
velocities: 300

Step 3 - Run the equilibration#

Hand the equilibration directory to the acemd CLI from your shell. This is the slow step - it needs a GPU and runs for hours - so don’t try to execute it inside a notebook.

acemd --input ./equil

ACEMD reads ./equil/input.yaml, runs minimisation + the requested ns of NPT, and writes the trajectory (output.xtc), the final state (output.coor, output.vel, output.xsc), and the checkpoint (restart.chk) back into ./equil/. The bundled ./equil/run.sh wrapper redirects ACEMD’s stdout/stderr into log.txt; if you call acemd --input ./equil directly the log goes to your terminal instead. When the command exits, the system is ready for production.

Step 4 - Set up production#

Production reuses the last frame of the equilibration as its starting state.

setup_production("./equil", "./prod", run="100ns", temperature=300)

setup_production() reads the final coordinates and box (output.coor + output.xsc) from ./equil/, copies them into ./prod/, and writes a production input.yaml configured for an unrestrained NVT run at the given temperature for the requested length. Velocities are regenerated at Maxwell-Boltzmann at temperature rather than copied from equilibration (set barostat=True to switch to NPT if you need pressure coupling in production).

Step 5 - Run production#

acemd --input ./prod

The trajectory lives in ./prod/output.xtc and is what you’d feed into the MSM analysis tutorials (or any other downstream analysis).

Running on a queue#

For longer campaigns - many starting structures, replicas, or adaptive epochs - drive ACEMD through HTMD’s job queues instead of running the acemd CLI by hand:

from htmd.ui import SlurmQueue

queue = SlurmQueue()         # or LocalGPUQueue() / LsfQueue() / PBSQueue()
queue.submit("./prod")       # ./prod is the directory setup_production produced
queue.wait()                 # blocks until all submitted runs finish

The queues run ACEMD on each directory you submit and stream the resulting trajectories back. This is the pattern the adaptive-sampling tutorials build on.

Parameters that matter#

Argument	Effect
`run="4ns"` (equilibration)	Equilibration length.
`run="100ns"` (production)	Production length. A single MD trajectory; pick to match the timescale you’re studying or use adaptive sampling for orders-of-magnitude reach.
`temperature=300`	Production temperature in K. Match the equilibration.
`setup_equilibration(..., extforces=[...])`	List of external-force dicts (e.g. positional restraints on a ligand). See `setup_equilibration()` for the schema.
`setup_equilibration(..., minimize=300)`	Override the minimisation step count. Useful when starting from a freshly-relaxed structure.

Gotchas#

The two acemd invocations dominate wall time; everything else (the build, the two setup_* calls) is seconds. Profile your run on a short equilibration before queuing up many production trajectories.
The output trajectory format follows the extension of trajectoryfile (default output.xtc); pass e.g. setup_production(..., trajectoryfile="output.dcd") to switch. MSM analysis through Sim and Metric reads both XTC and DCD natively.