Membrane simulations with cellular restraints

You will learn: how to keep a ligand on the correct side of the membrane during a long production run, using acemd.protocols.get_cellular_restraints() to derive a set of flat-bottomed restraints that follow the periodic image.

Prerequisites:

• HTMD installed.

• ACEMD installed.

• You’ve worked through Build a membrane-embedded protein and Run an MD simulation with ACEMD.

Why this matters

In a membrane-protein simulation the box wraps periodically in all three directions. The bilayer extends across the XY plane, with the protein crossing it once - but in the +z / -z direction, a ligand or peptide on the extracellular side can diffuse out of the top of the box and reappear at the bottom, suddenly sitting on the intracellular side. The same atom is one image away; the simulation doesn’t know the difference. For binding / unbinding studies this completely scrambles the kinetics.

acemd.protocols.get_cellular_restraints() derives a set of flat-bottomed positional restraints that confine a chosen molecule to its native compartment, optionally allowing it to enter but not cross the bilayer, or to cross only one way. The restraints are anchored to the lipid centre of mass so they track the membrane through translation but assume the box height stays fixed - hence “only use these in NVT” (the docstring spells this out).

The flow

1. Start from a pre-built membrane-embedded system (see tutorial 07 for the full build).

2. Compute the restraint list with get_cellular_restraints().

3. setup_equilibration() with conventional positional restraints on the ligand (because equilibration runs as NPT and the box height changes).

4. Run the equilibration via the acemd CLI - this produces the output.coor (final equilibrated coordinates) and output.xsc (final box) that production setup reads.

5. setup_production() with the cellular restraints for the NVT production run.

6. Run production via the acemd CLI.

Note

Steps 1-3 execute in the docs build. Step 5 (setup_production) is shown but not executed because it needs the equilibration’s output.coor + output.xsc, and steps 4 / 6 (the acemd CLI) run on a GPU - copy the commands into a script when you’re ready.

Setup

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

Step 1 - Load a built apelin receptor + membrane system

The membrane tutorial walks through the OPM-aligned 5VBL build end to end. We reuse its output: structure.prmtop + structure.pdb for the apelin receptor (a GPCR) embedded in a POPC + cholesterol bilayer with the agonist peptide bound on the extracellular side. The pre-built pair lives under tutorials/simulation/_systems/; the cell below copies it into a local ./build/ and loads it into a Molecule. The box dimensions come from the PDB’s CRYST1 record - get_cellular_restraints needs them.

built = Molecule("./build/structure.prmtop")
built.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 - Compute the cellular restraints

The apelin peptide sits in segment P0 (chain A of the original PDB). On the OPM-aligned structure, +z points toward the extracellular side of the membrane, so the peptide lives in the extracellular compartment. We want a restraint that keeps it there throughout production.

membrane_rel_z is the fraction of the box height at which the bilayer is centred. The default 0.5 only holds when the bilayer sits exactly mid-box; in practice the build’s z-padding above and below the membrane is rarely symmetric, so compute it from the actual lipid coordinates:

lipid_z = built.coords[built.atomselect("lipid"), 2, 0]
membrane_rel_z = float(lipid_z.mean() / built.box[2, 0])
print(f"membrane centred at z = {lipid_z.mean():.1f} / {built.box[2, 0]:.1f} "
      f"-> membrane_rel_z = {membrane_rel_z:.3f}")

membrane centred at z = 69.0 / 125.0 -> membrane_rel_z = 0.552

extracellular_sel = 'segid "P0"'   # apelin peptide
intracellular_sel = None           # no intracellular ligand in this system

restraints = get_cellular_restraints(
    built,
    extracellular_sel=extracellular_sel,
    intracellular_sel=intracellular_sel,
    membrane_rel_z=membrane_rel_z,
    extracellular_crossing=None,   # don't even enter the membrane
)
for r in restraints:
    print(r)

{'type': 'positionalrestraint', 'sel': 'segid "P0"', 'axes': 'z', 'fbwidth': 35.85469055175781, 'fbcenter': '(lipid or resname AR CHL DHA LAL MY OL PA PC PE PGR PGS PS SA SPM ST) and noh', 'fbcenteroffset': [0, 0, 38.13734436035156], 'setpoints': ['1@0']}

Each entry is a dict in ACEMD’s extforces format - the same shape setup_equilibration and setup_production accept. The restraint is a one-dimensional flat-bottomed potential along z, centred on the lipid centre of mass with an offset that places the “free zone” entirely in the extracellular compartment. Inside the free zone the ligand feels nothing; cross the boundary and a harmonic force kicks in to push it back.

