# Set up equilibration and production

**You will learn:** how to turn an HTMD build directory into a full equilibration → production pipeline using {py:func}`~acemd.protocols.setup_equilibration` and {py:func}`~acemd.protocols.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](https://software.acellera.com/htmd/) 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:

- {py:func}`~acemd.protocols.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`.
- {py:func}`~acemd.protocols.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:

- {py:func}`~acemd.protocols.setup_equilibration` reads the **builder** directory (`./build/`) — topology, parameters, starting coordinates — and writes an equilibration input directory (`./equil/`).
- {py:func}`~acemd.protocols.setup_production` reads the **equilibration directory** (`./equil/`) — final coordinates (`output.coor`), final box (`output.xsc`), and the topology / parameters that {py:func}`~acemd.protocols.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 {py:func}`~acemd.protocols.setup_production` is called.

### Step 1 — Generate the equilibration input

```python
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:

```bash
cd equil && ./run.sh
```

Or from Python:

```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:

```python
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

```bash
cd prod && ./run.sh
```

Or:

```python
acemd("./prod")
```

## Defaults you should know about

{py:func}`~acemd.protocols.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:

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

{py:func}`~acemd.protocols.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 {py:func}`~acemd.protocols.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`:

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

Add extra external forces on top of the defaults:

```python
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), {py:func}`~acemd.protocols.setup_equilibration` warns if `barostatconstratio` is unset. Turn it on so the membrane's xy aspect ratio is preserved during NPT:

```python
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 — {py:func}`~acemd.protocols.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}
{py:func}`~acemd.protocols.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:

```python
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.

## See also

- [Use positional restraints](use-positional-restraints.md) — what the auto-applied restraints actually do.
- [Input options](../reference/input-options.md).
- {py:func}`~acemd.protocols.setup_equilibration`, {py:func}`~acemd.protocols.setup_production`, {py:func}`~acemd.protocols.get_cellular_restraints`.
- [HTMD builder](https://software.acellera.com/htmd/) — for producing the `builddir` input.