Set up equilibration and production#

You will learn: how to turn an HTMD build directory into a full equilibration → production pipeline using setup_equilibration() and setup_production(), instead of writing input.yaml files by hand.

Prerequisites:

ACEMD installed.
A build directory containing a CHARMM (structure.psf + *.prm) or AMBER (structure.prmtop) topology plus a starting PDB. The HTMD builder produces this directly.

What the helpers do for you#

Each function writes a sensible input.yaml plus a run.sh launcher into an output directory:

setup_equilibration() — NPT, thermostat + barostat on, light minimization, default protein/nucleic-acid restraints that decay over the first half of the run, cis-peptide bond check, auto-box detection for vacuum systems, lipid auto-detection that warns about barostatconstratio.
setup_production() — NVT, picks up the final box and coordinates from the equilibration directory, no default restraints.

Standard pipeline#

Each helper runs before its corresponding ACEMD simulation and reads files from a previous directory to seed the next one:

setup_equilibration() reads the builder directory (./build/) — topology, parameters, starting coordinates — and writes an equilibration input directory (./equil/).
setup_production() reads the equilibration directory (./equil/) — final coordinates (output.coor), final box (output.xsc), and the topology / parameters that setup_equilibration() already copied there — and writes a production input directory (./prod/). Velocities are not carried over; production starts with fresh Maxwell-Boltzmann velocities at temperature.

That means the equilibration must have finished (so output.coor and output.xsc exist in ./equil/) before setup_production() is called.

Step 1 — Generate the equilibration input#

from acemd.protocols import setup_equilibration

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

This writes ./equil/input.yaml, a ./equil/run.sh launcher, and copies the topology / coordinates / parameters from ./build/ into ./equil/.

Step 2 — Run the equilibration#

From a shell:

cd equil && ./run.sh

Or from Python:

from acemd import acemd
acemd("./equil")

The run writes output.coor, output.vel, output.xsc, output.xtc, and restart.chk into ./equil/.

Step 3 — Generate the production input#

Once ./equil/output.* exists:

from acemd.protocols import setup_production

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

This writes ./prod/input.yaml and ./prod/run.sh, with the final coordinates and box copied from ./equil/output.coor and ./equil/output.xsc. New Maxwell-Boltzmann velocities at temperature are generated at the start of the production run.

Step 4 — Run the production#

cd prod && ./run.sh

Or:

acemd("./prod")

Defaults you should know about#

setup_equilibration() writes an input.yaml with:

thermostat: true at temperature (default 300 K — not ACEMD’s 298.15 K default).
Initial velocities: temperature — Maxwell-Boltzmann at that same value.
barostat: true (NPT).
minimize: 500.
restart: true, so re-running picks up from a checkpoint if one exists.
Positional restraints on:
- Protein Cα atoms — 1 kcal/mol/Å²
- Other protein heavy atoms — 0.1 kcal/mol/Å²
- Nucleic-acid backbone — 1 kcal/mol/Å²
- Other nucleic-acid heavy atoms — 0.1 kcal/mol/Å²
Each restraint decays linearly to 0 over the first half of the run by default; override with restraintdecay="5ns" or any time / step count.
boxsize is auto-derived from the build’s coordinates if you don’t pass one — the bounding box of the waters (plus 1 Å padding) for solvated systems, or of all atoms (plus 12 Å) for vacuum systems.

Turn the defaults off for systems that don’t need them:

setup_equilibration("./build", "./equil", run="10ns", defaultrestraints=False)

setup_production() writes an input.yaml with:

thermostat: true at temperature, barostat: false (NVT), restart: true.
No default restraints.
coordinates copied from ./equil/output.coor, boxsize from ./equil/output.xsc, plus the topology and parameters (whichever ones setup_equilibration() had copied into ./equil/).
Initial velocities: temperature — fresh Maxwell-Boltzmann, not read from ./equil/output.vel.

Override any input option#

Both functions accept any ACEMD input-file option as a keyword argument — it lands straight in the generated input.yaml:

setup_equilibration(
    "./build", "./equil",
    run="10ns",
    temperature=310,
    timestep=2,
)

Add extra external forces on top of the defaults:

setup_equilibration(
    "./build", "./equil",
    run="10ns",
    extforces=[{
        "type": "positionalRestraint",
        "sel": "resname LIG",
        "setpoints": ["1@0"],
    }],
)

Membrane systems#

When the builder output contains lipids (detected via the standard lipid residue names), setup_equilibration() warns if barostatconstratio is unset. Turn it on so the membrane’s xy aspect ratio is preserved during NPT:

setup_equilibration("./build", "./equil", run="10ns", barostatconstratio=True)

For membrane-ligand work — keeping a ligand on a chosen side of the bilayer instead of letting it drift through the periodic image — get_cellular_restraints() builds the right extforces block. The function assumes the membrane lies in the xy plane with +z pointing extracellular, and uses flat-bottom potentials that only kick in when the ligand exits the allowed slab.

Important

get_cellular_restraints() is for NVT production runs only — its docstring is explicit: “Only use these restraints in NVT simulations as they require the box size to remain constant.” During NPT equilibration the box changes shape and the absolute z coordinates of the restraint slab would drift. For equilibration, lean on the default positional restraints (which pin the ligand in place anyway) or pass an explicit positionalRestraint if you need something tighter.

Worked pipeline for a GPCR ligand that should stay on the extracellular side:

from acemd.protocols import (
    setup_equilibration,
    setup_production,
    get_cellular_restraints,
)
from moleculekit.molecule import Molecule
from acemd import acemd

# 1. Equilibrate under NPT. The default protein restraints keep the system
#    in place; barostatconstratio preserves the membrane's xy aspect ratio.
setup_equilibration(
    "./build", "./equil",
    run="20ns",
    barostatconstratio=True,
)
acemd("./equil")

# 2. Build the cellular restraints from the equilibrated system. Read the
#    structure plus the final coordinates / box so the restraint geometry
#    is computed against the post-equilibration state.
mol = Molecule("./equil/structure.prmtop")
mol.read("./equil/output.coor")
mol.read("./equil/output.xsc")

extforces = get_cellular_restraints(
    mol,
    extracellular_sel="resname LIG",
    intracellular_sel=None,           # no intracellular species to restrain
    membrane_rel_z=0.5,                # membrane sits in the middle of the box
    extracellular_crossing=None,       # don't even let it enter the membrane
)

# 3. Production under NVT with the cellular restraints applied.
setup_production(
    "./equil", "./prod",
    run="500ns",
    extforces=extforces,
)
acemd("./prod")

extracellular_crossing (and intracellular_crossing) accept None (block at the membrane surface), "enter" (allow entering the bilayer but not crossing), or "cross" (no restraint applied). See the function’s docstring for the full signature including the lipidsel override.

Gotchas#

builddir must hold the topology and starting coordinates. Auto-detected names (in priority order):
- coordinates: structure.pdb, then any *.pdb.
- structure: structure.prmtop, structure.psf, then any *.prmtop or *.psf.
- parameters: a file literally named parameters, then any *.prm.
Pass structure=, coordinates=, parameters= explicitly if yours don’t match.
The default temperature is 300 K, different from ACEMD’s thermostat default of 298.15 K. Pass temperature=298.15 to match.
Cis-peptide bonds emit a warning, not an error. A real cis peptide is rare; a flagged cis peptide is almost always a build mistake. Inspect before proceeding.
run.sh assumes a conda environment. The generated launcher does conda activate $CONDA_DEFAULT_ENV and runs acemd >log.txt 2>&1. If you don’t use conda, just call acemd directly inside the directory.