The extracellular_crossing knob controls what’s allowed:

Value	Behaviour
"cross"	Ligand may travel through the membrane and reach the intracellular side (still no PBC wrapping).
"enter"	Ligand may enter the membrane but not cross it - useful for partition-coefficient studies.
None	Ligand stays in the extracellular compartment - the default for “keep it on one side”.

intracellular_crossing is the mirror knob for any intracellular ligand. Pass None for either selection when you don’t need a restraint on that side.

Step 3 - Equilibration (NPT) with positional restraints

Equilibration runs NPT - the box height changes as the lipid bilayer relaxes - so cellular restraints aren’t safe here. Use simple positional restraints to keep the peptide near its starting coordinates while the rest of the system settles.

equil_restraints = [
    {
        "type": "positionalrestraint",
        "sel": f"{extracellular_sel} and noh",
        "setpoints": ["5@0"],
    }
]

setup_equilibration(
    "./build", "./equil",
    run="4ns",
    barostatconstratio=True,           # XY scale together, Z relaxes independently - required for membrane systems
    extforces=equil_restraints,
)
print(open("./equil/input.yaml").read())

barostat: true
barostatconstratio: true
boxsize:
- 87.0689926147461
- 90.40899658203125
- 122.00299835205078
coordinates: structure.pdb
extforces:
- sel: segid "P0" and noh
  setpoints:
  - 5@0
  type: positionalrestraint
- 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

2026-06-06 13:11:58,759 - acemd - INFO - Copying ./build/structure.pdb to ./equil/structure.pdb
2026-06-06 13:11:58,762 - acemd - INFO - Copying ./build/structure.prmtop to ./equil/structure.prmtop

The peptide’s heavy atoms (... and noh) are held with a force constant of 5 kcal/mol/Ų for the whole 4 ns - long enough for the bilayer and water box to settle without the peptide drifting. Hydrogens are left unrestrained so the peptide can still relax internally. The minimisation + protein-backbone defaults from setup_equilibration() still fire on top of this.

Step 4 - 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

When the command exits, ./equil/ contains:

• output.coor - the final equilibrated coordinates.

• output.xsc - the final equilibrated box.

• output.vel - the final velocities (production regenerates its own at the requested temperature, so this isn’t strictly needed downstream).

• output.xtc + restart.chk - the trajectory and resume checkpoint. The bundled ./equil/run.sh wrapper redirects ACEMD’s stdout/stderr into log.txt; calling acemd --input ./equil directly sends it to the terminal.

setup_production reads output.coor (starting coordinates) and output.xsc (the production box) from this directory, so you must run equilibration before setting up production. The topology (structure.prmtop) and force-field parameter files are reused from the equilibration directory as-is.

Step 5 - Production setup (NVT) with cellular restraints

Once equilibration has fixed the box dimensions, switch to NVT for production and replace the positional restraints with the cellular ones from step 2.

setup_production(
    "./equil", "./prod",
    run="100ns",
    temperature=300,
    extforces=restraints,   # the cellular restraints we computed above
)

setup_production already defaults barostat=False (production is NVT by default), so we don’t need to pass it — that’s exactly the constant-box regime the cellular restraints assume. The extforces list flows straight from get_cellular_restraints into the production input.yaml.

Step 6 - Run production

acemd --input ./prod

Gotchas

• get_cellular_restraints reads the box from mol.box, the lipid positions from mol.coords, and the membrane location from membrane_rel_z. Always compute membrane_rel_z from mol.coords[mol.atomselect("lipid"), 2, 0].mean() / mol.box[2, 0] rather than hardcoding 0.5 - the asymmetric water padding above and below the bilayer shifts the relative centre away from mid-box in practice. (mol.box is (3, n_frames), so the [2, 0] indexing picks the box-z of the first frame.)

• NVT only. The restraints are anchored to the lipid centre and to a fraction of the box height; an NPT run where the box rescales would silently misalign the free-zone with the lipids. Equilibrate the box dimensions under NPT with conventional positional restraints first, then switch to NVT for production.

• lipidsel defaults to moleculekit’s "lipid" macro plus a list of common lipid resnames (POPC, POPE, POPG, CHL1, etc.). If you’re using a custom lipid that isn’t in the macro, pass lipidsel="<your selection>" explicitly.

• Use barostatconstratio=True on setup_equilibration for membrane systems - without it ACEMD warns that the box scaling will distort the bilayer (lateral pressure should scale isotropically in XY while Z relaxes independently).

See also

• Run an MD simulation with ACEMD - the canonical equilibration → production flow this tutorial extends.

• Build a membrane-embedded protein - the long-form membrane build whose output we load here.