Metadata-Version: 2.5
Name: caddforge-md
Version: 0.13.0
Summary: Molecular dynamics workflows for the CADDForge toolbox.
Author-email: Manuel Sellner <manuel.sellner@idorsia.com>, Christoph Sager <christoph.sager@idorsia.com>
License: MIT
License-File: LICENSE
Requires-Python: >=3.10
Requires-Dist: caddforge-core>=0.21.0
Requires-Dist: caddforge-descriptors>=0.8.0
Requires-Dist: caddforge-docking>=0.3.0
Requires-Dist: matplotlib>=3.10.8
Requires-Dist: mdanalysis>=2.9.0
Requires-Dist: mdtraj>=1.10.3
Requires-Dist: parmed>=4.3.0
Requires-Dist: rdkit>=2022.9.5
Requires-Dist: scikit-learn>=1.1.0
Provides-Extra: gui
Requires-Dist: caddforge-gui>=0.2.0; extra == 'gui'
Requires-Dist: fastapi>=0.110; extra == 'gui'
Requires-Dist: python-multipart>=0.0.21; extra == 'gui'
Provides-Extra: minimize
Requires-Dist: numpy>=1.26; extra == 'minimize'
Requires-Dist: openmm>=8.2; extra == 'minimize'
Requires-Dist: scipy>=1.11; extra == 'minimize'
Provides-Extra: test
Requires-Dist: hypothesis>=6.0; extra == 'test'
Requires-Dist: pytest>=8.0; extra == 'test'
Requires-Dist: tomli>=2; (python_version < '3.11') and extra == 'test'
Description-Content-Type: text/markdown

# CADDForge MD Plugin

Install this plugin with `python -m pip install caddforge-md`. It depends directly
on `caddforge-core>=0.21.0`; Python imports and CLI commands still use `caddforge`.
For the complete public toolbox, including user-facing extras, install
`python -m pip install caddforge`. When migrating an existing environment, uninstall
the former `caddforge` core distribution before installing renamed core and updated plugins.

`caddforge-md` provides molecular dynamics workflows as plugins for CADDForge.

