Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions SKILLS.md
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ mdfactory/
cli.py # CLI entry point (cyclopts App)
config.py # Reads config_templates/config.ini, sets env vars
build.py # build_mixedbox(), build_bilayer(), build_lnp(), build_proteinbox(), build_protein_mixedbox()
workflows.py # run_build_from_file() — dispatches to build.py
workflows.py # run_build_from_dict() — dispatches to build.py
prepare.py # df_to_build_input_models() — CSV row -> BuildInput
check.py # check_bilayer_buildable()
setup/protein.py # PDBFixer cleanup, pdb2gmx, GROMACS version gate (amber14sb needs >= 2026)
Expand Down Expand Up @@ -149,14 +149,14 @@ mdfactory clean [--molecules] [--database]
### Building Systems
```python
from mdfactory.prepare import df_to_build_input_models
from mdfactory.workflows import run_build_from_file
from mdfactory.workflows import run_build_from_dict
from mdfactory.build import build_mixedbox, build_bilayer

# CSV -> models
models, errors = df_to_build_input_models(df)

# YAML -> build
run_build_from_file(Path("hash.yaml"))
# BuildInput (or raw dict) -> build in the working directory
run_build_from_dict(models[0])
```

### Models
Expand Down
6 changes: 3 additions & 3 deletions docs/content/docs/developer-guide/adding-parametrization.mdx
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ and returns a `GromacsSingleMoleculeParameterSet`. It must handle three kinds of

- **small molecules**: generate (or reuse cached) parameters
- **water** (`species.is_water`): use `spec.water_model`
- **ions** (`species.is_ion`, Na⁺/Cl⁻; moleculetypes from `NATIVE_ION_MOLECULETYPES`)
- **ions** (`species.is_ion`, Na⁺/Cl⁻; moleculetypes from `NATIVE_SPECIES` in `mdfactory/models/species.py`)

`native_forcefield` is set for protein builds: a `NativeForceField(name, directory, water_model)`
pointing at the bundled protein force field. Water and ions must then come from that directory
Expand Down Expand Up @@ -163,7 +163,6 @@ class GaffParametrization(_FamilySpec):
simulation_types: ClassVar[frozenset[str]] = frozenset({"mixedbox", "bilayer", "lnp"})
default_protein: ClassVar[Optional[str]] = None
default_water_model: ClassVar[str] = "tip3p"
small_molecules_model: ClassVar[type[BaseModel]] = GaffSmallMolecules