## Implemented Tools
- [`fix-structure`](#tool-fix-structure)
- [`center-and-align`](#tool-center-and-align)
- [`extract-frames`](#tool-extract-frames)
- [`cluster-representative`](#tool-cluster-representative)
- [`plain-md`](#tool-plain-md)
- [`openbpmd`](#tool-openbpmd)
- [`grand`](#tool-grand)
- [`inducia`](#tool-inducia)
- [`minimize`](#tool-minimize)

## Tool: `minimize`

Minimize one protein–ligand complex with OpenMM L-BFGS. This is local pose
cleanup: the reported physical-model energy includes implicit solvation and is
**not a binding affinity or a score for comparing different chemical systems**.
The tool emits finite poses with quality flags by default; strict mode can reject
flagged jobs while retaining their structures and diagnostics.

### Installation and execution

Use an environment containing AmberTools (`pdb4amber`, `antechamber`, `parmchk2`,
`tleap`), OpenMM >=8.2, ParmEd, RDKit, NumPy, SciPy, and this plugin. For example:

```bash
conda create -n caddforge-minimize -c conda-forge python=3.13 ambertools openmm=8.2 parmed rdkit numpy scipy
conda activate caddforge-minimize
python -m pip install -e '.[minimize]'
caddforge show-config md minimize
caddforge run md minimize --config examples/minimize/ligand_only.yaml
```

AmberTools is installed separately from the Python extra. A prepared Amber input
does not run AmberTools. Scientific imports are deferred until execution.
CLI and pipeline execution use core input staging, including transparent gzip
for declared input files. Checksums in chemistry/mapping records refer to the
decompressed bytes. Nested cofactor chemistry paths are ordinary uncompressed
configuration inputs and must be available on the execution node.

```python
from cadd_md.minimize import MinimizeConfig, MinimizeWorkflow

config = MinimizeConfig(
    complex_pdb_path="complex.pdb",
    ligand_path="ligand.sdf",
    lig_resname="LIG",
    mode="pocket_restrained",
    run_directory="runs/cleanup",
    registry_path="runs/job_registry.sqlite",
)
workflow = MinimizeWorkflow(config)
result = workflow.run(workflow.prepare_job())
print(result.status, result.metrics)
print(result.paths["minimized_complex"])
```

Paths follow the existing config convention: relative paths resolve from the
submission working directory. Use the core CLI/backend for Slurm submission;
the scientific workflow contains no scheduler-specific execution logic.
Choose a fresh `run_directory` for each run; existing minimization outputs are
preserved and cannot be overwritten by a new job.

### Inputs and physical model

Supply either `complex_pdb_path`, or `protein_pdb_path` plus a ligand pose in the
same coordinate frame. `ligand_path` is always required and supplies authoritative
SDF/MOL/MOL2 chemistry. For a complex PDB, coordinates come from that complex;
template coordinates are not used to align or move the pose. For separate inputs,
the ligand file must contain one explicitly hydrogenated 3D molecule. The target
residue is selected by `lig_resname`, which must occur exactly once.

Preparation preserves supplied atoms, protonation, ligand stereochemistry, and
disulfides, adds only missing hydrogens, and rejects missing protein heavy atoms.
Repair incomplete proteins upstream with `md/fix-structure`. Bare histidines need
unambiguous hydrogen evidence or an explicit `residue_states` entry. Optional
`protein_states_path` and `protein_links_path` must be provided together and match
the input structure checksum. Covalent target ligands are rejected.

Automatic preparation uses ff14SB, GAFF2/AM1-BCC, TIP3P retained waters, and OBC2
with `mbondi2` radii, ACE surface area, solute dielectric 1, solvent dielectric
78.5, and zero salt screening. It is nonperiodic with `NoCutoff`; no solvent box
or neutralizing ions are added. Retained waters are explicit solute particles in
the GB model, a hybrid approximation requiring application-specific validation.
All force-field choices, parameterization commands, checksums, and available
software versions are retained.
Both atom typing and missing-parameter lookup explicitly select GAFF2. Unresolved
`parmchk2` parameter warnings fail preparation and require reviewed parameters.

Generated GAFF2/AM1-BCC charges are checked and, when necessary, normalized
independently for each organic component. Minimize and OpenBPMD apply a uniform
per-atom offset only when the total correction is at most **0.005 e** and each
atom changes by at most **0.0001 e**. These are conservative engineering limits,
not universal scientific guarantees. Final component charges must sum to the
authoritative formal charge within **1e-8 e**, and mapped charges are verified
again after Amber topology construction. Invalid chemistry or larger residuals
fail preparation; supplied Amber charges are never normalized.

Preparation retains the original `*.raw.mol2`, the final MOL2, and a
`*.charge_validation.json` report with raw/final sums, correction, policy limits,
checksums, and the decision. `charge_validation` in the preparation summary links
the reports, including rejected corrections. Applied corrections also produce a
log warning. A residual within the charge-token rounding bound is distinguished
from a larger generation residual; the latter is not assumed to be rounding.
Changing `ligand_net_charge` or disabling quality checks is not a remedy for this
error. These rules also apply to pipeline and scatter jobs without new YAML
options.

Automatic preparation also retains Antechamber's `ANTECHAMBER_AC.AC` bond
evidence. The shared validator checks mapped atoms, coordinates, atom types,
connectivity (including hydrogen attachments), and total formal charge against
the authoritative ligand. A missing or inconsistent AC file fails preparation.

Amber MOL2 bond tokens are not always literal chemical bond orders: for example,
carboxylate C–O bonds can both be written as single bonds. Validation accepts
supported Amber encodings only in their permitted chemical environments and
requires one complete authoritative or RDKit resonance form to explain every
bond. Aromatic/Kekulé alternatives are supported; arbitrary bond-order changes,
tautomers, and changed protonation are not accepted. At most 1,000 candidate
forms are examined. An inconclusive capped search fails explicitly.

The charge report contains a `graph_validation` record with source checksums,
atom mapping, candidate count, interpreted bonds, and the acceptance or rejection
reason. The validators never rewrite bond records, atom types, coordinates, or
the authoritative ligand. Charge normalization runs only after graph validation.
These rules apply equally to minimize, OpenBPMD, organic cofactors, and pipeline
jobs; no configuration override is needed for supported representations.

Noncovalent organic cofactors require `cofactors` entries with `residue_index`
and `chemistry_path`; optional `atom_map_path` uses the shared explicit atom-map
schema against the original complex. Residue and atom indices are zero-based
indices in `preparation/inputs/input_identity.json`, not PDB serial numbers.
Automatic isolated ions are restricted to Na+, K+, and Cl- and require explicit
`isolated_ion_residue_indices`. These declarations must not be used to stand in
for a coordinated metal-site model.

Coordinated metals, metall cofactors, and other supported supplied components
use all three fields: `amber_topology_path`, `amber_coordinates_path`, and
`amber_mapping_path`, alongside the complete input complex. The topology must
be nonperiodic, contain all bonded/nonbonded parameters and positive GB
radii/screening factors, and have finite charges and positive masses. Unsupported
force forms and virtual sites fail explicitly. The supplied radii and charges
are used without substituting automatic values. Every non-isolated metal needs
a `metal_sites` entry containing `atom_index`, `oxidation_state`,
`coordinating_atom_indices`, and a descriptive `parameter_source`.

The prepared Amber mapping is JSON:

```json
{
  "schema_version": 1,
  "input_sha256": "SHA256 of the complex PDB",
  "topology_sha256": "SHA256 of the Amber topology",
  "coordinates_sha256": "SHA256 of the Amber coordinates",
  "input_atom_indices": [0, 1, 2, null],
  "input_residue_indices": [0]
}
```

The arrays follow Amber topology order and must cover every input atom/residue
exactly once. `null` is allowed only for an added hydrogen. Array lengths must
match the actual topology; the short arrays above illustrate the schema only.
Parameter sources and oxidation states are declarations, not an automatic
validation of a metal model's physical accuracy.

### Regions and modes

Pocket residues have any heavy atom within `pocket_radius_angstrom: 5.0` of any
ligand heavy atom. The shell contains other protein residues with any heavy atom
within `shell_radius_angstrom: 5.0` of **any pocket heavy atom**. Both selections
expand to whole residues and stay fixed during minimization. An explicit
`pocket_residues: ["A:42", "A:43:B"]` replaces automatic pocket selection;
`_` denotes a blank chain. Unknown, duplicate, ambiguous, or empty selections fail.

| `mode` | Pocket side chains | Pocket backbone | Shell |
|---|---|---|---|
| `ligand_only` (default) | Fixed | Fixed | Fixed |
| `sidechains_free_backbone_fixed` | Free | Fixed | Fixed |
| `sidechains_free_backbone_restrained` | Free | Restrained | Fixed |
| `pocket_restrained` | Free | Restrained | Strongly restrained |
| `pocket_free_shell_restrained` | Free | Free | Strongly restrained |

The target ligand is free and distant protein is fixed in every mode. Backbone
atoms are N/CA/C/O/OXT; terminal caps count as backbone. Hydrogens inherit their
bonded parent's region. Waters, ions, metals, and cofactors are fixed by default;
`component_mobility` permits explicit whole-residue `fixed`, `free`, or
`restrained` overrides keyed by input residue index.

Positional restraints use **E = ½ k |r − r₀|²**, with **k in kJ mol⁻¹ nm⁻²**.
Named strengths are `weak: 100`, `moderate: 1000`, and `very_strong: 10000`;
positive numeric strengths are accepted. `backbone_restraint_strength` defaults
to `moderate`, and `shell_restraint_strength` to `very_strong`. They are starting
regimes, not universally validated settings. Only heavy atoms are restrained
unless `restrain_hydrogens: true`. Fixed regions include their hydrogens.

Fixed particles have zero mass in a private minimization-system copy. The
exported physical system retains original masses. Bonds and water geometry are
flexible by default. `constraints: mobile_hbonds` constrains bonds involving
hydrogen only when both endpoints are mobile, using their force-field equilibrium
lengths and retaining the physical bonded terms. Constraints involving fixed
atoms are rejected. Disulfide connectivity is preserved across region boundaries.

### Convergence, diagnostics, and quality

Defaults are `tolerance_kj_mol_nm: 10.0`, `max_iterations: 1000`,
`constraint_tolerance: 0.00001`, `platform: CPU`, and `cpu_threads: 1`.
GPU platforms accept `precision: single|mixed|double` where supported. Explicitly
requested unavailable platforms or precision settings fail instead of falling back.

The iteration budget spans OpenMM constraint restarts. The trace distinguishes
the physical-plus-positional objective from OpenMM's temporary bond-constraint
penalty. Final convergence is independently checked using RMS force components
over mobile atoms, projected onto the allowed constraint tangent space when
needed, plus maximum relative constraint error. A normal OpenMM return alone
does not establish convergence. Status is `verified_converged`, `iteration_limit`,
`not_verified`, or `failed`; the record also reports a reached budget even if the
last iterate independently satisfies the convergence tests.
The implementation targets the [OpenMM 8.2 minimizer and reporter behavior](https://github.com/openmm/openmm/blob/8.2.0/openmmapi/src/LocalEnergyMinimizer.cpp).

Diagnostics include initial/final molecular-mechanics, GB, physical, positional,
and total objective energies; physical/total RMS and maximum forces for all,
mobile, and fixed atoms; iteration counts; constraint errors; heavy-atom ligand,
pocket, backbone, shell, receptor, and component RMSDs without alignment; and
maximum displacements with atom indices. Preparation displacement is separate.
Ordinary nonbonded clashes, 1–4 overlaps, and declared metal coordination are
reported separately. Bonded and angle-connected contacts are excluded.

Default review triggers under `quality_thresholds` are:

| Field | Default |
|---|---|
| `ligand_rmsd_angstrom` | 1.0 |
| `ligand_max_displacement_angstrom` | 2.0 |
| `pocket_backbone_rmsd_angstrom` | 0.5 |
| `receptor_max_displacement_angstrom` | 2.0 |
| `absolute_energy_per_atom_kj_mol` | 10000 |
| `maximum_force_kj_mol_nm` | 100000 |
| `allowed_vdw_overlap_angstrom` | 0.6 |
| `coordination_distance_change_angstrom` | 0.3 |
| `objective_increase_absolute_kj_mol` | 0.01 |
| `objective_increase_relative` | 1e-8 |

These are review triggers, not universal validity bounds; geometric overlaps,
especially 1–4 contacts, need chemical interpretation. Nonfinite states, invalid
chemistry, or any fixed-atom movement always fail. Other findings produce
machine-readable `quality_flags`. Set `strict_quality: true` to fail on every
flag, or provide `strict_quality_flags` with selected codes such as
`convergence_unverified`, `large_ligand_displacement`, or
`residual_nonbonded_clashes`. Finite strict-failure poses remain available.

### Artifacts and pipelines

Canonical results are saved in `results/result_record.json`. Output ports expose
`complex`, `protein`, `ligand`, `diagnostics`, `regions`, `atom_mapping`,
`output_mapping`, `preparation_summary`, `physical_system`, `prepared_positions`,
`minimized_positions`, `amber_topology`, `amber_coordinates`, `structure_manifest`,
`protein_states`, `protein_links`, `receptor_protein_states`,
`receptor_protein_links`, optional `iteration_trace`, `result_record`, and `output_dir`.
The complex handoff records match the complex checksum; `receptor_*` records
match the protein output. A covalently connected environment component is retained
in the receptor export to avoid cutting chemical bonds, and this is recorded.

The ligand SDF preserves chemistry, title, and user properties. PDB coordinates
are rounded to 0.001 Å and SDF coordinates to 0.0001 Å; NumPy coordinate artifacts
retain full precision in **nanometers**, in physical-system particle order. The
Amber coordinate artifact is the prepared starting configuration. Use the
separate minimized coordinate artifact for the final pose. The exported physical
system contains no zero-mass freezes or positional restraints and still needs an
appropriate MD protocol before dynamics.

See `examples/minimize/` for all five modes, component overrides, a supplied
metal-site configuration, and `pipeline_scatter.yaml`. The pipeline uses existing
scatter/collect/filter stages and preserves per-pose provenance. Filtering quality
flags does not establish pose accuracy; do not rank unrelated ligands by total
minimization energy.

Validation: `python -m pytest tests/test_minimize_*.py`. Scientific fixtures use
generated peptides/ligands and are marked `external`; optional GPU comparison
tests skip with an explicit reason when no usable GPU platform exists.

## Tool: `fix-structure`

### What it does
The `fix-structure` tool conservatively repairs protein hydrogens in a protein-only PDB, protein-ligand complex PDB, or existing POPC/POPE membrane system before OpenMM simulations.

It uses the following repair logic:
- preserves declared protein states and rejects contradictory or unresolved states
- records normalization and rejects modified-residue substitutions in side-chain preserve mode
- applies the explicit terminal policy independently of side-chain protonation
- adds missing protein/cap hydrogens with OpenMM `Modeller.addHydrogens()`
- optionally caps neutral uncapped termini with ACE/NME
- canonicalizes NME `CH3`/`HH31`–`HH33` atom names to OpenMM `C`/`H1`–`H3` before repair, recording renames in `.atom_aliases.json`; ACE names are unchanged
- removes `CONECT` references to atoms deleted by normalization, preserving surviving bonds and reserving numbered `TER` serials when reconstructing OXT
- automatically separates mixed Maestro-style PDB residue ids and reassembles split residue records before OpenMM parsing
- repairs POPC/POPE membrane lipid hydrogens with the supplied OpenMM force fields
- validates OpenMM system creation automatically when configured membrane lipids are present
- leaves ligands, cofactors, and other non-protein residues untouched
- preserves ligand/cofactor `CONECT` records and PDB formal-charge fields by default
- runs in node-local scratch during allocated Slurm execution and copies results back afterwards

### External dependencies
This workflow requires OpenMM in the runtime environment.

### CLI usage
```bash
caddforge run md fix-structure --config fix-structure.yaml
```

### Python API usage
```python
from caddforge.core.execution import load_result_for_job, submit_workflow
from caddforge.plugins import get_tool

spec = get_tool("md", "fix-structure")
cfg = spec.config_cls(
    input_pdb="complex.pdb",
    forcefields=["amber14-all.xml", "amber14/tip3p.xml"],
    ph=7.0,
    backend="local",
)

workflow = spec.workflow_cls(cfg)
job = submit_workflow(workflow, cfg, backend_name=cfg.backend)
result = load_result_for_job(job)
if result is not None:
    print(result.paths["output_pdb"])
```

### Example config
```yaml
backend: local
input_pdb: complex.pdb
forcefields:
  - amber14-all.xml
  - amber14/tip3p.xml
ph: 7.0
protonation_policy: preserve
residue_states: {}
termini_policy: charged
strict_all: false
membrane_lipid_resnames:
  - POPC
  - POPE
validate_openmm_system: auto
sanitize_mixed_residues: auto
fix_nonstandard_residues: auto
preserve_conect: true
quiet: false
output_dir_name: output
output_filename: null
```

### Tool-specific config options

#### Required inputs
- `input_pdb` (`Path`): Input protein or protein-ligand complex in PDB format.
- `forcefields` (`list[str]`): OpenMM force-field XML names/files. This is required and intentionally has no implicit default.

#### Repair options
- `protonation_policy` (`preserve` or `assign`, default `preserve`): Preserve named states and diagnostic hydrogens. Bare ASP/GLU/LYS/CYS names mean their standard states; missing hydrogens do not imply neutral lysine or thiolate. Bare HIS needs tautomer hydrogens or an explicit map. `assign` explicitly allows OpenMM to resolve otherwise unspecified states at `ph`; deliberate states remain pinned. A pH value alone does not opt into reassignment.
- `residue_states` (`dict[str, str]`, default `{}`): Explicit Amber states keyed by `chain:number[:insertion]`, for example `{"A:42": "HIE", "B:10:A": "ASH"}`. Use `_` for a blank chain. Keys refer to the original input; ambiguous keys fail. Supported states depend on the selected force field and include ASH/GLH, HID/HIE/HIP, LYS/LYN, and CYS/CYM/CYX. CYX requires an established disulfide partner. Acid proton aliases also swap the associated oxygen labels to preserve attachment geometry.
- `ph` (`float`, default `7.0`): pH used by OpenMM when it must choose unspecified variants.
- `termini_policy` (`charged`, `preserve`, or `cap`; default `charged`): Explicit terminal treatment, described below. Overrides `protonation_policy` only at terminal groups.
- `cap_neutral_termini` (`bool`, default `false`): Deprecated alias; `true` selects `cap` unless an explicitly conflicting terminal policy is supplied.
- `strict_all` (`bool`, default `false`): Fail if non-protein residues are unmatched by the supplied force fields.
- `membrane_lipid_resnames` (`list[str]`, default `["POPC", "POPE"]`): Residue names to treat as supported membrane lipids. These lipids are repaired with OpenMM templates instead of being treated as untouched ligands/cofactors.
- `validate_openmm_system` (`auto`, `always`, or `never`; default `auto`): Run an OpenMM `ForceField.createSystem()` validation after repair. `auto` validates only when configured membrane lipids are present, preserving the existing non-membrane behavior.
- `sanitize_mixed_residues` (`auto`, `always`, or `never`; default `auto`): Rewrite mixed Maestro-style PDB residue ids and reassemble split atom records before OpenMM parsing. `auto` only rewrites when one chain/residue id contains multiple residue names, such as POPC records interrupted by CL/T3P records using the same residue number.
- `fix_nonstandard_residues` (`auto`, `always`, or `never`; default `auto`): Map curated modified amino acids to standard parents, including chemical substitutions such as MSE→MET and SEP→SER. Such substitutions are forbidden under `protonation_policy: preserve`. Terminal treatment is independent of this setting. D-amino acids and ring-rearranged residues such as PCA are not automatically converted. Changes are recorded in `repair_metadata.nonstandard_residue_normalization`.
- `preserve_conect` (`bool`, default `true`): Preserve ligand/cofactor/untouched-residue `CONECT` records.
- `quiet` (`bool`, default `false`): Suppress progress and warning messages.

Under `protonation_policy: preserve`, modified-residue substitutions fail with a
report. Terminal changes follow `termini_policy` instead.
The repair records explicit CONECT/SSBOND cross-links and detects missing
disulfides from unique SG pairs at 1.8–2.3 Å. Ambiguous partners and disulfide
thiol hydrogens fail; recorded bonds may stretch within the 1.7–2.5 Å validation
range without being reassigned as free thiols. These limits check connectivity
integrity and are not simulation restraints.

For POPC/POPE membrane systems, include force fields that cover the protein, lipids, water, and ions. For Amber-style OpenMM setups, `amber14-all.xml` includes the Amber lipid templates and `amber14/tip3p.xml` covers TIP3P water and common ions. The tool detects both standard and 4-character PDB residue-name layouts for POPC/POPE. If validation is enabled, unmatched ligands/cofactors that are intentionally left untouched are reported in `repair_metadata.openmm_system_validation.ignored_unmatched_residues` and excluded from the repair validation subset; they still need parameters in the downstream MD workflow.

### Terminal chemistry

`termini_policy` controls terminal chemistry separately from side-chain protonation:

| Policy | Behavior |
| --- | --- |
| `charged` (default) | Prepare free N-termini as positively charged and free C-termini as negatively charged, independently of pH. Repair missing terminal OXT when the backbone permits reconstruction. |
| `preserve` | Preserve the input terminal state and validate against the supplied force fields. Ambiguous or incomplete terminal groups fail with a residue-specific explanation. |
| `cap` | Add ACE/NME caps to eligible neutral or incomplete termini. Already charged or capped termini remain unchanged. |

A protonated, uncapped C-terminus (`–COOH`) is neutral. A deprotonated C-terminus
(`–COO⁻`) is charged. A neutral NME cap (`–CONHCH₃`) is chemically different from
an uncapped neutral acid. Capping can affect both N- and C-termini.

Terminal policy takes precedence over `protonation_policy` **only at terminal
groups**. Thus the default `charged` policy can remove HXT even when
`protonation_policy: preserve`; unrelated residue substitutions remain forbidden.
Missing terminal hydrogens can make preservation ambiguous (for example, only
one H on an ordinary free N-terminus). Prepare the intended terminal state
explicitly rather than relying on pH to select it in preserve mode. Hydrogen
repair still fills unrelated missing hydrogens.

Stock OpenMM `amber14-all.xml` supports standard charged termini and ACE/NME caps,
but does not include a neutral uncapped `–COOH` template. To preserve that form,
supply a compatible force field or additional, scientifically validated templates:

```yaml
input_pdb: data/protein_with_explicit_terminal_states.pdb
forcefields:
  - amber14-all.xml
  - /absolute/path/to/your_validated_terminal_templates.xml
protonation_policy: preserve
termini_policy: preserve
```

The extra XML above is user-supplied, not bundled with CADDForge. Force-field
compatibility is checked; preserve mode never silently converts an unsupported
neutral terminus. Setting `validate_openmm_system: never` does not disable
protein/cap template matching or terminal-state validation. This feature concerns
terminal states of recognized protein residues, not arbitrary modified residues.

For standard charged preparation, set `termini_policy: charged`. For explicitly
requested caps, set `termini_policy: cap`; side-chain preservation can stay enabled.
Approximate reconstructed atoms and caps should be minimized before production MD.

**Migration:** `cap_neutral_termini: true` is a deprecated alias for
`termini_policy: cap`. Combining it with an explicit different policy fails.
`false` does not override the new policy. `fix_nonstandard_residues` now controls
modified-residue normalization only; setting it to `never` does not disable
terminal treatment. Users who previously relied on `protonation_policy: preserve`
to reject terminal conversions must now also set `termini_policy: preserve`.

The report at `paths.normalization` records requested/effective policy, original
terminal evidence, normalization edits or added caps, final decisions, and
validation status under `normalization.terminal_treatment`. It is retained on
failure. Successful results also expose `repair_metadata.terminal_treatment`.

### Pipeline example

```yaml
version: 1
name: prepare-protein
stages:
  - id: repair
    tool: md/fix-structure
    config:
      input_pdb: data/protein.pdb
      forcefields: [amber14-all.xml, amber14/tip3p.xml]
      protonation_policy: preserve
      termini_policy: charged
    outputs:
      protein: output_pdb
  - id: cleanup
    tool: md/minimize
    inputs:
      protein: repair.protein
    config:
      ligand_path: data/ligand.sdf
      mode: ligand_only
```

The input must be a protein-only receptor for this example; ligand chemistry is
supplied separately to `minimize`. The repair tool's existing input port
`structure` and output port `repaired_structure` are unchanged. The explicit
`protein: output_pdb` output alias above binds the same PDB for the next stage.

### Outputs
For each run, artifacts are written under `runs/<date>/<job-id>/`:
- `results/<output_dir_name>/<input_stem>_fixed.pdb` or `results/<output_dir_name>/<output_filename>`: repaired PDB
- `results/results.json`: workflow summary
- `results/result_record.json`: standardized CADDForge `Result`
- `stdout.log` / `stderr.log`: runtime logs

`Result.paths` also exposes `protein_states`, `protein_links`, `normalization`,
and `atom_aliases` JSON diagnostics when produced. The state report includes each
decision's source, input/output hydrogen inventory, and serialized-output hash.
Small diagnostics are retained separately under `results/chemistry_diagnostics/`,
including on repair or bulk-copy failure.

Pipeline artifact names:
- input: `structure`
- output: `repaired_structure`

## Tool: `center-and-align`

### What it does
The `center-and-align` tool centers, wraps, optionally aligns, and writes a processed MD trajectory through the standard CADDForge workflow/job interface.

It supports:
- topology and trajectory inputs readable by MDAnalysis
- trajectory outputs writable by MDAnalysis, inferred from `output_trajectory_filename`
- frame-0 structure outputs writable by MDAnalysis, inferred from `output_structure_filename`
- configurable center, keep, and align atom selections
- optional unwrapping, centering, wrapping, and rotational/translational fitting
- standardized CADDForge outputs (`results.json`, `result_record.json`)

### External dependencies
This workflow requires `MDAnalysis`, which is already a dependency of `caddforge-md`.

### CLI usage
```bash
caddforge run md center-and-align --config center-and-align.yaml
```

### Python API usage
```python
from caddforge.core.execution import load_result_for_job, submit_workflow
from caddforge.plugins import get_tool

spec = get_tool("md", "center-and-align")
cfg = spec.config_cls(
    topology_path="system.pdb",
    trajectory_path="traj.xtc",
    center_selection="protein or resname MOL",
    keep_selection="protein or resname MOL",
    align_selection="protein and backbone",
    output_trajectory_filename="traj-centered.xtc",
    output_structure_filename="traj-centered.pdb",
    backend="local",
)

workflow = spec.workflow_cls(cfg)
job = submit_workflow(workflow, cfg, backend_name=cfg.backend)
result = load_result_for_job(job)
if result is not None:
    print(result.status.value, result.metrics)
```

### Example config
```yaml
backend: local
topology_path: system.pdb
trajectory_path: traj.xtc
center_selection: protein or resname MOL
keep_selection: protein or resname MOL
align_selection: protein and backbone
unwrap: true
center: true
wrap: true
center_mode: geometry
wrap_compound: residues
align_weights: mass
output_dir_name: output
output_trajectory_filename: traj-centered.xtc
output_structure_filename: traj-centered.pdb
```

### Tool-specific config options

#### Required inputs
- `topology_path` (`Path`): Topology/structure file readable by MDAnalysis.
- `trajectory_path` (`Path`): Trajectory file readable by MDAnalysis.

#### Selection options
- `center_selection` (`str`, default `"all"`): Atom selection used for unwrapping, centering, and wrapping.
- `keep_selection` (`str`, default `"all"`): Atom selection written to output files.
- `align_selection` (`str | null`, default `null`): Optional atom selection aligned to frame 0 after centering. Set it to `protein and backbone` for the previous protein-ligand behavior.

#### Transformation options
- `unwrap` (`bool`, default `true`): Apply MDAnalysis `unwrap` to `center_selection`.
- `center` (`bool`, default `true`): Center `center_selection` in the periodic box.
- `wrap` (`bool`, default `true`): Apply MDAnalysis `wrap` after centering.
- `center_mode` (`"geometry" | "mass"`, default `"geometry"`): Centering mode passed to `center_in_box`.
- `wrap_compound` (`str`, default `"residues"`): Compound grouping passed to `wrap`.
- `align_weights` (`str | null`, default `"mass"`): Weights passed to `fit_rot_trans`.

#### Output options
- `output_dir_name` (`str`, default `"output"`): Output directory under `results/`.
- `output_trajectory_filename` (`str | null`, default `null`): Output trajectory filename. If omitted, the input trajectory suffix is reused, e.g. `production.dcd` becomes `production-centered.dcd`.
- `output_structure_filename` (`str | null`, default `null`): Output frame-0 structure filename. If omitted, defaults to `<trajectory_stem>-centered.pdb`.

Format behavior:
- Input format support follows `MDAnalysis.Universe(topology_path, trajectory_path)`.
- Output trajectory format support follows `MDAnalysis.Writer(output_trajectory_filename, n_atoms)`.
- Output structure format support follows `AtomGroup.write(output_structure_filename)`.
- Common choices include `.xtc`, `.dcd`, `.trr`, `.pdb`, and `.gro`, depending on the installed MDAnalysis writers and the information available in the topology.

### Outputs
For each run, artifacts are written under `runs/<date>/<job-id>/`:
- `results/<output_dir_name>/<output_trajectory_filename>`: centered/aligned trajectory
- `results/<output_dir_name>/<output_structure_filename>`: centered/aligned frame-0 structure
- `results/results.json`: workflow summary
- `results/result_record.json`: standardized CADDForge `Result`
- `stdout.log` / `stderr.log`: runtime logs

## Tool: `extract-frames`

### What it does
The `extract-frames` tool extracts selected frames from an MD trajectory through the standard CADDForge workflow/job interface.

It supports:
- topology and trajectory inputs readable by MDTraj
- configurable MDTraj atom selection for the primary extracted complex outputs
- frame selection by one explicit frame index, all frames from the last N ns, or a Python-style frame slice
- sliced trajectory output plus a first selected-frame complex structure for multi-frame selections
- structure-only output for single-frame selections
- optional one-PDB-per-frame complex export
- optional protein-only outputs when `ligand_resname` is provided
- optional corrected ligand SDF output when both `ligand_resname` and `ligand_template_sdf` are provided
- standardized CADDForge outputs (`results.json`, `result_record.json`)

### External dependencies
This workflow requires `mdtraj` and `rdkit`, which are already dependencies of `caddforge-md`.

### CLI usage
```bash
caddforge run md extract-frames --config extract-frames.yaml
```

### Python API usage
```python
from caddforge.core.execution import load_result_for_job, submit_workflow
from caddforge.plugins import get_tool

spec = get_tool("md", "extract-frames")
cfg = spec.config_cls(
    topology_path="traj-centered.pdb",
    trajectory_path="traj-centered.xtc",
    ligand_template_sdf="template.sdf",
    ligand_resname="UNK",
    last_ns=10.0,
    backend="local",
)

workflow = spec.workflow_cls(cfg)
job = submit_workflow(workflow, cfg, backend_name=cfg.backend)
result = load_result_for_job(job)
if result is not None:
    print(result.status.value, result.paths)
```

### Example config
```yaml
backend: local
topology_path: traj-centered.pdb
trajectory_path: traj-centered.xtc
ligand_template_sdf: template.sdf
ligand_resname: UNK
extract_selection: all
frame_index: null
last_ns: 10.0
frame_slice: null
write_per_frame_pdbs: false
output_dir_name: output
output_trajectory_filename: last-10ns.xtc
output_structure_filename: last-10ns_first_complex.pdb
protein_trajectory_filename: last-10ns_protein.xtc
protein_structure_filename: last-10ns_first_protein.pdb
ligand_sdf_filename: last-10ns_first_ligand.sdf
per_frame_dir_name: frames
```

For a Python-style slice, replace `last_ns` with:
```yaml
last_ns: null
frame_slice:
  start: 100
  stop: 1000
  step: 10
```

### Tool-specific config options

#### Required inputs
- `topology_path` (`Path`): Topology/structure file readable by MDTraj.
- `trajectory_path` (`Path`): Trajectory file readable by MDTraj.

#### Optional ligand inputs
- `ligand_resname` (`str | null`, default `null`): Ligand residue name used to write protein-only outputs and select ligand atoms.
- `ligand_template_sdf` (`Path | null`, default `null`): Template SDF used to write corrected ligand bond orders with trajectory coordinates. Requires `ligand_resname`.

#### Atom selection
- `extract_selection` (`str`, default `"all"`): MDTraj atom selection written to the primary extracted complex outputs. Split protein and ligand outputs still use `ligand_resname`.

#### Frame selection
Configure exactly one of:
- `frame_index` (`int | null`): Extract one frame by zero-based index.
- `last_ns` (`float | null`): Extract all frames whose trajectory time is within the last N ns.
- `frame_slice` (`{start, stop, step} | null`): Python-style frame slice.

#### Output options
- `write_per_frame_pdbs` (`bool`, default `false`): Also write one complex PDB per selected frame.
- `output_dir_name` (`str`, default `"output"`): Output directory under `results/`.
- `output_trajectory_filename` (`str | null`, default `null`): Sliced complex trajectory filename.
- `output_structure_filename` (`str | null`, default `null`): First selected-frame complex PDB filename.
- `protein_trajectory_filename`, `protein_structure_filename`, `ligand_sdf_filename`: Optional filename overrides for split outputs.
- `per_frame_dir_name` (`str`, default `"frames"`): Directory for optional per-frame PDBs.

### Outputs
For each run, artifacts are written under `runs/<date>/<job-id>/`:
- `results/<output_dir_name>/<output_trajectory_filename>`: sliced complex trajectory for multi-frame selections only
- `results/<output_dir_name>/<output_structure_filename>`: first selected-frame complex PDB
- `results/<output_dir_name>/<protein_trajectory_filename>`: optional sliced protein trajectory for multi-frame selections only
- `results/<output_dir_name>/<protein_structure_filename>`: optional first selected-frame protein PDB
- `results/<output_dir_name>/<ligand_sdf_filename>`: optional corrected ligand SDF for the first selected frame
- `results/<output_dir_name>/frame_manifest.json`: selected original frame indices, times, and outputs
- `results/results.json`: workflow summary
- `results/result_record.json`: standardized CADDForge `Result`
- `stdout.log` / `stderr.log`: runtime logs

## Tool: `cluster-representative`

### What it does
The `cluster-representative` tool clusters selected trajectory frames and writes medoids from the requested number of largest clusters (one by default).

It follows the existing standalone script behavior:
- use the last N ns or a Python-style frame slice
- cluster ligand heavy-atom coordinates by default when `ligand_resname` is provided
- write a configurable MDTraj atom selection to the representative complex PDB
- use complete-linkage agglomerative clustering with a configurable distance threshold
- select one medoid from each of the N largest clusters, ranked by population
- write representative complex PDB, optional protein-only PDB, and optional corrected ligand SDF
- standardize CADDForge outputs (`results.json`, `result_record.json`)

### External dependencies
This workflow requires `mdtraj`, `rdkit`, and `scikit-learn`.

### CLI usage
```bash
caddforge run md cluster-representative --config cluster-representative.yaml
```

### Python API usage
```python
from caddforge.core.execution import load_result_for_job, submit_workflow
from caddforge.plugins import get_tool

spec = get_tool("md", "cluster-representative")
cfg = spec.config_cls(
    topology_path="traj-centered.pdb",
    trajectory_path="traj-centered.xtc",
    ligand_template_sdf="template.sdf",
    ligand_resname="UNK",
    last_ns=10.0,
    distance_threshold_nm=0.2,
    n_representatives=3,
    backend="local",
)

workflow = spec.workflow_cls(cfg)
job = submit_workflow(workflow, cfg, backend_name=cfg.backend)
result = load_result_for_job(job)
if result is not None:
    print(result.status.value, result.metrics)
```

### Example config
```yaml
backend: local
topology_path: traj-centered.pdb
trajectory_path: traj-centered.xtc
ligand_template_sdf: template.sdf
ligand_resname: UNK
cluster_selection: null
extract_selection: all
last_ns: 10.0
frame_slice: null
distance_threshold_nm: 0.2
n_representatives: 1
representatives_csv_filename: representatives.csv
output_dir_name: output
output_complex_filename: largest_cluster_representative_complex.pdb
output_protein_filename: largest_cluster_representative_protein.pdb
output_ligand_sdf_filename: largest_cluster_representative_ligand.sdf
clustering_summary_filename: clustering_summary.json
```

For slice-based clustering, replace `last_ns` with:
```yaml
last_ns: null
frame_slice:
  start: 1000
  stop: null
  step: 5
```

### Tool-specific config options

#### Required inputs
- `topology_path` (`Path`): Topology/structure file readable by MDTraj.
- `trajectory_path` (`Path`): Trajectory file readable by MDTraj.

#### Ligand and clustering options
- `ligand_resname` (`str | null`, default `null`): Ligand residue name. If set and `cluster_selection` is omitted, clustering uses `resname <ligand_resname> and not element H`.
- `ligand_template_sdf` (`Path | null`, default `null`): Template SDF used to write corrected ligand bond orders with representative-frame coordinates. Requires `ligand_resname`.
- `cluster_selection` (`str | null`, default `null`): Optional explicit MDTraj atom selection for clustering. If omitted without `ligand_resname`, defaults to `protein and backbone`.
- `extract_selection` (`str`, default `"all"`): MDTraj atom selection written to the representative complex PDB.
- `distance_threshold_nm` (`float`, default `0.2`): Complete-linkage clustering cutoff in nm.
- `n_representatives` (positive integer, default `1`): Maximum number of largest-cluster medoids to write. If fewer clusters exist, all available clusters are written. Population ties use ascending cluster labels; medoid ties use first occurrence in the selected frame order. This does not change the clustering threshold or force a particular cluster count.

#### Frame selection
Configure exactly one of:
- `last_ns` (`float | null`): Use all frames whose trajectory time is within the last N ns.
- `frame_slice` (`{start, stop, step} | null`): Python-style frame slice.

#### Output options
- `output_dir_name` (`str`, default `"output"`): Output directory under `results/`.
- `output_complex_filename` (`str`): Representative complex PDB filename.
- `output_protein_filename` (`str`): Representative protein-only PDB filename.
- `output_ligand_sdf_filename` (`str`): Representative corrected ligand SDF filename.
- `clustering_summary_filename` (`str`): Clustering summary JSON filename.
- `representatives_csv_filename` (`str`, default `"representatives.csv"`): Representative table filename, ending in `.csv` and distinct from other outputs.

### Outputs
For each run, artifacts are written under `runs/<date>/<job-id>/`:
- `results/<output_dir_name>/<output_complex_filename>`: representative complex PDB
- `results/<output_dir_name>/<output_protein_filename>`: optional representative protein-only PDB
- `results/<output_dir_name>/<output_ligand_sdf_filename>`: optional corrected ligand SDF
- `results/<output_dir_name>/<clustering_summary_filename>`: cluster labels, sizes, representative index, and frame time
- `results/results.json`: workflow summary
- `results/result_record.json`: standardized CADDForge `Result`
- `stdout.log` / `stderr.log`: runtime logs

### Multiple representatives and provenance

Set `n_representatives: 3` to export the medoids of the three largest clusters.
With the default count of one, all existing filenames and singular results are preserved.
For requests greater than one, every structure filename gets a rank suffix before its extension:
for example, `complex.pdb` becomes `complex_rank_001.pdb`, `complex_rank_002.pdb`, and
`complex_rank_003.pdb`. This naming applies even if fewer clusters are available.
Each rank contains the complex and any configured protein/ligand outputs from the same frame.

The CSV contains one row per written representative in rank order:

| Column | Meaning |
| --- | --- |
| `rank` | One-based population rank among selected clusters |
| `cluster_label` | Label used in the full clustering summary |
| `cluster_size` | Number of clustered trajectory frames belonging to this cluster |
| `cluster_fraction` | Cluster size divided by the total number of clustered frames |
| `representative_sliced_index` | Zero-based index within the selected trajectory |
| `representative_original_frame_index` | Zero-based index in the original trajectory |
| `representative_time_ns` | Original trajectory timestamp in ns, including after slicing |
| `complex_structure`, `protein_structure`, `ligand_sdf` | Final output paths; optional missing outputs are empty |

The JSON summaries also contain `representatives`, `n_representatives_requested`, and
`n_representatives_written`. Canonical results expose the counts as metrics, the records
in `extra.representatives`, and the CSV as `paths.representatives_csv`. Existing singular
paths and representative metrics continue to describe rank 1. Scratch runs publish final
output paths in both CSV and JSON.

Pipeline outputs `complex`, `protein`, `ligand`, and `clustering_summary` retain their
existing roles. The new `representatives` port exposes the CSV table; `complexes`,
`proteins`, and `ligands` expose its ordered path collections. Bind, for example,
`cluster.complexes` to a downstream structure input and scatter over that input's
configuration field with `chunk_size: 1` to process one medoid per job. Protein and ligand
collections are empty when their exports are not configured.

## Tool: `plain-md`

### What it does
The `plain-md` tool runs a standard OpenMM simulation through the CADDForge workflow/job interface. It supports:
- generic system input as PDB, including protein-only systems
- optional parameterization of one or more non-standard molecules from SDF through OpenFF/SMIRNOFF
- hydrogen addition, TIP3P solvation, neutralization, and NaCl ionic strength setup
- prepared-system mode that skips hydrogen/solvent/ion addition and builds a snug periodic box around existing atoms
- PME nonbonded treatment, HBond constraints, Langevin middle integration, and isotropic or membrane Monte Carlo barostats
- checkpoint-based resume input and automatic checkpoint restart attempts
- standardized CADDForge outputs (`results.json`, `result_record.json`)

### External dependencies
The following runtime dependencies must be available in the execution environment:
- `openmm`
- `openff-toolkit`
- `openmmforcefields`

Recommended conda installation:
```bash
conda install -c conda-forge openmm openff-toolkit openmmforcefields
```

### CLI usage
```bash
caddforge run md plain-md --config plain-md.yaml
```

For Slurm, either set `backend: slurm` in the config or pass `--backend slurm`:
```bash
caddforge run md plain-md --config plain-md.yaml --backend slurm --slurm-time 24:00:00
```

Slurm execution note:
- The workflow runs the OpenMM simulation from a scratch work directory, preferring `SLURM_TMPDIR`.
- The submission directory is not used as the simulation CWD.
- Finished outputs are copied back to `runs/<date>/<job-id>/results/<output_dir_name>/`.

### Python API usage
```python
from caddforge.core.execution import load_result_for_job, submit_workflow
from caddforge.plugins import get_tool

spec = get_tool("md", "plain-md")
cfg = spec.config_cls(
    system_pdb_path="system.pdb",
    molecule_sdf_paths=["ligand.sdf", "cofactor.sdf"],  # optional
    simulation_time_ns=1.0,
    backend="local",  # or "slurm"
)

workflow = spec.workflow_cls(cfg)
job = submit_workflow(workflow, cfg, backend_name=cfg.backend)
result = load_result_for_job(job)
if result is not None:
    print(result.status.value, result.metrics)
```

### Example config
```yaml
backend: local
system_pdb_path: system.pdb
molecule_sdf_paths:
  - ligand.sdf
  - cofactor.sdf
ph: 7.4
system_preparation_mode: solvate
padding: 1.0
ionic_strength: 0.15
prepared_box_margin: 0.0
nonbonded_cutoff: 1.0
barostat: isotropic
pressure: 1.0
membrane_surface_tension: 0.0
membrane_xy_mode: isotropic
membrane_z_mode: free
barostat_frequency: 25
temperature: 300.0
step_size_fs: 4
simulation_time_ns: 1.0
trajectory_filename: traj.xtc
report_interval: 10000
checkpoint_interval: 5000
output_dir_name: output
forcefields: []
seed: null
checkpoint_path: null
original_trajectory_path: null
max_restarts: 10
```

### Tool-specific config options

#### Required inputs
- `system_pdb_path` (`Path`): System to simulate in PDB format.

#### Optional non-standard molecule inputs
- `molecule_sdf_paths` (`list[Path]`, default `[]`): SDF files for molecules that are not covered by the selected OpenMM force fields.

Input behavior:
- Leave `molecule_sdf_paths` empty for protein-only or otherwise fully force-field-covered systems.
- Add one SDF per ligand, cofactor, or other non-standard small molecule that needs SMIRNOFF template generation.
- For one molecule, `molecule_sdf_paths: ligand.sdf` is also accepted and normalized to a one-item list.
- Multi-record SDF files are accepted; all molecules loaded from each SDF are registered with the SMIRNOFF template generator.

#### Simulation setup options
- `ph` (`float`, default `7.4`): pH used for hydrogen addition.
- `system_preparation_mode` (`"solvate" | "prepared"`, default `"solvate"`): Use `"solvate"` for the existing hydrogen/solvent/ion setup, or `"prepared"` for already complete systems that should be simulated as provided.
- `padding` (`float`, default `1.0`): Solvent padding in nm, used only in `"solvate"` mode.
- `ionic_strength` (`float`, default `0.15`): NaCl concentration in molar, used only in `"solvate"` mode.
- `prepared_box_margin` (`float`, default `0.0`): Extra margin in nm added on each side of the coordinate bounding box in `"prepared"` mode.
- `nonbonded_cutoff` (`float`, default `1.0`): Nonbonded cutoff in nm.
- `barostat` (`"isotropic" | "membrane"`, default `"isotropic"`): Pressure-coupling mode. `"isotropic"` uses OpenMM `MonteCarloBarostat`; `"membrane"` uses `MonteCarloMembraneBarostat`.
- `pressure` (`float`, default `1.0`): Pressure in bar.
- `membrane_surface_tension` (`float`, default `0.0`): Surface tension in bar*nm for membrane barostat mode.
- `membrane_xy_mode` (`"isotropic" | "anisotropic"`, default `"isotropic"`): X/Y scaling mode for membrane barostat mode.
- `membrane_z_mode` (`"free" | "fixed" | "constant_volume"`, default `"free"`): Z scaling mode for membrane barostat mode.
- `barostat_frequency` (`int`, default `25`): Frequency in MD steps for attempted barostat moves.
- `temperature` (`float`, default `300.0`): Temperature in K.
- `forcefields` (`list[str]`, default `[]`): Optional OpenMM force field XML names. Empty uses `amber/ff14SB.xml`, `amber14/tip3p.xml`, `amber/tip3p_HFE_multivalent.xml`, and `amber/phosaa10.xml`.

Prepared-system mode uses the input atoms exactly: it does not add hydrogens, water, or ions. It computes an orthorhombic periodic box from the existing coordinates plus `prepared_box_margin`, shifts atoms inside that box, and fails early if any box length is too small for `2 * nonbonded_cutoff`.

Membrane barostat mode assumes the membrane lies in the XY plane and the membrane normal is Z. A typical membrane setup uses:

```yaml
barostat: membrane
membrane_surface_tension: 0.0
membrane_xy_mode: isotropic
membrane_z_mode: free
barostat_frequency: 25
```

#### Runtime options
- `step_size_fs` (`int`, default `4`): Integrator step size in fs.
- `simulation_time_ns` (`float`, default `1.0`): Simulation length in ns.
- `trajectory_filename` (`str`, default `"traj.xtc"`): Trajectory filename under the output directory.
- `report_interval` (`int`, default `10000`): State/trajectory reporting interval in steps.
- `checkpoint_interval` (`int`, default `5000`): Checkpoint writing interval in steps.
- `seed` (`int | null`, default `null`): Optional random seed for barostat, integrator, and initial velocities.
- `checkpoint_path` (`Path | null`, default `null`): Optional checkpoint to resume from.
- `original_trajectory_path` (`Path | null`, default `null`): Optional previous trajectory copied to the new output directory before appending during checkpoint resume. If omitted, checkpoint resume starts a new trajectory file instead of appending.
- `max_restarts` (`int`, default `10`): Maximum automatic restarts from the job-local checkpoint after runtime failure.
- `output_dir_name` (`str`, default `"output"`): Simulation output directory under `results/`.

### Outputs
For each run, artifacts are written under `runs/<date>/<job-id>/`:
- `results/<output_dir_name>/`: OpenMM simulation outputs:
  - `<system_stem>_prepared.pdb`
  - `<system_stem>_minimized.pdb`
  - trajectory file, by default `traj.xtc`
  - `checkpoint.chk`
  - `final.xml`
- `results/results.json`: workflow summary
- `results/result_record.json`: standardized CADDForge `Result`
  - exposes `prepared_structure`, `minimized_structure`, `trajectory`, `checkpoint`, and
    `output_dir` as pipeline-ready result paths
- `stdout.log` / `stderr.log`: command and runtime logs

## Tool: `openbpmd`

### What it does
The `openbpmd` tool runs an OpenBPMD workflow in three steps:
1. Preprocess a protein-ligand system into AMBER solvated inputs (`solvated.rst7`, `solvated.prm7`).
2. Run OpenBPMD simulations with standard CADDForge backend behavior (local/slurm).
3. Run averaged multi-replica analysis and write summary outputs.

This workflow follows the OpenBPMD implementation from the Gervasiolab repository:
- <https://github.com/Gervasiolab/OpenBPMD>

### References
- A. J. Clark, P. Tiwary, K. Borrelli, S. Feng, E. Miller, R. Abel, R. A. Friesner, B. J. Berne, J. Chem. Theory Comput. 2016, 12, 6, 2990–2998, DOI: <https://doi.org/10.1021/acs.jctc.6b00201>
- D. Lukauskis, M. L. Samways, S. Aureli, B. P. Cossins, R. D. Taylor, F. L. Gervasio, J. Chem. Inf. Model. 2022, 62, 23, 6209–6216, DOI: <https://doi.org/10.1021/acs.jctc.6b00201>

### External dependencies
The following external tools must be available in the runtime environment:
- AmberTools commands used in preprocessing (`pdb4amber`, `antechamber`, `parmchk2`, `tleap`)

Also install OpenMM via conda (required by OpenBPMD):
```bash
conda install -c conda-forge openmm
```

Bundled script note:
- `caddforge-md` vendors the upstream `openbpmd.py` script from the Gervasiolab OpenBPMD repository under `src/cadd_md/openbpmd/vendor/openbpmd.py` (MIT licensed).
- You can still override this via `openbpmd_script_path` if needed.

### CLI usage
```bash
caddforge run md openbpmd --config openbpmd.yaml
```

### Python API usage
```python
from caddforge.core.execution import load_result_for_job, submit_workflow
from caddforge.plugins import get_tool

spec = get_tool("md", "openbpmd")
cfg = spec.config_cls(
    complex_pdb_path="complex.pdb",
    ligand_path="ligand.sdf",
    lig_resname="MOL",
    nreps=10,
    hill_height=0.05,
    deposit_interval_ps=1.0,
    backend="local",
)

workflow = spec.workflow_cls(cfg)
job = submit_workflow(workflow, cfg, backend_name=cfg.backend)
result = load_result_for_job(job)
if result is not None:
    print(result.status.value, result.metrics)
```

### Example config
```yaml
backend: local
complex_pdb_path: complex.pdb
ligand_path: ligand.sdf
lig_resname: MOL
ligand_net_charge: null
input_system_kind: dry
hydrogen_policy: preserve
solvent_padding_angstrom: 10.0
ionic_strength_molar: 0.0
nreps: 10
hill_height: 0.05
deposit_interval_ps: 1.0
log_sync_interval_seconds: 30
output_dir_name: output
analysis_dir_name: analysis
```

### Tool-specific config options

#### Input modes
Provide `ligand_path` together with either `complex_pdb_path` (protein and exactly
one ligand) or `protein_pdb_path` (protein only). The authoritative ligand must be
one explicit-hydrogen 3D molecule in SDF/SD, MOL, or MOL2 format. Its graph, charge,
stereochemistry, and individual atom coordinates must match the complex. OpenBPMD
neither generates a replacement conformer nor aligns mismatched coordinates.
PDB-only ligand perception is no longer supported; supply the original chemistry
template with coordinates from the same selected frame.

`input_system_kind: dry` rejects inherited waters and ions. A crystal `CRYST1`
record alone does not declare a simulation box. `pocket_hydrated` requires a
nonperiodic complex and its `structure_manifest_path` hydration certificate,
including matching checksum, water inventory, and sphere metadata. Full periodic
boxes and uncertified solvated complexes are unsupported inputs. Inducia supplies
the certified pocket, selected-frame ligand SDF, and protein chemistry records.

Runtime, preparation, and validation outputs are copied without filesystem
metadata operations and verified by SHA256. This supports CIFS destinations that
reject metadata updates. An incomplete runtime or preparation copy marks the job
failed and preserves scratch evidence; a preceding scientific failure remains in
the final error. Failure to retain the separate validation-report copy is logged
and preserves scratch and the report path without changing the job status.

#### OpenBPMD run options (mirroring `openbpmd.py`)
- `lig_resname` (`str`, default `"MOL"`): ligand residue name (`-lig_resname`)
- `nreps` (`int`, default `10`): number of OpenBPMD repeats (`-nreps`)
- `hill_height` (`float`, default `0.05`): metadynamics hill height (`-hill_height`)
- `deposit_interval_ps` (`float`, default `1.0`): bias deposit interval (`--deposit_interval_ps`)
- `log_sync_interval_seconds` (`int`, default `30`): interval in seconds for best-effort runtime-log
  and structured progress synchronization; accepted range is 1 through 3600.
- `output_dir_name` (`str`, default `"output"`): simulation output directory under `results/`
- `openbpmd_script_path` (`Path | null`, default `null`): optional path to OpenBPMD script. If `null`/omitted, the bundled vendored script is used.

#### Preprocessing options
- `ligand_net_charge` (`int | null`, default `null`): optional formal-charge assertion.
  The charge is derived from the authoritative ligand, including when the assertion
  is zero. Generated charges use the bounded normalization policy described above;
  the final component sum must agree within 1e-8 e.
- `hydrogen_policy` (`preserve` or `reduce`, default `preserve`): protein-only
  `pdb4amber` runs without `--reduce` by default. Explicit Reduce still must preserve
  pinned states and supplied coordinates. Internal Reduce failure is fatal even
  when `pdb4amber` exits zero; its diagnostics are retained.
- `solvent_padding_angstrom` (`float`, default `10.0`): positive padding for the one
  newly built TIP3P box.
- `ionic_strength_molar` (`float`, default `0.0`): nonnegative NaCl pair concentration.
  Pair count is round-half-up(concentration × Avogadro constant × actual box volume
  in liters); neutralizing ions are additional. The report distinguishes requested
  concentration, realized salt-pair concentration, total ionic strength, and counterions.

Ion placement replaces only fresh bulk waters, with a recorded deterministic seed
and minimum separations. Retained pocket waters keep their identities and supplied
coordinates. Fresh waters inside the certified sphere are removed so solvation
cannot refill sampled vacancies. Amber's common solvation translation is removed
from the whole system before writing the final restart. The final topology/restart
is reloaded and checked for coordinates, box, state and bond inventory, ligand,
ion counts, and charge; no large PDB round trip is used.

OpenBPMD checks PDB identifiers and connectivity before preparation and after `pdb4amber`.
Atom serials must be unique decimal values; residue addresses and bond endpoints must be
unambiguous. Oversized PDBs fail at this boundary. Splitting retains chain boundaries and
remapped covalent links, with reversible identifiers; splitting a covalently attached ligand
is rejected until that parameterization is supported.

### Outputs
For each run, artifacts are written under `runs/<date>/<job-id>/`:
- `results/preparation/`: durable Amber topology/restart, preparation commands and
  outputs, input mappings, state/link and ligand evidence, and Reduce diagnostics
  when present. Published as `Result.paths["preparation_dir"]`, including failures.
- `results/preparation/preparation_summary.json`: versioned preparation policy,
  requested/realized solvent composition, validation evidence, and output checksums;
  published as `Result.paths["preparation_summary"]`. Scratch is retained if copying
  these diagnostics fails.
- `results/structure_validation/`: durable input/export identity and validation reports, also
  published as `Result.paths["structure_validation"]` on preparation failure
- `results/output/`: OpenBPMD simulation outputs (replica folders, etc.)
- `results/output/runtime_summary.json`: executed step counts, software/source versions,
  input hashes and runtime-coordinate hashes; published as `Result.paths["runtime_summary"]`
  by the bundled script. Failed dynamics retain this summary and partial outputs.
- `results/output/minimized_system.npz` and `equil_system.npz`: full-precision positions
  and matching cell vectors in nanometers, bound to the original parameter topology's
  atom order by SHA256. These replace runtime PDB intermediates, including for systems
  above 99,999 atoms. Load with `cadd_md.openbpmd.runtime_coordinates.read_coordinates`
  and `openmm.app.AmberPrmtopFile("results/preparation/solvated.prm7").topology`. Analysis uses that
  topology directly. External custom scripts must provide their own compatible diagnostics.
- `results/analysis/`: averaged multi-replica analysis outputs:
  - `results.csv`
  - `multi_rep_timeseries.csv`
  - `multi_rep_CompScore.png`
  - `multi_rep_PoseScore.png`
  - `multi_rep_ContactScore.png`
- `results/results.json`: workflow summary
- `results/result_record.json`: standardized CADDForge `Result`
- `results/progress.json`: best-effort current stage and replica progress snapshot
- `stdout.log` / `stderr.log`: command and runtime logs
Protein chemistry can be supplied explicitly with the paired `protein_states_path`
and `protein_links_path` JSON records. Their checksum, heavy-atom identities, and
link geometry must match the input. OpenBPMD maps them through splitting and
`pdb4amber` renumbering, restores Amber hydrogen aliases and cross-link CONECT
records, and rejects unsupported non-disulfide links. The installed AmberTools
TER-number convention is adapted only after strict atom/connectivity validation;
the raw output and export mapping are retained in `structure_validation`.

## Tool: `grand`

### What it does
The `grand` tool runs Grand Canonical Monte Carlo (GCMC) water sampling in a ligand-centred pocket region and then selects a single representative protein-ligand-water structure from the resulting simulation, all within one job.

It runs two stages in sequence:

1. GCMC water-sampling stage:
   - prepares the input protein-ligand complex (repairs hydrogens within preserved protein states, solvates with padding and ionic strength, restrains solute heavy atoms, and adds ghost waters)
   - runs GCMC insertion/deletion sampling within a ligand-centred spherical region using the configured pocket sphere radius, initialization passes, and maximum production passes
   - writes the final-state checkpoint before any trajectory readback or post-processing
   - produces a trajectory, a per-frame ghost-water state file, a run checkpoint, and a hydrated-complex structure
   - writes a stage summary recording the produced output paths, the trajectory frame count, and a success/failure indication
   - treats trajectory post-processing as non-strict: a post-processing failure is recorded in the summary while the completed checkpoint and trajectory are preserved and the run continues

2. Representative-selection stage:
   - loads the produced trajectory over the configured frame analysis interval
   - uses the per-frame ghost-water state file to exclude, for each analysed frame, that frame's inactive ghost waters from all downstream selection and output
   - deterministically selects exactly one representative frame, so repeated runs on the same trajectory and configuration select the same representative frame index
   - writes a representative complex structure containing the ligand residue, the protein atoms, and the active non-ghost water residues present at the selected frame, plus an optional protein-only structure
   - writes an optional representative ligand SDF whose bond orders and chemistry come from the template and whose coordinates come from the selected frame
   - writes a selection summary recording the selected representative frame index and a success/failure status

It runs on both the local and Slurm backends. Under an allocated Slurm backend it executes in a node-local scratch directory and copies results back to the shared job directory only after the run finishes, while periodically and resiliently syncing the runtime log files so progress is visible during long runs.

Trajectory postprocessing reconstructs whole bonded components before imaging
them near the ligand and aligning protein C-alpha atoms. Recorded bonds, including
cross-chain links, take precedence over residue-name classification; supported
Amber protonation states and terminal caps remain part of their bonded protein.
Ghost waters are displaced only after imaging. The serialized aligned view must
match the raw solute's bonded geometry within `0.0002 nm`; failures identify the
frame, bond, and maximum residual. The raw-periodic fallback is explicitly checked
against its raw frame and reconstructed into whole molecules for structure exports.

### External dependencies
The following runtime dependencies must be available in the execution environment:
- a GCMC water-sampling engine that performs grand canonical insertion/deletion moves
- an MD / force-field runtime used to prepare the system, build force fields, and integrate dynamics

These engine libraries are imported lazily inside the GCMC stage, so the rest of the tool (configuration, validation, selection, and result handling) can be exercised without them installed. Representative selection additionally uses `mdtraj` and `rdkit`, which are already dependencies of `caddforge-md`.

### CLI usage
```bash
caddforge run md grand --config grand_conf.yaml
```

Print a templated example configuration with:
```bash
caddforge show-config md grand
```

For Slurm, either set `backend: slurm` in the config or pass `--backend slurm`:
```bash
caddforge run md grand --config grand_conf.yaml --backend slurm
```

Slurm execution note:
- The workflow runs both stages from a scratch work directory, preferring `SLURM_TMPDIR`.
- The submission directory is not used as the working directory.
- Runtime logs are synced back periodically during the run so progress is visible.
- Finished outputs are copied back to `runs/<date>/<job-id>/results/<output_dir_name>/` only after the run completes.

### Python API usage
```python
from caddforge.core.execution import load_result_for_job, submit_workflow
from caddforge.plugins import get_tool

spec = get_tool("md", "grand")
cfg = spec.config_cls(
    complex_pdb_path="complex.pdb",
    ligand_template_sdf="ligand_template.sdf",
    ligand_resname="LIG",
    sphere_radius_angstrom=10.0,
    n_ghost_waters=150,
    max_production_passes=100000,
    backend="local",  # or "slurm"
)

workflow = spec.workflow_cls(cfg)
job = submit_workflow(workflow, cfg, backend_name=cfg.backend)
result = load_result_for_job(job)
if result is not None:
    print(result.status.value, result.metrics)
    print(result.paths["representative_complex"])
```

### Example config
```yaml
backend: local
complex_pdb_path: complex.pdb
ligand_template_sdf: ligand_template.sdf
ligand_resname: LIG
sphere_radius_angstrom: 10.0
n_ghost_waters: 150
nonbonded_cutoff_nm: 0.9
temperature_kelvin: 300.0
initialization_passes: 7500
max_production_passes: 100000
ph: 7.4
padding_nm: 1.2
ionic_strength_molar: 0.15
ligand_selection: null
align_selection: null
frame_analysis_interval: 1
log_sync_interval_seconds: 30
output_dir_name: output
trajectory_filename: gcmc_trajectory.dcd
ghost_water_filename: ghost_water_state.txt
checkpoint_filename: gcmc_checkpoint.chk
hydrated_complex_filename: hydrated_complex.pdb
gcmc_summary_filename: gcmc_summary.json
representative_complex_filename: representative_complex.pdb
representative_pocket_complex_filename: representative_pocket_complex.pdb
representative_protein_filename: representative_protein.pdb
representative_ligand_filename: representative_ligand.sdf
selection_summary_filename: selection_summary.json
```

### Tool-specific config options

#### Required inputs
- `complex_pdb_path` (`Path`): Input protein-ligand complex in PDB format. Must exist and have a `.pdb` suffix.
- `ligand_template_sdf` (`Path`): One connected ligand chemistry template, as a 2D or 3D SDF with explicit or implicit hydrogens. Supplies connectivity, bond orders, formal charges, isotopes, protonation and specified stereochemistry. Template coordinates and atom order need not match the complex. Must exist and have a `.sdf` or `.sd` suffix.
- `ligand_resname` (`str`): Residue name identifying the ligand in the input complex, used by both stages. Must be 1 to 4 characters with no whitespace or path separators.

Input paths are normalized to absolute paths after validation, so execution from a scratch directory resolves inputs correctly.

GRAND binds chemistry to the ligand in `complex_pdb_path` using element identity,
connectivity evidence and bound bond geometry. It preserves every supplied ligand
coordinate and adds only missing hydrogens from that geometry. It never aligns the
template or uses template coordinates for hydrogen placement. Missing heavy atoms,
contradictory protonation/stereochemistry, covalent ligand attachments and genuinely
ambiguous chemical assignments fail before sampling. Equivalent symmetry and
resonance assignments are handled deterministically. The search is bounded; exceeding
the limit produces an explicit diagnostic rather than choosing an arbitrary match.

No existing config change is required. If mapping is ambiguous, inspect
`ligand_binding.json` and check the intended protonation and input connectivity.
An optional `ligand_atom_map_path` (`Path | null`, default `null`; pipeline input
`ligand_atom_map`) can resolve an otherwise valid ambiguity. Its JSON schema is:

```json
{
  "schema_version": 1,
  "complex_sha256": "<SHA-256 of the input PDB bytes>",
  "template_sha256": "<SHA-256 of the original SDF bytes>",
  "atom_pairs": [
    {"complex_atom_serial": 12, "template_atom_index": 0},
    {"complex_atom_serial": 13, "template_atom_index": 1}
  ]
}
```

The pairs must cover every ligand heavy atom exactly once. PDB serials are the
original serial numbers; SDF indices are zero-based original record atom indices.
Supplied hydrogens may also be mapped when explicitly present in the SDF. The map
must be regenerated when either input file changes and cannot override geometry,
connectivity, charge or stereochemical validation.

Optional paired `protein_states_path` and `protein_links_path` records preserve
upstream protein chemistry. Without records, `residue_states` can specify named
states using `chain:number[:insertion]` keys; names and explicit hydrogens supply
the remaining evidence. GRAND preserves these decisions during hydrogen repair.
Use the hydrogen-repair tool's explicit `assign` policy before GRAND if automatic
protonation is needed. State/link records are embedded in the lossless topology
companion and published with the mapped representative PDB.

#### GCMC water-sampling options
- `sphere_radius_angstrom` (`float`, default `10.0`): Radius in angstrom of the ligand-centred GCMC sampling region. Must be greater than 0.
- `n_ghost_waters` (`int`, default `150`): Number of ghost waters added for GCMC insertion/deletion. Must be at least 1.
- `nonbonded_cutoff_nm` (`float`, default `0.9`): Nonbonded interaction cutoff in nm. Must be greater than 0.
- `temperature_kelvin` (`float`, default `300.0`): Simulation temperature in K. Must be greater than 0.
- `initialization_passes` (`int`, default `7500`): Number of GCMC equilibration passes before production sampling. Must be at least 1.
- `max_production_passes` (`int`, default `100000`): Maximum number of GCMC production passes. Must be at least 1.
- `ph` (`float`, default `7.4`): protonation metadata; it does not override preserved states or opt into automatic assignment.
- `padding_nm` (`float`, default `1.2`): Solvent box padding in nm. Must be greater than 0.
- `ionic_strength_molar` (`float`, default `0.15`): Ionic strength in molar used when solvating. Must be greater than or equal to 0.

#### Representative-selection options
- `ligand_selection` (`str | null`, default `null`): Atom selection identifying the ligand for selection output. When omitted, it is derived from `ligand_resname` as `resname <ligand_resname>`.
- `align_selection` (`str | null`, default `null`): Optional atom selection used to align frames. Blank values are normalized to `null`.
- `frame_analysis_interval` (`int`, default `1`): Analyse every Nth trajectory frame during representative selection. Must be a positive integer.

#### Log-syncing options
- `log_sync_interval_seconds` (`int`, default `30`): Interval in seconds between periodic log-sync copies from scratch to the job directory during allocated-runtime execution. Must be within 1 to 3600 seconds inclusive.

#### Output options
- `output_dir_name` (`str`, default `"output"`): Output directory under `results/`. Must be a plain name with no path separator.
- `trajectory_filename` (`str`, default `"gcmc_trajectory.dcd"`): Produced GCMC trajectory filename.
- `ghost_water_filename` (`str`, default `"ghost_water_state.txt"`): Per-frame ghost-water state filename.
- `checkpoint_filename` (`str`, default `"gcmc_checkpoint.chk"`): Final-state GCMC checkpoint filename.
- `hydrated_complex_filename` (`str`, default `"hydrated_complex.pdb"`): Hydrated-complex structure filename.
- `gcmc_summary_filename` (`str`, default `"gcmc_summary.json"`): GCMC stage summary JSON filename.
- `representative_complex_filename` (`str`, default `"representative_complex.pdb"`): Representative complex PDB filename.
- `representative_pocket_complex_filename` (`str`, default `"representative_pocket_complex.pdb"`): Compact protein/ligand/pocket-water PDB filename.
- `representative_protein_filename` (`str`, default `"representative_protein.pdb"`): Representative protein-only PDB filename.
- `representative_ligand_filename` (`str`, default `"representative_ligand.sdf"`): Representative ligand SDF filename.
- `selection_summary_filename` (`str`, default `"selection_summary.json"`): Selection stage summary JSON filename.

All output directory and filename fields must be plain names with no `/` or `\` path separator.

### Outputs
For each run, artifacts are written under `runs/<date>/<job-id>/`:
- `results/<output_dir_name>/`: stage outputs, by default:
  - `gcmc_trajectory.dcd`: the GCMC trajectory
  - `gcmc_topology.json`: lossless topology companion with stable atom/residue indices, bonds,
    original labels, atom-order checksum, and reference box vectors; available through the
    `trajectory_topology` output port and result path
  - `gcmc_topology.pdb`: legacy viewing output; large-system identifiers may wrap, so this
    file is not used for trajectory selection
  - `*.validation.json`: input, topology, and representative-export validation reports
  - `resolved_ligand.sdf`: explicit-hydrogen chemistry in the original bound pose, exposed as `resolved_ligand`
  - `ligand_binding.json`: input checksums, mapping, hydrogen counts and stage-specific failure diagnostics, exposed as `ligand_binding`
  - `input_structure.identity.json`: input identities captured before preparation renumbering
  - `ghost_water_state.txt`: per-frame inactive ghost-water residue indices (one row per trajectory frame)
  - `gcmc_checkpoint.chk`: final-state checkpoint
  - `hydrated_complex.pdb`: hydrated-complex structure
  - `gcmc_summary.json`: GCMC stage summary
  - `representative_complex.pdb`: legacy representative with all active waters; omitted with an export-capacity diagnostic when it cannot fit PDB
  - `representative_pocket_complex.pdb`: protein, one ligand, and complete active waters inside the sampled sphere, with no periodic box
  - `representative_pocket_complex.hydration.json`: sphere, water/ion inventory, source frame, archive checksum, and reversible export identities
  - `gcmc_raw_trajectory.dcd` and `periodic_archive.json`: recoverable periodic coordinates and the matching topology, ghost rows, box vectors, checksums, and raw/view frame correspondence
  - `representative_protein.pdb`: optional representative protein-only structure
  - `representative_ligand.sdf`: optional template-corrected representative ligand
  - `selection_summary.json`: selection stage summary
- `results/results.json`: human-readable workflow summary
- `results/result_record.json`: standardized CADDForge `Result`
  - `metrics` include `n_frames_sampled` and `representative_frame_index`
  - `paths` expose `representative_complex`, `trajectory`, `checkpoint`, `output_dir`, and, when produced, `representative_protein` and `representative_ligand`
- `results/progress.json`: best-effort current stage, pass count, acceptance rate, and ETA snapshot
- `stdout.log` / `stderr.log`: runtime logs

The scientific handoff uses `representative_pocket_complex` and `pocket_manifest`
(pipeline ports `pocket_complex` and `pocket_manifest`). Water membership uses the
selected raw frame's periodic cell and the same unweighted ligand-heavy-atom
center used by sampling, with oxygen distance at or below the radius. Water counts
refer to molecules; atom indices are recorded separately. Ordinary Na/K/Cl ions,
remote active waters, and inactive ghosts are excluded. Unknown structural
components, including metals and cofactors without a supported contract, fail
with an inventory. Missing ghost rows fail; an empty pocket is valid.

New topology and periodic archive records use schema 2 and contain relative,
checksum-verified references to the resolved ligand and binding report. Representative
export uses this saved chemistry with the selected frame coordinates. The archive
can be relocated and reselected without the original template; keep all referenced
artifacts together. Legacy schema 1 topology/archive records retain their original
mapping semantics and require their original template for ligand export. Pocket
certificate versions are independent of topology/archive versions.

Pocket manifest schema 2 records `bonded_components_near_ligand` imaging. Each
bonded protein component is reconstructed from raw periodic bond vectors and
placed near the ligand; complete waters are imaged into that same frame. The
protein-only representative uses the pocket coordinates too. This remains valid
when optional trajectory-view processing fails and selection uses raw frames.
Schema 1 pocket certificates must be regenerated by reselecting a valid periodic
archive; they are rejected as downstream scoring inputs.

`raw_periodic_trajectory`, `trajectory_topology`, `ghost_water_states`, and
`periodic_archive` preserve the full system, including ions and ghosts. Remove
the water residues listed in the selected raw frame's ghost row to obtain its
active system. Indices refer to the topology companion, never PDB labels. The
`aligned_trajectory` is a separate box-free view; it is absent on postprocessing
failure. The legacy `trajectory` key retains its documented raw fallback, whose
role is recorded explicitly in the archive. Archive file references are relative
and survive scratch copy-back. Inducia requires the pocket contract; older GRAND
results must be reselected from valid raw artifacts and cannot fall back to a
bulk-water PDB.
## Compressed inputs

Declared structure, molecule, topology, trajectory, checkpoint, and restart inputs may use a
trailing `.gz`. CADDForge expands wrappers in runtime scratch before the MD workflow runs; plain
inputs and result formats keep their existing behavior.

## Tool: `inducia`

GRAND accepts chemistry-only ligand templates. Inducia's known-binder template
still requires a positioned 3D pose for PheSA, and OpenBPMD's authoritative ligand
still requires explicit hydrogens and the same frame as its complex. GRAND exports
that strict downstream representation automatically.

### What it does

`inducia` is a composite induced-fit docking and binding-pose stability workflow built from the
open-source PheSA, gnina, GRAND, and OpenBPMD tools. One Inducia job accepts exactly one 3D query
ligand. Use CADDForge pipeline scatter for a multi-record SDF library.

The two pose-generation modes are:

- `template`: align the query to a correctly positioned known binder with 21 fixed PheSA
  pharmacophore-point weights, cluster the alignments, detect clashing protein residues, and
  minimize each representative with gnina using those residues as flexible.
- `auto`: perform full flexible gnina docking in a configured search box. Flexible residues
  are selected within a radius around the query and capped at a configurable count.

Every surviving gnina pose proceeds through mandatory GRAND water sampling. GRAND receives the
paired gnina complex and normalized ligand-pose SDF so both inputs use the same coordinate frame.
Hydrated complexes are then ordered by descending gnina CNNscore, where higher is better, and at most
`openbpmd.max_input_complexes` are passed to OpenBPMD. The final pose ranking is independent of that
shortlist: it uses ascending OpenBPMD CompScore, where lower is better.

### External dependencies and compute requirements

Install compatible releases of `caddforge`, `caddforge-descriptors`, and `caddforge-docking`; the
minimum versions are declared in `pyproject.toml`. Scientific execution additionally requires:

- PheSA: the `pyocl`/OpenChemLib runtime supplied by `pyocl-idor>=0.8.3` through
  `caddforge-descriptors`, requiring Java 21 or newer and a configured `JAVA_HOME`.
- Docking: a `gnina` executable, or another executable path supplied through
  `gnina.gnina_executable`.
- GRAND: `grand`, OpenMM, OpenMMTools, OpenMMForceFields, the OpenFF toolkit and units packages, and
  MDTraj.
- OpenBPMD preprocessing: the AmberTools commands `pdb4amber`, `antechamber`, `parmchk2`, and
  `tleap`. The OpenBPMD Python script is bundled, while `openbpmd.openbpmd_script_path` can select a
  compatible replacement.

PheSA alignment and most preparation work can run on CPUs. Gnina CNN scoring, GRAND, and OpenBPMD
normally benefit substantially from GPU execution. `gnina.device` selects a device and
`gnina.no_gpu` requests CPU-only gnina execution. Local child concurrency is deliberately bounded
by `max_local_workers` to avoid accidental GPU oversubscription. For cluster execution, configure
`backend: slurm` and the common Slurm fields described below; required software must be available
through the submitted environment.

### CLI usage

Generate the complete default configuration schema:

```bash
caddforge show-config md inducia
```

Run either checked-in single-ligand example after replacing its placeholder paths:

```bash
caddforge run md inducia --config examples/inducia_template.yaml
caddforge run md inducia --config examples/inducia_auto.yaml
```

Scatter a multi-record query SDF into one Inducia job per ligand:

```bash
caddforge pipeline run examples/inducia_pipeline_scatter.yaml
```

### Python API usage

```python
from caddforge.core.execution import load_result_for_job, submit_workflow
from caddforge.plugins import get_tool

spec = get_tool("md", "inducia")
config = spec.config_cls.from_yaml("examples/inducia_template.yaml")
workflow = spec.workflow_cls(config)
job = submit_workflow(
    workflow,
    config,
    backend_name=config.backend,
    input_ports=spec.plugin.input_ports,
)
result = load_result_for_job(job)
if result is not None:
    print(result.status.value)
    if result.status.value == "COMPLETED":
        print(result.metrics["best_comp_score"])
        print(result.paths["best_complex"])
```

### Inputs and mode selection

Both modes require:

- `mode`: exactly `template` or `auto`.
- `receptor_pdb_path`: existing receptor PDB. Pipeline input port: `receptor`.
- `query_ligand_sdf_path`: existing SDF/SD containing exactly one readable molecule with one 3D
  conformer. Pipeline input port: `query_ligand`. A multi-record file is valid only at a pipeline
  scatter boundary, which creates one single-ligand job per shard.

Template mode additionally requires:

- `template_ligand_sdf_path`: one 3D known binder already in the target binding coordinates.
  Pipeline input port: `template_ligand`. It is forbidden in auto mode.

Auto mode requires exactly one of these search-space strategies:

- `autobox_ligand_path`: one positioned 3D ligand used by gnina to define the search box. Pipeline
  input port: `autobox_ligand`.
- `box_center` together with `box_size`: three finite XYZ values for the center and three positive
  dimensions in angstrom. Neither field may be supplied without the other.

`autobox_ligand_path`, `box_center`, and `box_size` are forbidden in template mode because the known
binder is the autobox reference.

### Protocol and child-job behavior

Template mode runs these phases in order:

1. Submit PheSA alignment children for weights 0.00 through 1.00 in increments of 0.05.
2. Cluster successful alignments by direct heavy-atom RMSD without superposition.
3. Keep the highest-total-similarity representative per cluster, breaking ties by lower PheSA
   weight, and detect its clashing protein residues.
4. Exclude covalently linked residues and Ala/Gly/Pro from flexibility, then run one gnina
   minimization child per representative. A representative with no permitted flexible residues uses the rigid
   complex builder and has the same normalized complex/protein/ligand artifact contract.

Auto mode runs one full gnina docking child and retains every successful output mode up to
`gnina.num_modes`. Inducia selects side-chain heavy atoms strictly within
`gnina.auto_flex_radius_angstrom` of any query-ligand atom, excludes linked
residues and Ala/Gly/Pro, then retains at most
`gnina.auto_max_flexible_residues`. The resulting explicit `flexres` list
is passed to GNINA. Empty lists use rigid complex assembly in both modes.

Disulfides are established once from explicit records or unique SG pairs at
1.8–2.3 Å; both partners are protected. Known other covalent links receive the
same protection. Clash reports still include protected residues, with exclusions
and partners recorded separately. GNINA outputs must retain each link endpoint
and valid geometry before entering GRAND; validated links are rebuilt using the
output's current atom serials. Per-pose validation reports include measured
distances and rejection reasons.

The auto radius matches GNINA's side-chain-heavy-to-ligand-atom cutoff.
Inducia also uses that distance for maximum-count ranking, whereas GNINA's
internal ranking includes receptor side-chain hydrogens. Ties are resolved by
residue identity, making selection stable across hydrogen repair.
See [GNINA's selection implementation](https://github.com/gnina/gnina/blob/master/gninasrc/lib/flexinfo.cpp)
(checked 2026-09-09). The existing conservative rejection of blank-chain and
insertion-code flexible identifiers remains in place.

Both modes then run one mandatory GRAND child for every gnina survivor, using the paired gnina
complex and normalized ligand-pose SDF. GRAND survivors with finite CNNscores are ordered from
highest to lowest with pose ID as the deterministic tie-break. The first
`openbpmd.max_input_complexes` are submitted as OpenBPMD children, or every valid survivor when the
limit is `null` or larger than the survivor count. Finally, finite CompScores are ordered from
lowest to highest with pose ID as the tie-break.

Child jobs live below the parent `stages/` directory and their IDs are recorded by phase in parent
metadata. Local phases use a worker pool no larger than `max_local_workers`. Slurm phases submit all
siblings before waiting for any result. Backend, registry, resource, environment, and job-name
settings propagate to every child, and child results are read through canonical CADDForge result
records.

The parent writes timestamped phase, submission, status, heartbeat, child-progress, completion, and
failure summaries to `stdout.log`. Follow them while a job runs with
`caddforge logs <job-id> --follow`. Detailed command output remains in each child job's own
`stdout.log` and `stderr.log`; the parent includes those paths and concise progress summaries.
Structured parent events are published as `results/inducia_progress.jsonl`.

For allocated Slurm runs, parent events are written to node-local scratch first and synchronized to
the shared run directory periodically. Child progress snapshots use the same scratch-first,
best-effort publication pattern. Temporary registry or shared-filesystem failures are treated as
monitoring outages and retried; they do not change a scientific job's outcome. The parent result's
`extra.progress_log_sync` value records whether final log publication was complete, partial, or
unavailable at shutdown.

### Configuration options

Scientific settings are grouped in the strict nested mappings `phesa`, `clustering`,
`clash_detection`, `gnina`, `grand`, and `openbpmd`. Operational reporting settings use the strict
`progress` mapping. Unknown keys are rejected at every level.

#### Runtime and execution

- `run_directory`: optional parent run directory; otherwise CADDForge creates its standard dated
  job directory.
- `registry_path`: optional job-registry SQLite path shared by parent and child jobs.
- `backend`: `local` or `slurm`.
- `job_name`: optional readable parent name and child-name prefix.
- `export_results`: enable the core result-export behavior.
- `max_local_workers` (default `1`): maximum concurrent children in a local phase.
- `progress.interval_seconds` (default `30`): polling and best-effort shared-log synchronization
  interval in seconds; accepted range is 1 through 3600.
- `slurm_partition`, `slurm_time`, and `slurm_account`: optional Slurm scheduling values.
- `sbatch_options`: additional `sbatch` arguments.
- `slurm_nodes`, `slurm_ntasks`, and `slurm_cpus_per_task`: Slurm resource counts.
- `slurm_setup_commands`: commands used to prepare a shared Slurm environment.
- `slurm_env_pack_path`: optional packed environment; mutually exclusive with the wheelhouse mode.
- `slurm_env_cache_enabled` and `slurm_env_cache_dir`: control compute-node environment caching;
  a configured cache directory must be absolute on the compute node.
- `slurm_wheelhouse_path` and `slurm_requirements_path`: paired offline environment inputs; both
  must be set together and are mutually exclusive with `slurm_env_pack_path`.

#### General protocol and publication

- `residue_states` (default `{}`): explicit receptor Amber states keyed by `chain:number[:insertion]`.
  Otherwise named variants and diagnostic hydrogens define the states. Unresolved histidine
  tautomers must be prepared explicitly before docking. GNINA outputs retain the decisions in
  JSON records; GRAND repairs missing hydrogens with those pinned states, and representative
  exports carry mapped records into OpenBPMD.
- `ligand_resname` (default `LIG`): one-to-three-character ligand residue name passed to GRAND and
  OpenBPMD.
- `output_dir_name` (default `output`): publication directory below `results/`.
- `ranking_table_filename` (default `ranking_table.csv`) and `pose_manifest_filename` (default
  `pose_manifest.csv`): final CSV names.
- `best_complex_filename`, `best_protein_filename`, and `best_ligand_filename`: stable filenames
  for the rank-one structures.

All output names must be plain filenames or directory names without path separators. The mode and
input fields `mode`, `receptor_pdb_path`, `query_ligand_sdf_path`, `template_ligand_sdf_path`,
`autobox_ligand_path`, `box_center`, and `box_size` follow the rules in the preceding input section.

#### PheSA, clustering, and clash detection

- `phesa.max_confs` (default `200`): maximum conformers generated by each PheSA child. The 21
  protocol weights are fixed and cannot be overridden.
- `clustering.rmsd_cutoff_angstrom` (default `1.0`): direct-coordinate Butina clustering cutoff.
- `clash_detection.allowed_vdw_overlap_angstrom` (default `0.4`): tolerated heavy-atom van der
  Waals overlap before the protein residue becomes a flexibility candidate.

#### Gnina

- `gnina.gnina_executable` (default `gnina`): executable name or path.
- `gnina.auto_flex_radius_angstrom` (default `3.5`) and
  `gnina.auto_max_flexible_residues` (default `12`): auto-mode flexibility controls.
- `gnina.num_modes` (default `9`): auto docking output-mode count; template minimization
  always enforces one mode.
- `gnina.exhaustiveness`, `gnina.min_rmsd_filter`, and `gnina.scoring`: optional docking search and
  classical scoring settings.
- `gnina.cnn_scoring` (default `rescore`), `gnina.cnn`, and `gnina.cnn_model`: CNN scoring mode,
  built-in model name, and optional model file. CNN rescoring cannot be disabled because CNNscore is
  required for the OpenBPMD shortlist.
- `gnina.cpu`, `gnina.seed`, and `gnina.add_hydrogens`: optional CPU count, random seed, and gnina
  hydrogen-addition mode.
- `gnina.device` and `gnina.no_gpu`: optional GPU device or CPU-only flag; they are mutually
  exclusive.
- `gnina.pose_sort_order` (default `CNNscore`): gnina output ordering. Inducia independently uses
  the parsed pose-quality CNNscore for its shortlist.
- `gnina.autobox_add` (default `4.0`) and `gnina.autobox_extend`: autobox buffer controls.
- `gnina.extra_gnina_args`: additional non-protocol arguments. Inducia rejects attempts to override
  owned flags such as minimization, flexible residues, ligand/receptor paths, mode count, and CNN
  scoring.

#### GRAND

- `grand.sphere_radius_angstrom` (default `10.0`) and `grand.n_ghost_waters` (default `150`): GCMC
  sphere and ghost-water count.
- `grand.nonbonded_cutoff_nm` (default `0.9`), `grand.temperature_kelvin` (default `300.0`), and
  `grand.ph` (default `7.4`): nonbonded, temperature, and protonation metadata settings.
  GRAND preserves pinned states; changing this pH alone does not reassign protonation.
- `grand.initialization_passes` (default `7500`) and `grand.max_production_passes` (default
  `100000`): equilibration and production pass limits.
- `grand.padding_nm` (default `1.2`) and `grand.ionic_strength_molar` (default `0.15`): solvent-box
  padding and ionic strength.
- `grand.align_selection`: optional atom selection used when aligning sampled frames.
- `grand.frame_analysis_interval` (default `1`): analyze every Nth GRAND trajectory frame during
  representative selection.

#### OpenBPMD

- `openbpmd.max_input_complexes` (default `null`): maximum GRAND complexes submitted to OpenBPMD.
  A positive integer sets the limit; `null` passes every survivor with a finite CNNscore.
- `openbpmd.ligand_net_charge` (default `null`): optional assertion against the
  selected-frame ligand's formal charge; remove old hardcoded zeros for charged ligands.
- `openbpmd.solvent_padding_angstrom` (default `10.0`): padding for the fresh box.
  Inducia forwards `grand.ionic_strength_molar` as the new salt-pair target and
  always uses certified pocket input with hydrogen preservation.
- `openbpmd.openbpmd_script_path`: optional compatible OpenBPMD script; `null` uses the bundled
  script.
- `openbpmd.nreps` (default `10`), `openbpmd.hill_height` (default `0.05`), and
  `openbpmd.deposit_interval_ps` (default `1.0`): replica count and metadynamics bias settings.

### Failure semantics

Early PheSA and GNINA failures are published in `Result.extra.failure_summaries`
even when no final ranking was created. Existing stage CSV manifests now include
failure kind and available diagnostic paths; matching `.failures.json` files retain
every individual cause, pose/weight, child job ID, and manifest link. Parent errors
show at most three distinct causes with counts and bounded text. Missing canonical
child results are identified explicitly; available metadata and logs remain linked,
and absent evidence is not interpreted as a scientific cause. Failed GNINA producer
rows remain visible alongside accepted modes. These reports are available through
the optional `template_failures`, `gnina_failures`, and `downstream_failures` ports.

Input existence and absolute-path checks remain in configuration validation.
A file that disappears later is a distinct execution failure with its absolute
worker input path in the retained diagnostic record. Child scheduling still
submits the whole Slurm phase before waiting and isolates failed siblings.

Pose-level PheSA, clash-detection, gnina, GRAND, and OpenBPMD failures do not cancel successful
siblings. Every known pose remains in `pose_manifest.csv` with its stage, error, child job ID, and
available artifacts. A valid hydrated complex outside `openbpmd.max_input_complexes` is
`not_selected`, not failed. Missing, nonnumeric, NaN, or infinite CNNscore prevents OpenBPMD
submission for that pose; an invalid CompScore prevents final ranking for only that pose.

The parent job fails when no pose has a finite CompScore, when every mandatory upstream child fails,
or when publication/copy-back integrity fails. Partial stage and final manifests remain available
for diagnosis. A successful parent therefore always has at least one fully hydrated,
OpenBPMD-scored pose and all declared required output artifacts. Allocated runtimes publish through
scratch, copy file contents without preserving filesystem metadata for CIFS compatibility, verify
copied files by SHA256, and clean the scratch tree only after verified retention.
Failure handling preserves the original error and records secondary publication
or recovery errors in `extra.recovery_errors`. If retention remains incomplete,
`paths.retained_scratch_dir` identifies the preserved work directory. Publication
paths already relocated to durable storage can safely pass through recovery again.

### Outputs and pipeline artifacts

The plugin publishes every declared output port:

- `complex`: stable PDB copy of the best CompScore-ranked hydrated complex (`best_complex`).
- `pocket_complex`: the same best compact pocket structure, explicitly identified by `best_pocket_complex`.
- `pocket_manifest`: its hydration and identity manifest (`best_pocket_manifest`).
- `protein`: stable PDB copy of its protein (`best_protein`).
- `ligand`: stable SDF copy of its ligand (`best_ligand`).
- `ranking_table`: concise CSV with scored poses first in ascending CompScore order, followed by
  unselected and failed poses.
- `pose_manifest`: complete pose and child-job provenance CSV.
- `template_failures`: optional JSON view of PheSA and template-preparation failures.
- `gnina_failures`: optional JSON view of child failures and rejected GNINA poses.
- `downstream_failures`: optional JSON view of hydration and OpenBPMD failures.
- `complexes`: ordered path collection read from the manifest's `final_complex` column.
- `proteins`: ordered path collection read from `final_protein`.
- `ligands`: ordered path collection read from `final_ligand`.
- `output_dir`: persistent publication directory.

The canonical `results/result_record.json` also exposes logs, structured progress events, the
human-readable `results/inducia_summary.json`, `stages/`, and available template, gnina, and
downstream manifests.
Metrics include pose counts at every phase, the best gnina CNNscore, best CompScore, and best pose
ID. `pose_manifest.csv` records PheSA scores and clusters, flexible residues, all child IDs and
statuses, hydrated paths, CNNscore shortlist rank, OpenBPMD scores, final paths, and failures.

### Examples

- [Complete template-mode configuration](examples/inducia_template.yaml)
- [Complete auto-mode configuration](examples/inducia_auto.yaml)
- [Template-mode pipeline scatter over a query library](examples/inducia_pipeline_scatter.yaml)

All example paths are placeholders. Replace them with local inputs before running the commands
above. For auto autobox docking, replace `box_center` and `box_size` with `null` and set
`autobox_ligand_path` to one positioned 3D ligand.


### Remediation validation and migration

Scientific preparation policy 1 records the preservation, authoritative-ligand,
single-solvation and explicit-ion rules in preparation summaries. GRAND topology,
periodic archive, pocket certificate and runtime-coordinate formats each carry
an independent schema version. Summaries record installed Python/package and
AmberTools versions plus source hashes, so a source checkout can be distinguished
from the installed distribution. Input and artifact checksums tie the records to
the actual files.

For old GRAND archives, retain the original template used for their atom mapping.
For new GRAND runs, a chemistry-only template is sufficient. OpenBPMD inputs still
require an explicit-hydrogen ligand in the complex pose's coordinate frame.
Remove a manual charge override to derive formal charge automatically, or
keep an integer only as an assertion. Regenerate pocket certificates from a valid
raw trajectory/topology/ghost/cell archive; a bulk-water PDB alone cannot certify
sampled hydration. If that archive is missing, repeat hydration from a validated
pose. Poses with broken disulfides require regeneration from the intact receptor,
then hydration and scoring. Historical scores remain historical and are not
validated retroactively by these changes.

The generated scientific regression requires AmberTools, OpenFF, GRAND, OpenMM,
MDTraj, ParmEd, MDAnalysis and CUDA. Activate the scientific environment and set
`AMBERHOME`/`PATH` before Python starts, then run only the generated modules:

```bash
PYTHONPATH=src python -m pytest -q tests/test_openbpmd_real_preparation.py tests/test_grand_template_preparation.py
PYTHONPATH=src python -m pytest -q tests/test_inducia_generated_science.py
```

The second test supplies generated upstream docking records to both Inducia modes
and runs real GRAND sampling, Amber preparation, shortened OpenBPMD dynamics and
scoring, and parent publication, with local and scratch storage. It checks chemical
and artifact invariants as well as finite scores. It does not run PheSA/GNINA or
establish production sampling convergence, pose stability, or restored success
on the original receptor systems. Separate generated tests cover protein states,
disulfides, protected docking, ligand chemistry, large identifiers, and failure
retention. The default CI suite excludes external/slow scientific tests.

## Continuous integration and releases

`.github/workflows/release.yml` builds a wheel and source archive, validates
metadata, and tests the installed wheel on Linux with Python 3.10 and 3.12.
The lightweight test selection is recorded in `scripts/ci_tests.json`;
external scientific engines are validated separately. Python 3.10 tests use
the conditional `tomli` backport from the `[test]` extra. The wheel test runner
preserves the declared pytest markers while excluding checkout runtime sources.

Pushing a matching `v*` release tag publishes the tested artifacts to PyPI.
Branch pushes, pull requests and manual runs never publish. Configure a
PyPI Trusted Publisher for the `idorsia` repository, `release.yml` workflow
and `pypi` environment before pushing release tags. Runtime CADDForge
dependencies must already be available publicly; publish core first,
and publish docking/descriptors before MD.