family: Literal["gaff"] = "gaff"
small_molecules: Optional[GaffSmallMolecules] = None
Expand All @@ -180,7 +179,8 @@ class GaffParametrization(_FamilySpec):
... # optional: per-type restrictions on the resolved components
```

Then add it to `FAMILY_SPECS` and to the `ParametrizationSpec` union with a `Tag("gaff")`.
The default small-molecule spec is the model named in the `small_molecules` annotation. Add the
family to `FAMILY_SPECS` and to the `ParametrizationSpec` union with a `Tag("gaff")`.
If the family needs its own nonbonded settings, add a new `NonbondedScheme` value and its overlay
in `mdfactory/run_schedules/gromacs.py`: an entry in `NONBONDED` (electrostatics and van der
Waals keys) and one in `DISPCORR` (the dispersion correction per simulation type). Extend the
Expand Down
49 changes: 7 additions & 42 deletions docs/content/docs/user-guide/parametrization.mdx
Original file line number Diff line number Diff line change
Expand Up @@ -18,28 +18,15 @@ every simulation type; there is no global default.

| Family | Protein | Small molecules | Water | Ions | Nonbonded scheme | Simulation types |
|--------|---------|-----------------|-------|------|------------------|------------------|
| `charmm` | `charmm36m` | CGenFF | `tip3p` | CHARMM (`ions.itp`) | `charmm` (Force-switch) | all five |
| `amber` | `amber14sb` | OpenFF Sage | `tip3p` | Joung–Cheatham (`ions_tip3p.itp`) | `sage` (Potential-switch, `DispCorr EnerPres`) | `proteinbox`, `protein_mixedbox` |
| `charmm` | `charmm36m` | CGenFF | `tip3p` | CHARMM (`ions.itp`) | `charmm` | all five |
| `amber` | `amber14sb` | OpenFF Sage | `tip3p` | Joung–Cheatham (`ions_tip3p.itp`) | `sage` | `proteinbox`, `protein_mixedbox` |
| `openff` | n/a | OpenFF Sage | `opc3` | from the Sage `.offxml` | `sage` | `mixedbox`, `bilayer`, `lnp` |

The **nonbonded scheme** sets the nonbonded block of the generated run schedule, identical in all
four stages and for every simulation type:

- `charmm`: PME with `rcoulomb = 1.2`, `vdwtype = cutoff`, `vdw-modifier = Force-switch` from
`rvdw-switch = 1.0` to `rvdw = 1.2`, as all CHARMM36/CGenFF parameter files specify.
`DispCorr = EnerPres` for `mixedbox` and `lnp`, `no` for `bilayer` and the protein types
(CHARMM lipids and proteins were parametrized without a long-range correction).
- `sage`: PME with `rcoulomb = 0.9`, `vdwtype = cutoff`, `vdw-modifier = Potential-switch`
from `rvdw-switch = 0.8` to `rvdw = 0.9`, and `DispCorr = EnerPres`. These are the cutoff and
switch of the Sage `.offxml` files, which OpenFF treats as part of the force field.

See [Simulation Settings](/docs/user-guide/simulation-settings) for the rest of the protocol and
how to configure it.

<Callout type="info" title="Box size with the CHARMM cutoff">
The 1.2 nm CHARMM cutoff needs a box edge of at least about 2.6 nm (twice the pair-list
radius). The repository examples are far larger.
</Callout>
The **nonbonded scheme** sets the nonbonded block of the generated run schedule: `charmm` uses the
cutoffs and force switch that all CHARMM36/CGenFF parameter files specify, `sage` those of the
Sage `.offxml` files, which OpenFF treats as part of the force field. The values are listed once,
under [the generated protocol](/docs/user-guide/simulation-settings#the-generated-protocol) in
Simulation Settings, together with the coupling and the rest of the protocol.

## Syntax

Expand Down Expand Up @@ -248,28 +235,6 @@ Each parametrized species produces a `GromacsSingleMoleculeParameterSet`:
| `parameter_itp` | Path to the extra parameters ITP (`None` for native water/ions) |
| `forcefield_itp` | Path to the force field's `forcefield.itp` (CGenFF and native sets) |

## Migrating from the old keys

The flat `parametrization: cgenff | smirnoff | pdb2gmx` values and the `parametrization_config`
block were removed. They are rejected with a validation error:

| Old | New |
|-----|-----|
| `parametrization: cgenff` | `parametrization: charmm` |
| `parametrization: smirnoff` | `parametrization: openff` |
| `parametrization: pdb2gmx` | `parametrization: charmm` |
| `parametrization_config.forcefield` (protein) | `parametrization.protein` |
| `parametrization_config.water_model: opc3.offxml` | `parametrization.water_model: opc3` |
| `parametrization_config.forcefield` / `sage_forcefield` (OpenFF) | `parametrization.small_molecules.forcefield` |
| `parametrization_config.charge_method` | `parametrization.small_molecules.charge_method` |
| `parametrization_config.ignore_hydrogens` / `merge_all` | `system.protein.ignore_hydrogens` / `system.protein.merge_all` |

This is a clean break: build hashes change, and old build directories, run-database rows
(`input_data`) and pushed data cannot be parsed into the new `BuildInput`. Old OpenFF cache
entries under `parameter_store/gromacs/smirnoff/` are not reused; OpenFF parameters (including
NAGL/AM1-BCC charges) are regenerated once under `gromacs/openff/`. The CGenFF cache under
`gromacs/cgenff/` is unchanged and still valid.

## Troubleshooting

### OpenFF issues
Expand Down
7 changes: 4 additions & 3 deletions docs/content/docs/user-guide/troubleshooting.mdx
Original file line number Diff line number Diff line change
Expand Up @@ -165,11 +165,12 @@ the `--stages` filter — only files belonging to requested stages are deleted.
</Accordion>
<Accordion title="grompp: cut-off longer than half the box">

`charmm` builds use a 1.2 nm cutoff (Force-switch from 1.0 nm), so every box edge must be at
least about 2.6 nm. grompp reports *"The cut-off length is longer than half the shortest box
`charmm` builds use the longer nonbonded cutoff of the two schemes (see
[Simulation Settings](/docs/user-guide/simulation-settings#the-generated-protocol)), so every box
edge must be at least about 2.6 nm. grompp reports *"The cut-off length is longer than half the shortest box
vector or longer than the smallest box diagonal element"* for smaller boxes. Increase the
molecule count (mixedbox), the lipid count or `z_padding` (bilayer), or the padding (LNP,
protein types). `sage` builds (`openff`, `amber`) use 0.9 nm.
protein types). `sage` builds (`openff`, `amber`) have a shorter cutoff and rarely hit this.

</Accordion>
<Accordion title="Settings changes: duplicate inputs or refused re-renders">
Expand Down
9 changes: 4 additions & 5 deletions mdfactory/analysis/simulation.py
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@
}

# Register system_chemistry for all simulation types
from .utils import system_chemistry # noqa: E402
from .utils import read_build_metadata, system_chemistry # noqa: E402

for _sim_type, _registry in ANALYSIS_REGISTRY.items():
_registry["system_chemistry"] = system_chemistry
Expand Down Expand Up @@ -231,10 +231,9 @@ def metadata(self) -> dict[str, Any]:

"""
metadata = self.build_input.metadata
build_metadata_path = self.path / "build_metadata.json"
if build_metadata_path.is_file():
with open(build_metadata_path) as file_handle:
metadata["build_metadata"] = json.load(file_handle)
build_metadata = read_build_metadata(self.path)
if build_metadata is not None:
metadata["build_metadata"] = build_metadata
return metadata

@property
Expand Down
33 changes: 17 additions & 16 deletions mdfactory/analysis/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@

import pandas as pd

from mdfactory.models.input import BuildInput
from mdfactory.models.input import BuildInput, protein_anchor
from mdfactory.utils.utilities import load_yaml_file

from .constants import STATUS_ORDER
Expand Down Expand Up @@ -408,6 +408,14 @@ def get_chemistry_extractor(
raise ValueError(f"Unknown mode '{mode}'. Use 'all', 'lnp', or 'custom'.")


def read_build_metadata(directory: str | Path | None) -> dict[str, Any] | None:
"""Return the parsed ``build_metadata.json`` of a build directory, or None if absent."""
if directory is None:
return None
path = Path(directory) / "build_metadata.json"
return json.loads(path.read_text()) if path.is_file() else None


def system_chemistry(simulation, **kwargs) -> pd.DataFrame:
"""Extract species composition as a long-format DataFrame.

Expand Down Expand Up @@ -436,25 +444,18 @@ def system_chemistry(simulation, **kwargs) -> pd.DataFrame:

"""
build_input = simulation.build_input
achieved_species = []
build_metadata_path = getattr(simulation, "path", None)
if build_metadata_path is not None:
build_metadata_path = Path(build_metadata_path) / "build_metadata.json"
if build_metadata_path.is_file():
build_metadata = json.loads(build_metadata_path.read_text())
achieved_species = build_metadata.get("species", [])
build_metadata = read_build_metadata(getattr(simulation, "path", None))
achieved_species = (build_metadata or {}).get("species", [])
rows = []
if build_input.simulation_type == "protein_mixedbox":
# The protein anchors the box but is not part of system.species (packing
# is solution-only). List it as a species row here to match proteinbox
# output; it consumes no achieved_species index.
protein = build_input.system.protein
if (anchor := protein_anchor(build_input)) is not None:
# List the protein as a species row to match proteinbox output; it
# consumes no achieved_species index.
rows.append(
{
"resname": protein.resname,
"resname": anchor["resname"],
"smiles": None,
"count": protein.count,
"fraction": protein.fraction,
"count": anchor["count"],
"fraction": anchor["fraction"],
"simulation_type": build_input.simulation_type,
"concentration": None,
"resolved_count": None,
Expand Down
Loading
Loading