From 58335ed3f0dc604c2c223951cd4a674d20f1cfc4 Mon Sep 17 00:00:00 2001
From: Maximilian Scheurer <11406414+maxscheurer@users.noreply.github.com>
Date: Tue, 10 Nov 2026 18:02:26 +0100
Subject: [PATCH 1/2] chore: open PR for issue 68
From 20b15711c500d4d4c580b8646889ef6afc8082c6 Mon Sep 17 00:00:00 2001
From: Maximilian Scheurer <11406414+maxscheurer@users.noreply.github.com>
Date: Wed, 11 Nov 2026 13:30:11 +0100
Subject: [PATCH 2/2] Clean up legacy compat, dead code and duplication from
the protein/parametrization/mdp stack
Behaviour-neutral: hashes, generated mdps, topologies and cache layouts are unchanged.
- Remove Settings.DEFAULT_CONFIG, the docs migration section, identity(),
run_build_from_file, write-only GromacsProteinParameterSet fields and the
charmm36m-ljpme registry entry; make default-only parameters constants
- Require forcefield/chains/forcefield_types/out_path; report a missed LNP
density target through loguru
- Share topology walkers (iter_topology_entries, iter_includes), the native
species table, unit constants, the simulation-type definition and the
build_metadata reader
- Pass the minimum GROMACS version as a tuple; drop the small_molecules_model
ClassVar and the OpenFF protein override
- Factor the builder helpers (engine functions, solvation, bounding box,
protein topology staging) out of build.py
- Consolidate test fixtures and merge duplicated spec rejection/frozen tests
---
SKILLS.md | 8 +-
.../adding-parametrization.mdx | 6 +-
.../docs/user-guide/parametrization.mdx | 49 +--
.../docs/user-guide/troubleshooting.mdx | 7 +-
mdfactory/analysis/simulation.py | 9 +-
mdfactory/analysis/utils.py | 33 +-
mdfactory/build.py | 181 ++++++-----
mdfactory/models/composition.py | 5 +-
mdfactory/models/input.py | 29 +-
mdfactory/models/parametrization.py | 43 +--
mdfactory/models/simulation_settings.py | 4 -
mdfactory/models/species.py | 36 ++-
mdfactory/orchestration/apps.py | 39 +--
mdfactory/orchestration/simulate.py | 12 +-
mdfactory/orchestration/stages.py | 29 +-
mdfactory/parametrize.py | 114 +++----
mdfactory/prepare.py | 23 +-
mdfactory/settings.py | 33 --
mdfactory/setup/protein.py | 120 +++-----
mdfactory/simulation/openmm_utils.py | 77 ++---
mdfactory/tests/conftest.py | 63 +++-
mdfactory/tests/test_amber_gmx_build.py | 1 +
mdfactory/tests/test_csv_backend.py | 30 +-
mdfactory/tests/test_e2e_regressions.py | 18 +-
mdfactory/tests/test_models.py | 27 +-
.../tests/test_orchestration_gmx_version.py | 26 +-
mdfactory/tests/test_parametrization.py | 22 +-
.../tests/test_parametrization_schema.py | 64 +++-
mdfactory/tests/test_protein_mixedbox.py | 84 ++---
mdfactory/tests/test_proteinbox.py | 173 ++++++-----
mdfactory/tests/test_push.py | 35 +--
mdfactory/tests/test_push_analysis.py | 79 +----
mdfactory/tests/test_push_artifacts.py | 65 +---
mdfactory/tests/test_run_schedules.py | 9 +-
mdfactory/tests/test_simulation_settings.py | 11 +-
.../tests/test_sync_config_local_paths.py | 2 +-
mdfactory/tests/test_topology_utilities.py | 98 +++++-
mdfactory/utils/setup_utilities.py | 31 +-
mdfactory/utils/topology_utilities.py | 290 +++++++++---------
mdfactory/utils/units.py | 12 +
mdfactory/workflows.py | 19 --
41 files changed, 967 insertions(+), 1049 deletions(-)
create mode 100644 mdfactory/utils/units.py
diff --git a/SKILLS.md b/SKILLS.md
index 08e939b..8ef2820 100644
--- a/SKILLS.md
+++ b/SKILLS.md
@@ -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)
@@ -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
diff --git a/docs/content/docs/developer-guide/adding-parametrization.mdx b/docs/content/docs/developer-guide/adding-parametrization.mdx
index 8a50519..ad38d07 100644
--- a/docs/content/docs/developer-guide/adding-parametrization.mdx
+++ b/docs/content/docs/developer-guide/adding-parametrization.mdx
@@ -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
@@ -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
@@ -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
diff --git a/docs/content/docs/user-guide/parametrization.mdx b/docs/content/docs/user-guide/parametrization.mdx
index 9b78ac3..ddbe67e 100644
--- a/docs/content/docs/user-guide/parametrization.mdx
+++ b/docs/content/docs/user-guide/parametrization.mdx
@@ -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.
-
-
-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.
-
+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
@@ -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
diff --git a/docs/content/docs/user-guide/troubleshooting.mdx b/docs/content/docs/user-guide/troubleshooting.mdx
index 7ebe46a..ff96fa4 100644
--- a/docs/content/docs/user-guide/troubleshooting.mdx
+++ b/docs/content/docs/user-guide/troubleshooting.mdx
@@ -165,11 +165,12 @@ the `--stages` filter — only files belonging to requested stages are deleted.
-`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.
diff --git a/mdfactory/analysis/simulation.py b/mdfactory/analysis/simulation.py
index 4ac517e..7dcf987 100644
--- a/mdfactory/analysis/simulation.py
+++ b/mdfactory/analysis/simulation.py
@@ -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
@@ -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
diff --git a/mdfactory/analysis/utils.py b/mdfactory/analysis/utils.py
index 811850d..075b731 100644
--- a/mdfactory/analysis/utils.py
+++ b/mdfactory/analysis/utils.py
@@ -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
@@ -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.
@@ -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,
diff --git a/mdfactory/build.py b/mdfactory/build.py
index bd0d801..4f59c45 100644
--- a/mdfactory/build.py
+++ b/mdfactory/build.py
@@ -20,17 +20,59 @@
from .check import check_bilayer_buildable
from .models.composition import BilayerComposition, LNPComposition, MixedBoxComposition
from .models.input import BuildInput
-from .models.species import SingleMoleculeSpecies
+from .models.species import NATIVE_SPECIES, SingleMoleculeSpecies
from .run_schedules import write_mdp_files
from .setup.solvation import ionize, solvate
from .utils.setup_utilities import (
- N_AVOGADRO,
+ LOOSE_PACKING_FACTOR,
create_bilayer_from_model,
create_mixed_box_universe,
)
+from .utils.units import A3_TO_L, N_AVOGADRO
from .utils.utilities import working_directory
+def _native_species(smiles: str, count: int) -> SingleMoleculeSpecies:
+ """Return the water or ion species of *smiles* with its native residue name."""
+ return SingleMoleculeSpecies(count=count, smiles=smiles, resname=NATIVE_SPECIES[smiles].resname)
+
+
+def _resolve_engine_functions(inp: BuildInput):
+ """Return ``(parametrize_function, topology_function)`` for the build's family and engine."""
+ parametrize_function = resolve_parametrizer(inp.parametrization, inp.engine)
+ topology_function = DISPATCH_TOPOLOGY_BUILD.get(inp.engine, None)
+ if topology_function is None:
+ raise NotImplementedError(f"Cannot build topology for engine {inp.engine}.")
+ return parametrize_function, topology_function
+
+
+def _solvate_with_water(u: mda.Universe, **solvate_kwargs):
+ """Solvate *u* keeping its box; return the solvated universe and the water species.
+
+ Raises
+ ------
+ ValueError
+ If solvation added no water.
+ """
+ from loguru import logger
+
+ u_solvated = solvate(u, **solvate_kwargs)
+ u_solvated.dimensions = u.dimensions
+ n_water = len(u_solvated.select_atoms("water").residues)
+ logger.info(f"Solvation added {n_water} water molecules.")
+
+ if n_water == 0:
+ raise ValueError("Solvation did not add any water molecules.")
+ return u_solvated, _native_species("O", n_water)
+
+
+def _bounding_box(positions: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+ """Return the minimum corner, maximum corner and center of the bounding box of *positions*."""
+ lower = positions.min(axis=0)
+ upper = positions.max(axis=0)
+ return lower, upper, (lower + upper) / 2
+
+
def _compute_charge_from_universe(u: mda.Universe, species: list) -> int:
"""Compute total system charge from actual residue counts in the universe.
@@ -116,16 +158,13 @@ def create_mixedbox_by_compression(inp: BuildInput, require_target_density: bool
from loguru import logger
loose = inp.system.model_copy()
- loose.target_density = inp.system.target_density * 0.1
+ loose.target_density = inp.system.target_density * LOOSE_PACKING_FACTOR
u, _ = create_mixed_box_universe(loose)
logger.info("Successfully created universe for mixed box system.")
# parametrize all species in the box for the FF/engine combination
- parametrize_function = resolve_parametrizer(inp.parametrization, inp.engine)
- topology_function = DISPATCH_TOPOLOGY_BUILD.get(inp.engine, None)
- if topology_function is None:
- raise NotImplementedError(f"Cannot build topology for engine {inp.engine}.")
+ parametrize_function, topology_function = _resolve_engine_functions(inp)
parameters = [parametrize_function(spec) for spec in inp.system.species]
logger.info("Parametrization done.")
@@ -168,21 +207,7 @@ def build_bilayer(inp: BuildInput):
logger.info("Bilayer system equilibrated in vacuum.")
# 5. solvate
# TODO: solvation settings in model?
- u_solvated = solvate(u_equil, prune_in_z=True)
- u_solvated.dimensions = u_equil.dimensions
- n_water = len(u_solvated.select_atoms("water").residues)
-
- logger.info(f"Solvation added {n_water} water molecules.")
-
- if n_water == 0:
- raise ValueError("Solvation did not add any water molecules.")
-
- # TODO: refactor for more reusability?
- water_spec = SingleMoleculeSpecies(
- count=n_water,
- smiles="O",
- resname="SOL",
- )
+ u_solvated, water_spec = _solvate_with_water(u_equil, prune_in_z=True)
u_ionized, ion_species = ionize_solvated_system(
inp.system.ionization, u_solvated, inp.system.charge
@@ -223,16 +248,13 @@ def create_bilayer_by_compression(inp):
# 1. parametrize lipid species
# parametrize all species in the box for the FF/engine combination
- parametrize_function = resolve_parametrizer(inp.parametrization, inp.engine)
+ parametrize_function, topology_function = _resolve_engine_functions(inp)
# preliminary parameters (only lipids)
parameters = [parametrize_function(spec) for spec in inp.system.species]
logger.info("Parametrization done.")
# 2. write temporary toplogy file
- topology_function = DISPATCH_TOPOLOGY_BUILD.get(inp.engine, None)
- if topology_function is None:
- raise NotImplementedError(f"Cannot build topology for engine {inp.engine}.")
with working_directory("bilayer_squeeze", create=True, cleanup=True) as wd:
topology_function(u, inp.system.species, parameters, "bilayer")
@@ -299,9 +321,7 @@ def build_lnp(inp: BuildInput):
# Step 4: Solvate with spherical exclusion
# Calculate box size and center for solvation
- mins = u_lnp.atoms.positions.min(axis=0)
- maxs = u_lnp.atoms.positions.max(axis=0)
- center = (mins + maxs) / 2
+ mins, maxs, _ = _bounding_box(u_lnp.atoms.positions)
# Set box dimensions with padding
box_size = maxs - mins + 2 * inp.system.padding
@@ -314,17 +334,9 @@ def build_lnp(inp: BuildInput):
# Remove water inside the LNP core using remove_sphere
# Use a slightly smaller radius to avoid water clashing with shell
exclusion_radius = inp.system.core_radius * 0.9
- u_solvated = solvate(
+ u_solvated, water_spec = _solvate_with_water(
u_lnp, prune_in_z=False, remove_sphere=(center[0], center[1], center[2], exclusion_radius)
)
- u_solvated.dimensions = u_lnp.dimensions
- n_water = len(u_solvated.select_atoms("water").residues)
- logger.info(f"Solvation added {n_water} water molecules.")
-
- if n_water == 0:
- raise ValueError("Solvation did not add any water molecules.")
-
- water_spec = SingleMoleculeSpecies(count=n_water, smiles="O", resname="SOL")
# Step 5: Ionize
# Compute charge from actual residue counts (not input counts) since transformations
@@ -393,8 +405,8 @@ def create_lnp_core(inp: BuildInput):
)
# Use the existing mixedbox compression workflow for the unit cell. The lipid
- # core often cannot reach the nominal density within the compression budget;
- # keep the pre-existing behavior of continuing with a warning instead of failing.
+ # core often cannot reach the nominal density within the compression budget, so
+ # compression continues with a warning instead of failing.
u_unit, parameters, parametrize_function, topology_function = create_mixedbox_by_compression(
core_inp, require_target_density=False
)
@@ -665,8 +677,8 @@ def ionize_solvated_system(ion_config, u_solvated, total_charge, solvent_volume_
if n_water == 0:
raise ValueError("Cannot ionize system without water.")
if solvent_volume_a3 is not None:
- # 1 ų = 1e-27 L, so n_salt = concentration[mol/L] * V[L] * N_A.
- v_liters = solvent_volume_a3 * 1e-27
+ # n_salt = concentration[mol/L] * V[L] * N_A.
+ v_liters = solvent_volume_a3 * A3_TO_L
n_ions = int(round(c_ions * v_liters * N_AVOGADRO))
else:
M_water = 55.55 # mol/l
@@ -694,9 +706,9 @@ def ionize_solvated_system(ion_config, u_solvated, total_charge, solvent_volume_
additional_species = []
if num_na > 0:
- additional_species.append(SingleMoleculeSpecies(count=num_na, smiles="[Na+]", resname="NA"))
+ additional_species.append(_native_species("[Na+]", num_na))
if num_cl > 0:
- additional_species.append(SingleMoleculeSpecies(count=num_cl, smiles="[Cl-]", resname="CL"))
+ additional_species.append(_native_species("[Cl-]", num_cl))
return u_ionized, additional_species
@@ -708,7 +720,9 @@ def _ion_counts(ion_species) -> tuple[int, int]:
up by residue name; indexing by position would report a lone Cl- as Na+.
"""
counts = {spec.resname: spec.count for spec in ion_species}
- return counts.get("NA", 0), counts.get("CL", 0)
+ return counts.get(NATIVE_SPECIES["[Na+]"].resname, 0), counts.get(
+ NATIVE_SPECIES["[Cl-]"].resname, 0
+ )
def _center_in_cubic_box(u: mda.Universe, padding: float) -> float:
@@ -718,13 +732,10 @@ def _center_in_cubic_box(u: mda.Universe, padding: float) -> float:
geometry), so an asymmetric molecule keeps the requested clearance to each
face. Returns the cubic box edge length in Angstroms.
"""
- bbox_min = u.atoms.positions.min(axis=0)
- bbox_max = u.atoms.positions.max(axis=0)
- extent = bbox_max - bbox_min
- box_size = float(extent.max() + 2 * padding)
+ bbox_min, bbox_max, bbox_center = _bounding_box(u.atoms.positions)
+ box_size = float((bbox_max - bbox_min).max() + 2 * padding)
u.dimensions = [box_size, box_size, box_size, 90, 90, 90]
box_center = np.array([box_size / 2, box_size / 2, box_size / 2])
- bbox_center = (bbox_min + bbox_max) / 2
u.atoms.translate(box_center - bbox_center)
return box_size
@@ -768,6 +779,19 @@ def _prepare_protein(protein, spec):
return params
+def _copy_protein_topology(params, dest_dir: Path, topology_name: str) -> Path:
+ """Copy the pdb2gmx topology as ``dest_dir/topology_name``, with its chain and posre itps.
+
+ The topology references its per-chain and position-restraint itps by bare
+ name, so they must sit beside it. Returns the copied topology's path.
+ """
+ topology = dest_dir / topology_name
+ shutil.copy(params.topology_file, topology)
+ for include_file in params.topology_include_files:
+ shutil.copy(include_file, dest_dir / include_file.name)
+ return topology
+
+
def _relax_protein_system(u, topology, itps, ff_dir_name, protein_indices, **relax_kwargs):
"""Relax ``u`` with OpenMM in a scratch ``relaxation/`` dir, protein restrained.
@@ -844,13 +868,7 @@ def build_proteinbox(inp: BuildInput):
logger.info(f"Protein centered in cubic box: {box_size:.1f} A")
# Solvate (reuse existing solvation)
- u_solvated = solvate(u, prune_in_z=False)
- u_solvated.dimensions = u.dimensions
- n_water = len(u_solvated.select_atoms("water").residues)
- logger.info(f"Solvation added {n_water} water molecules.")
-
- if n_water == 0:
- raise ValueError("Solvation did not add any water molecules.")
+ u_solvated, _ = _solvate_with_water(u, prune_in_z=False)
# Ionize (reuse existing ionization)
u_ionized, ion_species = ionize_solvated_system(
@@ -861,10 +879,7 @@ def build_proteinbox(inp: BuildInput):
# Update topology [ molecules ] section. topol.top references its per-chain and
# posre .itp files by bare name, so they must sit alongside it.
- topology_dest = Path("topology.top")
- shutil.copy(params.topology_file, topology_dest)
- for include_file in params.topology_include_files:
- shutil.copy(include_file, Path(include_file.name))
+ topology_dest = _copy_protein_topology(params, Path(), "topology.top")
ff_dir_name = bundle_forcefield_into_topology(topology_dest)
update_topology_molecules(topology_dest, n_water_final, num_na, num_cl)
@@ -909,7 +924,7 @@ def _reimage_protein_centered(u: mda.Universe, n_protein_atoms: int, protein_ref
delta = current_displacement - reference_displacement
protein_whole = protein - box * np.round(delta / box)
- bbox_center = (protein_whole.min(axis=0) + protein_whole.max(axis=0)) / 2.0
+ bbox_center = _bounding_box(protein_whole)[2]
shift = box / 2.0 - bbox_center
pos[:n_protein_atoms] = protein_whole
pos += shift
@@ -920,10 +935,8 @@ def _reimage_protein_centered(u: mda.Universe, n_protein_atoms: int, protein_ref
solution.positions = solution.wrap(compound="residues")
-def _check_no_protein_clashes(
- u: mda.Universe, n_protein_atoms: int, cutoff: float = PROTEIN_MIXEDBOX_CLASH_CUTOFF_A
-) -> None:
- """Raise if any protein atom is closer than ``cutoff`` Å to a solution atom (periodic)."""
+def _check_no_protein_clashes(u: mda.Universe, n_protein_atoms: int) -> None:
+ """Raise if a protein atom is closer than the clash cutoff to a solution atom (periodic)."""
if n_protein_atoms >= u.atoms.n_atoms:
return
from MDAnalysis.lib.distances import capped_distance
@@ -931,18 +944,20 @@ def _check_no_protein_clashes(
protein_positions = u.atoms.positions[:n_protein_atoms]
solution_positions = u.atoms.positions[n_protein_atoms:]
pairs, _ = capped_distance(
- protein_positions, solution_positions, max_cutoff=cutoff, box=u.dimensions
+ protein_positions,
+ solution_positions,
+ max_cutoff=PROTEIN_MIXEDBOX_CLASH_CUTOFF_A,
+ box=u.dimensions,
)
if len(pairs) > 0:
raise ValueError(
- f"{len(pairs)} protein-solution atom pairs are closer than {cutoff} Å after "
- "re-imaging; the protein was likely mis-imaged across the PBC boundary."
+ f"{len(pairs)} protein-solution atom pairs are closer than "
+ f"{PROTEIN_MIXEDBOX_CLASH_CUTOFF_A} Å after re-imaging; the protein was "
+ "likely mis-imaged across the PBC boundary."
)
-def _check_no_solution_clashes(
- u: mda.Universe, n_protein_atoms: int, cutoff: float = PROTEIN_MIXEDBOX_CLASH_CUTOFF_A
-) -> None:
+def _check_no_solution_clashes(u: mda.Universe, n_protein_atoms: int) -> None:
"""Raise when atoms from different solution residues overlap under PBC."""
if n_protein_atoms >= u.atoms.n_atoms:
return
@@ -950,7 +965,10 @@ def _check_no_solution_clashes(
solution = u.atoms[n_protein_atoms:]
pairs, _ = capped_distance(
- solution.positions, solution.positions, max_cutoff=cutoff, box=u.dimensions
+ solution.positions,
+ solution.positions,
+ max_cutoff=PROTEIN_MIXEDBOX_CLASH_CUTOFF_A,
+ box=u.dimensions,
)
if len(pairs) == 0:
return
@@ -960,7 +978,7 @@ def _check_no_solution_clashes(
if n_clashes:
raise ValueError(
f"{n_clashes} solution-solution atom pairs from different residues are "
- f"closer than {cutoff} Å."
+ f"closer than {PROTEIN_MIXEDBOX_CLASH_CUTOFF_A} Å."
)
@@ -998,12 +1016,9 @@ def _stage_protein_topology(params, protein_ff_dir, dest_dir, species, parameter
"""
from .parametrize import generate_gromacs_topology_with_protein
- top_src = dest_dir / "topol_protein.top"
+ top_src = _copy_protein_topology(params, dest_dir, "topol_protein.top")
ff_dir = dest_dir / protein_ff_dir.name
- shutil.copy(params.topology_file, top_src)
shutil.copytree(protein_ff_dir, ff_dir, dirs_exist_ok=True)
- for include_file in params.topology_include_files:
- shutil.copy(include_file, dest_dir / include_file.name)
topology = dest_dir / "topology.top"
generate_gromacs_topology_with_protein(
top_src, ff_dir, species, parameters, counts, "protein_mixedbox", out_path=topology
@@ -1034,7 +1049,6 @@ def build_protein_mixedbox(inp: BuildInput):
from .setup.protein import bundle_forcefield_into_topology, extract_mass_from_topology
from .simulation.openmm_utils import compress_box
from .utils.setup_utilities import (
- PROTEIN_MIXEDBOX_DILUTION_FACTOR,
concentration_reference_volume_a3,
create_protein_mixedbox_universe,
cubic_box_edge_for_density,
@@ -1076,9 +1090,8 @@ def build_protein_mixedbox(inp: BuildInput):
# Reject a box that cannot hold the protein with the requested padding.
u_protein = mda.Universe(str(params.structure_file))
- extent = float(
- (u_protein.atoms.positions.max(axis=0) - u_protein.atoms.positions.min(axis=0)).max()
- )
+ lower, upper, _ = _bounding_box(u_protein.atoms.positions)
+ extent = float((upper - lower).max())
if box_edge < extent + 2 * system.padding:
raise ValueError(
f"Protein extent {extent:.1f} Å needs a box of at least "
@@ -1121,7 +1134,7 @@ def build_protein_mixedbox(inp: BuildInput):
# Pack the solution around the centered protein in a dilute cubic box (0.1×
# density idiom), then compress it down to the target box.
- dilute_edge = box_edge / PROTEIN_MIXEDBOX_DILUTION_FACTOR ** (1.0 / 3.0)
+ dilute_edge = box_edge / LOOSE_PACKING_FACTOR ** (1.0 / 3.0)
size_kw = (
{"target_volume": box_volume_a3}
if fixed_box
@@ -1235,7 +1248,7 @@ def build_protein_mixedbox(inp: BuildInput):
final_box_volume_a3 = float(u_relaxed.trajectory.ts.volume)
final_concentration_volume_a3 = reference_volume(final_box_volume_a3)
- final_concentration_volume_l = final_concentration_volume_a3 * 1e-27
+ final_concentration_volume_l = final_concentration_volume_a3 * A3_TO_L
build_report = {
"simulation_type": "protein_mixedbox",
"requested_box_volume_a3": box_volume_a3,
diff --git a/mdfactory/models/composition.py b/mdfactory/models/composition.py
index 88055f2..b46bb93 100644
--- a/mdfactory/models/composition.py
+++ b/mdfactory/models/composition.py
@@ -15,6 +15,7 @@
model_validator,
)
+from ..utils.units import A3_TO_CM3, N_AVOGADRO
from .species import (
LipidSpecies,
ProteinSpecies,
@@ -220,9 +221,9 @@ def calculate_counts(self, radius: float, min_count: int = 5) -> None:
"""
V_angstrom3 = (4 / 3) * np.pi * radius**3
- V_cm3 = V_angstrom3 * 1e-24 # ų to cm³
+ V_cm3 = V_angstrom3 * A3_TO_CM3
mass_g = self.target_density * V_cm3
- mass_dalton = mass_g * 6.022e23 # g to Daltons
+ mass_dalton = mass_g * N_AVOGADRO
# Weighted average mass
fractions = np.array([s.fraction for s in self.species])
diff --git a/mdfactory/models/input.py b/mdfactory/models/input.py
index cab0fe3..0045cee 100644
--- a/mdfactory/models/input.py
+++ b/mdfactory/models/input.py
@@ -18,7 +18,7 @@
ProteinBoxComposition,
ProteinMixedBoxComposition,
)
-from .parametrization import ParametrizationSpec
+from .parametrization import PROTEIN_SIMULATION_TYPES, ParametrizationSpec, SimulationType
from .simulation_settings import SimulationSettings
# Map system_type to the corresponding model class
@@ -32,13 +32,25 @@
}
+def protein_anchor(build_input: Any) -> dict[str, Any] | None:
+ """Return the species entry (resname, count, fraction) of a protein_mixedbox's protein.
+
+ The protein anchors a protein_mixedbox but is not part of ``system.species``,
+ which holds the packed solution only. Returns None for every other type.
+ """
+ if build_input.simulation_type != "protein_mixedbox":
+ return None
+ protein = build_input.system.protein
+ return {"resname": protein.resname, "count": protein.count, "fraction": protein.fraction}
+
+
class BuildInput(BaseModel):
"""Represent a complete simulation build specification with composition and parametrization."""
# Unknown keys are rejected so a misspelled or removed key is never dropped silently.
model_config = ConfigDict(extra="forbid")
- simulation_type: Literal["mixedbox", "bilayer", "lnp", "proteinbox", "protein_mixedbox"]
+ simulation_type: SimulationType
system: (
MixedBoxComposition
| BilayerComposition
@@ -130,16 +142,7 @@ def metadata(self) -> dict[str, Any]:
# The protein anchors the box but is not part of system.species
# (packing is solution-only); list it as a species so metadata
# matches proteinbox.
- protein = self.system.protein
- species_composition.insert(
- 0,
- {
- "resname": protein.resname,
- "count": protein.count,
- "fraction": protein.fraction,
- "concentration": None,
- },
- )
+ species_composition.insert(0, {**protein_anchor(self), "concentration": None})
return {
"hash": self.hash,
@@ -247,7 +250,7 @@ def validate_protein_chain_config(self) -> "BuildInput":
Compression additionally requires a multi-chain complex to be merged so
the barostat cannot scale chain centers independently.
"""
- if self.simulation_type not in ("proteinbox", "protein_mixedbox"):
+ if self.simulation_type not in PROTEIN_SIMULATION_TYPES:
return self
protein = self.system.protein
if (
diff --git a/mdfactory/models/parametrization.py b/mdfactory/models/parametrization.py
index 7d00c44..5598d99 100644
--- a/mdfactory/models/parametrization.py
+++ b/mdfactory/models/parametrization.py
@@ -9,7 +9,7 @@
import json
from pathlib import Path
-from typing import Annotated, Any, ClassVar, Literal, Optional, Self
+from typing import Annotated, Any, ClassVar, Literal, Optional, Self, get_args
from pydantic import (
AfterValidator,
@@ -30,9 +30,10 @@
#: potential switch from 0.8 nm. Adding a scheme means adding its nonbonded overlay.
NonbondedScheme = Literal["charmm", "sage"]
+ALL_SIMULATION_TYPES = frozenset(get_args(SimulationType))
PROTEIN_SIMULATION_TYPES = frozenset({"proteinbox", "protein_mixedbox"})
# Simulation types that parametrize small molecules (a proteinbox has none).
-SMALL_MOLECULE_SIMULATION_TYPES = frozenset({"mixedbox", "bilayer", "lnp", "protein_mixedbox"})
+SMALL_MOLECULE_SIMULATION_TYPES = ALL_SIMULATION_TYPES - {"proteinbox"}
# The CHARMM proteinbox pipeline solvates with 3-site spc216 coordinates and uses
# CHARMM Force-switch LJ in every mdp, so only these water models are supported.
@@ -131,13 +132,21 @@ class _FamilySpec(BaseModel):
simulation_types: ClassVar[frozenset[str]]
default_protein: ClassVar[Optional[str]]
default_water_model: ClassVar[str]
- small_molecules_model: ClassVar[type[BaseModel]]
-
family: str
protein: Optional[str] = None
small_molecules: Optional[BaseModel] = None
water_model: WaterModelName = None
+ @classmethod
+ def _default_small_molecules(cls) -> BaseModel:
+ """Return the family's default small-molecule spec (from the field's annotation)."""
+ (model,) = (
+ arg
+ for arg in get_args(cls.model_fields["small_molecules"].annotation)
+ if arg is not type(None)
+ )
+ return model()
+
@property
def nonbonded_scheme(self) -> NonbondedScheme:
"""Return the nonbonded scheme of the family's generated run schedule.
@@ -204,7 +213,7 @@ def resolve(self, simulation_type: str) -> Self:
"family": self.family,
"protein": (self.protein or self.default_protein) if has_protein else None,
"small_molecules": (
- (self.small_molecules or self.small_molecules_model()).model_dump()
+ (self.small_molecules or self._default_small_molecules()).model_dump()
if has_small_molecules
else None
),
@@ -227,12 +236,9 @@ class CharmmParametrization(_FamilySpec):
silently mismatch the mdp settings.
"""
- simulation_types: ClassVar[frozenset[str]] = frozenset(
- {"mixedbox", "bilayer", "lnp", "proteinbox", "protein_mixedbox"}
- )
+ simulation_types: ClassVar[frozenset[str]] = ALL_SIMULATION_TYPES
default_protein: ClassVar[Optional[str]] = "charmm36m"
default_water_model: ClassVar[str] = "tip3p"
- small_molecules_model: ClassVar[type[BaseModel]] = CgenffSmallMolecules
family: Literal["charmm"] = Field("charmm", description="Force-field family.")
protein: Optional[str] = Field(
@@ -309,7 +315,6 @@ class AmberParametrization(_FamilySpec):
simulation_types: ClassVar[frozenset[str]] = PROTEIN_SIMULATION_TYPES
default_protein: ClassVar[Optional[str]] = "amber14sb"
default_water_model: ClassVar[str] = "tip3p"
- small_molecules_model: ClassVar[type[BaseModel]] = OpenFFSmallMolecules
family: Literal["amber"] = Field("amber", description="Force-field family.")
protein: Optional[str] = Field(
@@ -349,15 +354,11 @@ def _check_amber_components(self) -> Self:
class OpenFFParametrization(_FamilySpec):
"""OpenFF family: Sage small molecules with OpenFF water and ions (no protein)."""
- simulation_types: ClassVar[frozenset[str]] = frozenset({"mixedbox", "bilayer", "lnp"})
+ simulation_types: ClassVar[frozenset[str]] = ALL_SIMULATION_TYPES - PROTEIN_SIMULATION_TYPES
default_protein: ClassVar[Optional[str]] = None
default_water_model: ClassVar[str] = "opc3"
- small_molecules_model: ClassVar[type[BaseModel]] = OpenFFSmallMolecules
family: Literal["openff"] = Field("openff", description="Force-field family.")
- protein: Optional[str] = Field(
- None, description="Not supported: the openff family has no protein force field."
- )
small_molecules: Optional[OpenFFSmallMolecules] = Field(
None, description="Small-molecule force field (OpenFF Sage)."
)
@@ -375,15 +376,6 @@ def ions(self) -> Optional[str]:
# Ion parameters come from the small-molecule .offxml.
return self.small_molecules.forcefield if self.small_molecules is not None else None
- @model_validator(mode="after")
- def _check_no_protein(self) -> Self:
- if self.protein is not None:
- raise ValueError(
- "The openff family has no protein force field; use family 'amber' "
- "for protein systems."
- )
- return self
-
FAMILY_SPECS: dict[str, type[_FamilySpec]] = {
"charmm": CharmmParametrization,
@@ -466,7 +458,4 @@ class GromacsProteinParameterSet(BaseModel):
topology_file: Path
structure_file: Path
topology_include_files: list[Path]
- chains: list[str]
- forcefield: str
- water_model: str
total_charge: int
diff --git a/mdfactory/models/simulation_settings.py b/mdfactory/models/simulation_settings.py
index 7970769..7e5222f 100644
--- a/mdfactory/models/simulation_settings.py
+++ b/mdfactory/models/simulation_settings.py
@@ -301,10 +301,6 @@ def resolve(self, simulation_type: str) -> Self:
data["production"]["length"] = DEFAULT_PRODUCTION_LENGTH_NS[simulation_type]
return type(self).model_validate(data)
- def identity(self) -> dict[str, Any]:
- """Return the fields that enter the build hash."""
- return {name: getattr(self, name) for name in sorted(self.IDENTITY_FIELDS)}
-
def _as_mapping(settings: "SimulationSettings | Mapping[str, Any] | None") -> dict[str, Any]:
if settings is None:
diff --git a/mdfactory/models/species.py b/mdfactory/models/species.py
index ed0ba98..cc831d9 100644
--- a/mdfactory/models/species.py
+++ b/mdfactory/models/species.py
@@ -7,7 +7,7 @@
import warnings
from functools import cached_property
from pathlib import Path
-from typing import Literal, Optional
+from typing import Literal, NamedTuple, Optional
from pydantic import BaseModel, Field, model_validator
from rdkit import Chem
@@ -38,6 +38,24 @@ def charge(self) -> int:
raise NotImplementedError("Charge property not implemented for base Species class.")
+class NativeSpecies(NamedTuple):
+ """Naming of a water or ion that the force-field directory provides itself."""
+
+ #: Residue name, which is also the GROMACS moleculetype name.
+ resname: str
+ #: Atom names in packing order (the same in CHARMM36m and ``amber14sb.ff``).
+ atom_names: tuple[str, ...]
+
+
+#: Native water and monatomic ions, keyed by SMILES.
+NATIVE_SPECIES: dict[str, NativeSpecies] = {
+ "O": NativeSpecies("SOL", ("OW", "HW1", "HW2")),
+ "[Na+]": NativeSpecies("NA", ("NA",)),
+ "[Cl-]": NativeSpecies("CL", ("CL",)),
+}
+WATER_RESNAME = NATIVE_SPECIES["O"].resname
+
+
class SingleMoleculeSpecies(Species):
"""Represent a single-molecule species identified by its canonical SMILES string."""
@@ -59,12 +77,11 @@ def is_water(self):
@property
def is_ion(self):
"""Return True if the species is a monatomic ion."""
- # TODO: expose/central location?
- allowed_ions = {"[Na+]", "[Cl-]"}
- ret = (
- self.smiles in allowed_ions and self.rdkit_molecule.GetAtoms()[0].GetFormalCharge() != 0
+ return (
+ self.smiles in NATIVE_SPECIES
+ and not self.is_water
+ and self.rdkit_molecule.GetAtoms()[0].GetFormalCharge() != 0
)
- return ret
@cached_property
def rdkit_molecule(self) -> "Chem.rdchem.Mol":
@@ -180,12 +197,11 @@ def check_fraction_or_count(self) -> "SolutionSpecies":
)
if (self.count is None) == (self.concentration is None):
raise ValueError("SolutionSpecies requires exactly one of 'count' or 'concentration'.")
- native_resnames = {"O": "SOL", "[Na+]": "NA", "[Cl-]": "CL"}
- expected_resname = native_resnames.get(self.smiles)
- if expected_resname is not None and self.resname.upper() != expected_resname:
+ native = NATIVE_SPECIES.get(self.smiles)
+ if native is not None and self.resname.upper() != native.resname:
raise ValueError(
f"Native water/ion species {self.smiles!r} must use resname "
- f"'{expected_resname}', not '{self.resname}'."
+ f"'{native.resname}', not '{self.resname}'."
)
return self
diff --git a/mdfactory/orchestration/apps.py b/mdfactory/orchestration/apps.py
index 4936cfc..14438a6 100644
--- a/mdfactory/orchestration/apps.py
+++ b/mdfactory/orchestration/apps.py
@@ -389,7 +389,7 @@ def _get_gromacs_detect_script() -> str:
)
-def _build_gmx_version_check(min_gmx_version: str) -> str:
+def _build_gmx_version_check(min_gmx_version: tuple[int, int]) -> str:
"""Return a bash snippet that exits unless ``$GMX_BIN`` is new enough.
The check runs on the compute node against the exact binary the job will
@@ -397,24 +397,19 @@ def _build_gmx_version_check(min_gmx_version: str) -> str:
Parameters
----------
- min_gmx_version : str
- Minimum GROMACS version, e.g. ``"2026"`` or ``"2026.1"``.
+ min_gmx_version : tuple of int
+ Minimum GROMACS version as ``(year, minor)``, e.g. ``(2026, 0)``.
Returns
-------
str
Bash block with a leading blank line and trailing newline.
- Raises
- ------
- ValueError
- If *min_gmx_version* is not of the form ``[.]``.
-
"""
- major, _, minor = min_gmx_version.partition(".")
- if not major.isdigit() or (minor and not minor.isdigit()):
- raise ValueError(f"Invalid minimum GROMACS version: {min_gmx_version!r}")
- threshold = int(major) * 1000 + int(minor or 0)
+ from ..setup.protein import format_gmx_version
+
+ threshold = min_gmx_version[0] * 1000 + min_gmx_version[1]
+ min_gmx_version = format_gmx_version(min_gmx_version)
return f"""
# Require GROMACS >= {min_gmx_version} for the force field in the topology
GMX_VERSION=$("$GMX_BIN" --version 2>/dev/null \\
@@ -535,7 +530,7 @@ def _build_grompp_script(
ref_file: str = "",
cpt_file: str = "",
maxwarn: int = 1,
- min_gmx_version: str = "",
+ min_gmx_version: tuple[int, int] | None = None,
) -> str:
"""Build the bash script for ``gmx grompp``.
@@ -560,11 +555,11 @@ def _build_grompp_script(
Checkpoint file for velocities (-t flag).
maxwarn : int
Maximum number of warnings to allow.
- min_gmx_version : str, optional
- Minimum GROMACS version the topology's force field needs, e.g.
- ``"2026"``. When given, the script exits with an error before grompp
- if the detected binary is older. Empty (default) means no check, and
- the script is unchanged.
+ min_gmx_version : tuple of int, optional
+ Minimum GROMACS version ``(year, minor)`` the topology's force field
+ needs, e.g. ``(2026, 0)``. When given, the script exits with an error
+ before grompp if the detected binary is older. None (default) means no
+ check, and the script is unchanged.
Returns
-------
@@ -574,7 +569,7 @@ def _build_grompp_script(
"""
ref_flag = f"-r {ref_file}" if ref_file else ""
cpt_flag = f"-t {cpt_file}" if cpt_file else ""
- version_check = _build_gmx_version_check(min_gmx_version) if min_gmx_version else ""
+ version_check = _build_gmx_version_check(min_gmx_version) if min_gmx_version is not None else ""
return f"""
set -euo pipefail # Exit on error, undefined vars, pipe failures
@@ -723,7 +718,7 @@ def run_grompp(
ref_file: str = "",
cpt_file: str = "",
maxwarn: int = 1,
- min_gmx_version: str = "",
+ min_gmx_version: tuple[int, int] | None = None,
stdout=None,
stderr=None,
inputs=None,
@@ -748,8 +743,8 @@ def run_grompp(
Checkpoint file for velocities (-t flag).
maxwarn : int
Maximum number of warnings to allow.
- min_gmx_version : str, optional
- Minimum GROMACS version required by the topology; empty for none.
+ min_gmx_version : tuple of int, optional
+ Minimum GROMACS version ``(year, minor)`` required by the topology.
stdout : str, optional
File path for stdout.
stderr : str, optional
diff --git a/mdfactory/orchestration/simulate.py b/mdfactory/orchestration/simulate.py
index def5bbf..9446685 100644
--- a/mdfactory/orchestration/simulate.py
+++ b/mdfactory/orchestration/simulate.py
@@ -344,8 +344,9 @@ def clean_simulation_outputs(
def _log_dry_run_plan(work_plan: list[dict], config: "ExecutorConfig") -> list[dict]:
"""Log the dry-run work plan with resolved gmx commands per stage."""
+ from ..setup.protein import format_gmx_version, min_gmx_version_for_topology
from .apps import _build_grompp_script, _build_mdrun_script
- from .stages import _extract_resource_hints, min_gmx_version_for_sim_dir
+ from .stages import _extract_resource_hints
logger.info("=" * 60)
logger.info("DRY RUN - No jobs will be submitted")
@@ -366,9 +367,12 @@ def _log_dry_run_plan(work_plan: list[dict], config: "ExecutorConfig") -> list[d
work_dir = str(Path(item["sim_dir"]).resolve())
stage_config = get_stage_config_or_none(config, item["stages"][0])
hints = _extract_resource_hints(stage_config)
- min_gmx_version = min_gmx_version_for_sim_dir(item["sim_dir"])
- if min_gmx_version:
- logger.info(f" Requires GROMACS >= {min_gmx_version} (checked before grompp)")
+ min_gmx_version = min_gmx_version_for_topology(Path(item["sim_dir"]) / "topology.top")
+ if min_gmx_version is not None:
+ logger.info(
+ f" Requires GROMACS >= {format_gmx_version(min_gmx_version)} "
+ "(checked before grompp)"
+ )
for stage_name in item["stages"]:
spec = STAGE_BY_NAME[stage_name]
diff --git a/mdfactory/orchestration/stages.py b/mdfactory/orchestration/stages.py
index 4b44789..8b9a8bf 100644
--- a/mdfactory/orchestration/stages.py
+++ b/mdfactory/orchestration/stages.py
@@ -12,6 +12,8 @@
from pathlib import Path
from typing import TYPE_CHECKING, Any, NamedTuple
+from ..run_schedules.gromacs import STAGE_FILES
+
if TYPE_CHECKING:
from parsl import AppFuture
@@ -74,7 +76,7 @@ class StageSpec:
StageSpec(
name="EM",
deffnm="min",
- mdp_file="em.mdp",
+ mdp_file=STAGE_FILES["minimization"],
gro_in="system.pdb",
gro_out="min.gro",
tpr_file="min.tpr",
@@ -88,7 +90,7 @@ class StageSpec:
StageSpec(
name="NVT",
deffnm="nvt",
- mdp_file="nvt.mdp",
+ mdp_file=STAGE_FILES["nvt"],
gro_in="min.gro",
gro_out="nvt.gro",
tpr_file="nvt.tpr",
@@ -102,7 +104,7 @@ class StageSpec:
StageSpec(
name="NPT",
deffnm="npt",
- mdp_file="npt.mdp",
+ mdp_file=STAGE_FILES["npt"],
gro_in="nvt.gro",
gro_out="npt.gro",
tpr_file="npt.tpr",
@@ -116,7 +118,7 @@ class StageSpec:
StageSpec(
name="Production",
deffnm="prod",
- mdp_file="md.mdp",
+ mdp_file=STAGE_FILES["production"],
gro_in="npt.gro",
gro_out="",
tpr_file="prod.tpr",
@@ -188,19 +190,6 @@ def _extract_resource_hints(stage_config: Any) -> ResourceHints:
return ResourceHints(ntasks=ntasks, disable_gpu=disable_gpu, gmx_binary=gmx_binary)
-def min_gmx_version_for_sim_dir(sim_dir: Path) -> str:
- """Return the minimum GROMACS version a simulation directory's topology needs.
-
- Derived from the force field ``topology.top`` includes (e.g.
- ``amber14sb.ff/forcefield.itp`` needs 2026). Returns ``""`` when there is no
- requirement (CHARMM and small-molecule builds) or no topology.
- """
- from ..setup.protein import format_gmx_version, min_gmx_version_for_topology
-
- required = min_gmx_version_for_topology(Path(sim_dir) / "topology.top")
- return format_gmx_version(required) if required is not None else ""
-
-
def run_stage(
spec: StageSpec,
sim_dir: Path,
@@ -287,8 +276,10 @@ def run_stage(
grompp_kwargs["cpt_file"] = spec.prereq_cpt
# Force fields that need a newer GROMACS (amber14sb) make grompp check the
# node's binary first; other topologies keep the unchanged script.
- min_gmx_version = min_gmx_version_for_sim_dir(sim_dir)
- if min_gmx_version:
+ from ..setup.protein import min_gmx_version_for_topology
+
+ min_gmx_version = min_gmx_version_for_topology(Path(sim_dir) / "topology.top")
+ if min_gmx_version is not None:
grompp_kwargs["min_gmx_version"] = min_gmx_version
grompp_fut = grompp_app(**grompp_kwargs)
diff --git a/mdfactory/parametrize.py b/mdfactory/parametrize.py
index e5a22fb..02654a4 100644
--- a/mdfactory/parametrize.py
+++ b/mdfactory/parametrize.py
@@ -15,7 +15,7 @@
GromacsSingleMoleculeParameterSet,
OpenFFSmallMolecules,
)
-from .models.species import SingleMoleculeSpecies
+from .models.species import NATIVE_SPECIES, WATER_RESNAME, SingleMoleculeSpecies
from .settings import settings
from .utils.topology_utilities import (
check_atomtype_collisions,
@@ -23,6 +23,8 @@
collect_forcefield_atomtype_definitions,
count_contiguous_strings,
extract_reusable_parts_from_cgenff_gmx_top,
+ iter_includes,
+ iter_topology_entries,
merge_extra_parameter_itps,
read_itp_atoms,
rename_atomtypes,
@@ -98,9 +100,6 @@ def _build_smirnoff_parameter_itp(top_path: Path, atomtypes_path: Path, out_path
fb.write(atom_fb.read())
-NATIVE_ION_MOLECULETYPES = {"[Na+]": "NA", "[Cl-]": "CL"}
-
-
@dataclass(frozen=True)
class NativeForceField:
"""A GROMACS force-field directory that provides native water and ion topologies.
@@ -161,29 +160,8 @@ def parametrize_cgenff_gromacs(
forcefield_itp = native.directory / "forcefield.itp"
# check for water/ions and return standard parameters
- if species.is_water:
- return GromacsSingleMoleculeParameterSet(
- moleculetype="SOL",
- smiles=species.smiles,
- component="water",
- forcefield=native.name,
- water_model=native.water_model,
- topology_merge="native",
- itp=native.directory / f"{native.water_model}.itp",
- parameter_itp=None,
- forcefield_itp=forcefield_itp,
- )
- if species.is_ion:
- return GromacsSingleMoleculeParameterSet(
- moleculetype=NATIVE_ION_MOLECULETYPES[species.smiles],
- smiles=species.smiles,
- component="ion",
- forcefield=native.name,
- topology_merge="native",
- itp=native.directory / "ions.itp",
- parameter_itp=None,
- forcefield_itp=forcefield_itp,
- )
+ if species.is_water or species.is_ion:
+ return _native_parameter_set(species, native, water_specific_ions=False)
# parametrize molecule with cgenff
pardir = settings.parameter_store / "gromacs" / "cgenff"
@@ -313,7 +291,7 @@ def atomtype_prefix(self) -> str:
def is_complete(self) -> bool:
"""Return True if every artifact exists and the completion marker was written."""
- artifacts = (self.itp, self.atomtypes, self.params, self.top, self.provenance)
+ artifacts = (self.itp, self.atomtypes, self.params, self.top)
return self.marker.is_file() and all(path.is_file() for path in artifacts)
@@ -344,12 +322,12 @@ def openff_water_cache_entry(species: SingleMoleculeSpecies, water_model: str) -
"""Return the cache entry of OpenFF water, keyed on the water model."""
key = {"format": OPENFF_CACHE_FORMAT_VERSION, "component": "water", "water_model": water_model}
workdir = _openff_root() / "water" / f"{_readable_stem(water_model)}_{_key_digest(key)}"
- return OpenFFCacheEntry(workdir.resolve(), "SOL", key, species.hash)
+ return OpenFFCacheEntry(workdir.resolve(), WATER_RESNAME, key, species.hash)
def openff_ion_cache_entry(species: SingleMoleculeSpecies, forcefield: str) -> OpenFFCacheEntry:
"""Return the cache entry of an OpenFF ion, keyed on the force field."""
- ion = NATIVE_ION_MOLECULETYPES[species.smiles]
+ ion = NATIVE_SPECIES[species.smiles].resname
key = {"format": OPENFF_CACHE_FORMAT_VERSION, "component": "ion", "forcefield": forcefield}
keydir = f"{_readable_stem(forcefield)}_{_key_digest(key)}"
workdir = _openff_root() / "ions" / keydir / ion
@@ -428,26 +406,27 @@ def _ensure_openff_entry(
entry.marker.touch()
-def _native_openff_parameter_set(
- species: SingleMoleculeSpecies, native: NativeForceField
+def _native_parameter_set(
+ species: SingleMoleculeSpecies, native: NativeForceField, *, water_specific_ions: bool
) -> GromacsSingleMoleculeParameterSet:
- """Return water or ion parameters taken from a GROMACS-native force-field directory."""
- if species.is_water:
- moleculetype = "SOL"
- component = "water"
- itp = native.directory / f"{native.water_model}.itp"
+ """Return water or ion parameters taken from a GROMACS-native force-field directory.
+
+ ``water_specific_ions`` selects ``ions_.itp`` (Amber's Joung-Cheatham
+ ions, recorded with their water model) over the CHARMM ``ions.itp``.
+ """
+ is_water = species.is_water
+ if is_water:
+ itp_name = f"{native.water_model}.itp"
else:
- moleculetype = NATIVE_ION_MOLECULETYPES[species.smiles]
- component = "ion"
- itp = native.directory / f"ions_{native.water_model}.itp"
+ itp_name = f"ions_{native.water_model}.itp" if water_specific_ions else "ions.itp"
return GromacsSingleMoleculeParameterSet(
- moleculetype=moleculetype,
+ moleculetype=NATIVE_SPECIES[species.smiles].resname,
smiles=species.smiles,
- component=component,
+ component="water" if is_water else "ion",
forcefield=native.name,
- water_model=native.water_model,
+ water_model=native.water_model if is_water or water_specific_ions else None,
topology_merge="native",
- itp=itp,
+ itp=native.directory / itp_name,
parameter_itp=None,
forcefield_itp=native.directory / "forcefield.itp",
)
@@ -483,7 +462,7 @@ def parametrize_openff_gromacs(
"""
if native_forcefield is not None and (species.is_water or species.is_ion):
- return _native_openff_parameter_set(species, native_forcefield)
+ return _native_parameter_set(species, native_forcefield, water_specific_ions=True)
# Lazy imports to avoid loading OpenFF when not needed
import numpy as np # noqa: PLC0415
@@ -504,13 +483,13 @@ def to_interchange(molecule, forcefield_name: str, **kwargs):
if species.is_water:
entry = openff_water_cache_entry(species, spec.water_model)
forcefield_name = f"{spec.water_model}.offxml"
- moleculetype = "SOL"
+ moleculetype = WATER_RESNAME
def build():
molecule = species.openff_molecule
# Name the molecule SOL regardless of the user's resname (e.g. WAT),
# so the generated ITP is always SOL_SOL.itp.
- molecule.name = "SOL"
+ molecule.name = WATER_RESNAME
return to_interchange(molecule, forcefield_name), molecule.name
_ensure_openff_entry(entry, species, moleculetype, build)
@@ -712,26 +691,19 @@ def _read_include_layout(
"""
ff_include_idx = None
water_include_idx = None
- system_idx = None
- molecules_idx = None
- for i, line in enumerate(lines):
- stripped = line.strip()
- if stripped.startswith("#include") and '"' in stripped:
- include_path = stripped.split('"')[1]
- if ".ff/" not in include_path:
- continue
- if "forcefield.itp" in include_path:
- if ff_include_idx is None:
- ff_include_idx = i
- elif ff_include_idx is not None and water_include_idx is None:
- water_include_idx = i
- elif stripped.startswith("[") and stripped.endswith("]"):
- section = stripped[1:-1].strip().lower()
- if section == "system" and system_idx is None:
- system_idx = i
- elif section == "molecules" and molecules_idx is None:
- molecules_idx = i
- return ff_include_idx, water_include_idx, system_idx, molecules_idx
+ for i, include_path in iter_includes(lines):
+ if ".ff/" not in include_path:
+ continue
+ if "forcefield.itp" in include_path:
+ if ff_include_idx is None:
+ ff_include_idx = i
+ elif ff_include_idx is not None and water_include_idx is None:
+ water_include_idx = i
+ header_idx: dict[str, int] = {}
+ for entry in iter_topology_entries(lines, headers=True):
+ if entry.is_header:
+ header_idx.setdefault(entry.section, entry.index)
+ return ff_include_idx, water_include_idx, header_idx.get("system"), header_idx.get("molecules")
def generate_gromacs_topology_with_protein(
@@ -741,7 +713,7 @@ def generate_gromacs_topology_with_protein(
parameters: list[GromacsSingleMoleculeParameterSet],
counts: list[int],
system_name: str,
- out_path: Path = Path("topology.top"),
+ out_path: Path,
) -> None:
"""Write one GROMACS topology for a pdb2gmx protein packed with solution molecules.
@@ -786,9 +758,9 @@ def generate_gromacs_topology_with_protein(
Actual molecule count for each species, aligned with ``species``.
system_name : str
Name written into the ``[ system ]`` section (inherited from pdb2gmx).
- out_path : Path, optional
+ out_path : Path
Destination topology path. Its parent receives the merged parameter file
- and copied small-molecule ITPs. Default ``topology.top`` in the cwd.
+ and copied small-molecule ITPs.
"""
out_path = Path(out_path)
@@ -831,7 +803,7 @@ def generate_gromacs_topology_with_protein(
check_atomtype_collisions(prm_files, forcefield_types)
merged = merge_extra_parameter_itps(prm_files, drop_defaults=True)
else:
- validate_cgenff_atomtype_compatibility(prm_files, ff_dir, forcefield_types)
+ validate_cgenff_atomtype_compatibility(prm_files, forcefield_types)
merged = merge_extra_parameter_itps(
prm_files, exclude_atomtypes=forcefield_types, drop_defaults=True
)
diff --git a/mdfactory/prepare.py b/mdfactory/prepare.py
index 3864792..c56044f 100644
--- a/mdfactory/prepare.py
+++ b/mdfactory/prepare.py
@@ -11,11 +11,9 @@
import pandas as pd
from .models.input import BuildInput
+from .models.parametrization import PROTEIN_SIMULATION_TYPES
from .models.simulation_settings import SimulationSettings, merge_simulation_settings
-# Simulation types whose system carries a protein with a pdb_path to resolve.
-PROTEIN_PDB_SIMULATION_TYPES = frozenset({"proteinbox", "protein_mixedbox"})
-
# CSV column holding the protein PDB path (flat, dotted notation).
PROTEIN_PDB_PATH_COLUMN = "system.protein.pdb_path"
@@ -212,6 +210,13 @@ def df_to_build_input_models(
return ret, errors
+def _resolve_relative(path: Any, base_dir: Path) -> Any:
+ """Return *path* made absolute against *base_dir*; blanks and absolute paths pass through."""
+ if pd.isna(path) or Path(path).is_absolute():
+ return path
+ return str((Path(base_dir) / path).resolve())
+
+
def resolve_protein_pdb_path(dct: Any, base_dir: Path) -> None:
"""Rewrite a relative protein ``pdb_path`` in a build dict to be absolute.
@@ -228,15 +233,13 @@ def resolve_protein_pdb_path(dct: Any, base_dir: Path) -> None:
Directory that relative paths are resolved against
"""
- if not isinstance(dct, dict) or dct.get("simulation_type") not in PROTEIN_PDB_SIMULATION_TYPES:
+ if not isinstance(dct, dict) or dct.get("simulation_type") not in PROTEIN_SIMULATION_TYPES:
return
protein = (dct.get("system") or {}).get("protein")
if not isinstance(protein, dict):
return
- pdb_path = protein.get("pdb_path")
- if pdb_path is None or Path(pdb_path).is_absolute():
- return
- protein["pdb_path"] = str((Path(base_dir) / pdb_path).resolve())
+ if protein.get("pdb_path") is not None:
+ protein["pdb_path"] = _resolve_relative(protein["pdb_path"], base_dir)
def load_build_input_csv(
@@ -272,9 +275,7 @@ def load_build_input_csv(
df = df.loc[:, ~df.columns.str.contains("^Unnamed")]
if PROTEIN_PDB_PATH_COLUMN in df.columns:
df[PROTEIN_PDB_PATH_COLUMN] = df[PROTEIN_PDB_PATH_COLUMN].apply(
- lambda p: (
- p if pd.isna(p) or Path(p).is_absolute() else str((csv_path.parent / p).resolve())
- )
+ _resolve_relative, base_dir=csv_path.parent
)
models, errors = df_to_build_input_models(df, simulation_settings)
return df, models, errors
diff --git a/mdfactory/settings.py b/mdfactory/settings.py
index ea45dd0..23f8476 100644
--- a/mdfactory/settings.py
+++ b/mdfactory/settings.py
@@ -82,39 +82,6 @@ def _get_defaults() -> dict[str, dict[str, str]]:
},
}
- # Keep DEFAULT_CONFIG as class attribute for backward compat in tests
- DEFAULT_CONFIG = {
- "cgenff": {"SILCSBIODIR": ""},
- "gromacs": {"GMX_PATH": "", "FORCEFIELD_DIR": ""},
- "storage": {"PARAMETERS": ""},
- "database": {"TYPE": "sqlite"},
- "databases": {
- "USE_RUN_DB": "true",
- "USE_ANALYSIS_DB": "true",
- },
- "sqlite": {
- "RUN_DB_PATH": "",
- "ANALYSIS_DB_PATH": "",
- },
- "csv": {
- "RUN_DB_PATH": "",
- "ANALYSIS_DB_PATH": "",
- },
- "foundry": {
- "BASE_PATH": "/Group Functions/mdfactory",
- "ANALYSIS_NAME": "analysis",
- "RUN_DB_PATH": "/Group Functions/mdfactory/runs",
- "ANALYSIS_DB_PATH": "/Group Functions/mdfactory/analysis",
- "ARTIFACT_DB_PATH": "/Group Functions/mdfactory/artifacts",
- },
- "slurm": {
- "ACCOUNT": "",
- "PARTITION_CPU": "",
- "PARTITION_GPU": "",
- "DEFAULT_QOS": "",
- },
- }
-
def __init__(self):
self.config = configparser.ConfigParser()
self._config_base_dir = Path.cwd()
diff --git a/mdfactory/setup/protein.py b/mdfactory/setup/protein.py
index d99490a..0af650a 100644
--- a/mdfactory/setup/protein.py
+++ b/mdfactory/setup/protein.py
@@ -17,7 +17,9 @@
from loguru import logger
from ..models.parametrization import GromacsProteinParameterSet
+from ..models.species import NATIVE_SPECIES
from ..settings import settings
+from ..utils.topology_utilities import iter_includes, iter_topology_entries
_DISULFIDE_CUTOFF_ANGSTROM = 2.2
_CYSTEINE_RESNAMES = {"CYS", "CYM", "CYX"}
@@ -28,11 +30,6 @@
"dirname": "charmm36-feb2026_cgenff-5.0.ff",
"sha256": "c2b0c70e21cced150528ee25d205e5bd3367a5e55e85c2ff0403911980528151",
},
- "charmm36m-ljpme": {
- "url": "https://mackerell.umaryland.edu/download.php?filename=CHARMM_ff_params_files/charmm36-feb2026_ljpme_cgenff-5.0.ff.tgz",
- "dirname": "charmm36-feb2026_ljpme_cgenff-5.0.ff",
- "sha256": "b2221dcba73066b6b5c30a52b9894216d03ea6c3e68206d9ce43857dd8116c21",
- },
}
@@ -54,10 +51,10 @@
_HISTIDINE_RESNAMES = {"HIS", "HID", "HIE", "HIP", "HSD", "HSE", "HSP"}
-def _resolve_protonation_residue_name(residue_name: str, forcefield: str | None) -> str:
+def _resolve_protonation_residue_name(residue_name: str, forcefield: str) -> str:
"""Return the residue name to write into the PDB for a protonation override."""
residue_name = residue_name.upper()
- if forcefield and forcefield.lower().startswith("charmm"):
+ if forcefield.lower().startswith("charmm"):
residue_name = _CHARMM_PROTONATION_ALIASES.get(residue_name, residue_name)
if len(residue_name) != 3:
@@ -192,7 +189,6 @@ def _resolve_cysteine_reference(
def _build_disulfide_prompt_input(
pdb_path: Path,
disulfide_bonds: list[tuple[str, str]],
- cutoff_angstrom: float = _DISULFIDE_CUTOFF_ANGSTROM,
) -> str:
"""Build deterministic yes/no input for pdb2gmx disulfide prompts."""
sg_atoms = _read_cysteine_sg_atoms(pdb_path)
@@ -210,7 +206,7 @@ def _build_disulfide_prompt_input(
candidate_pairs = [
(key_a, key_b)
for key_a, key_b in itertools.combinations(sg_atoms, 2)
- if math.dist(sg_atoms[key_a], sg_atoms[key_b]) <= cutoff_angstrom
+ if math.dist(sg_atoms[key_a], sg_atoms[key_b]) <= _DISULFIDE_CUTOFF_ANGSTROM
]
candidate_identities = {frozenset(pair) for pair in candidate_pairs}
@@ -219,7 +215,7 @@ def _build_disulfide_prompt_input(
formatted_missing = sorted(tuple(sorted(pair)) for pair in missing_pairs)
raise ValueError(
"Requested disulfide bonds were not detected as close CYS SG pairs "
- f"within {cutoff_angstrom:.1f} A in {pdb_path}: {formatted_missing}"
+ f"within {_DISULFIDE_CUTOFF_ANGSTROM:.1f} A in {pdb_path}: {formatted_missing}"
)
answers = ["y" if frozenset(pair) in requested else "n" for pair in candidate_pairs]
@@ -270,7 +266,7 @@ def _apply_protonation_states(
pdb_path: Path,
protonation_states: dict[str, str],
output_dir: Path,
- forcefield: str | None = None,
+ forcefield: str,
) -> Path:
"""Rename residues in the PDB so pdb2gmx applies the requested protonation states.
@@ -428,17 +424,13 @@ def min_gmx_version_for_topology(top_path: Path) -> tuple[int, int] | None:
if not top_path.is_file():
return None
requirements = []
- with open(top_path) as fb:
- for line in fb:
- stripped = line.strip()
- if not stripped.startswith("#include") or '"' not in stripped:
- continue
- parts = Path(stripped.split('"')[1]).parts
- if len(parts) < 2 or parts[-1] != "forcefield.itp" or not parts[-2].endswith(".ff"):
- continue
- required = min_gmx_version_for_forcefield(parts[-2])
- if required is not None:
- requirements.append(required)
+ for _, include in iter_includes(top_path.read_text().splitlines()):
+ parts = Path(include).parts
+ if len(parts) < 2 or parts[-1] != "forcefield.itp" or not parts[-2].endswith(".ff"):
+ continue
+ required = min_gmx_version_for_forcefield(parts[-2])
+ if required is not None:
+ requirements.append(required)
return max(requirements) if requirements else None
@@ -757,19 +749,11 @@ def _parse_protein_molecule_names(top_path: Path) -> list[str]:
pdb2gmx names each protein moleculetype ``Protein_chain_`` (or ``Protein``
when chains are merged), listing each once. Water and ion entries are ignored.
"""
- names: list[str] = []
- in_molecules = False
- with open(top_path) as f:
- for line in f:
- stripped = line.strip()
- if not stripped or stripped.startswith(";"):
- continue
- if stripped.startswith("["):
- in_molecules = stripped.strip("[] ").lower() == "molecules"
- continue
- if in_molecules and stripped.split()[0].startswith("Protein"):
- names.append(stripped.split()[0])
- return names
+ return [
+ entry.code.split()[0]
+ for entry in iter_topology_entries(Path(top_path).read_text().splitlines())
+ if entry.section == "molecules" and entry.code.startswith("Protein")
+ ]
def _run_gmx(
@@ -804,7 +788,7 @@ def run_pdb2gmx(
disulfide_bonds: list[tuple[str, str]],
protonation_states: dict[str, str],
output_dir: Path,
- chains: list[str] | None = None,
+ chains: list[str],
) -> GromacsProteinParameterSet:
"""Run gmx pdb2gmx to generate GROMACS topology from a protein PDB.
@@ -826,9 +810,9 @@ def run_pdb2gmx(
Residue-specific protonation state overrides.
output_dir : Path
Directory for output files.
- chains : list[str] | None
+ chains : list[str]
Declared chain identifiers. When given, the chains pdb2gmx produces must
- match exactly. When omitted, the PDB must yield a single protein chain.
+ match exactly. When empty, the PDB must yield a single protein chain.
Returns
-------
@@ -836,7 +820,6 @@ def run_pdb2gmx(
Paths to the generated topology, structure, and per-chain include files.
"""
- chains = chains or []
gmx_bin = str(check_gmx_available())
check_forcefield_available(forcefield)
resolved_ff = resolve_forcefield(forcefield)
@@ -946,9 +929,6 @@ def run_pdb2gmx(
topology_file=topology_file,
structure_file=structure_file,
topology_include_files=include_files,
- chains=produced_chains,
- forcefield=forcefield,
- water_model=water_model,
total_charge=total_charge,
)
@@ -960,24 +940,20 @@ def _sum_atoms_column(path: Path, column: int, follow_includes: bool = True) ->
level; force-field library includes (``.ff/``) are skipped.
"""
total = 0.0
- in_atoms = False
- for line in path.read_text().splitlines():
- stripped = line.strip()
- if stripped.startswith("#include"):
- include = stripped.split('"')[1] if '"' in stripped else None
- if follow_includes and include and ".ff/" not in include:
- itp_path = path.parent / include
- if itp_path.is_file():
- total += _sum_atoms_column(itp_path, column, follow_includes=False)
- elif stripped.startswith("["):
- in_atoms = stripped.strip("[] ").lower() == "atoms"
- elif in_atoms and stripped and not stripped.startswith(";"):
- parts = stripped.split()
- if len(parts) > column:
- try:
- total += float(parts[column])
- except ValueError:
- pass
+ lines = path.read_text().splitlines()
+ if follow_includes:
+ for _, include in iter_includes(lines):
+ if ".ff/" not in include and (path.parent / include).is_file():
+ total += _sum_atoms_column(path.parent / include, column, follow_includes=False)
+ for entry in iter_topology_entries(lines):
+ if entry.section != "atoms":
+ continue
+ parts = entry.code.split()
+ if len(parts) > column:
+ try:
+ total += float(parts[column])
+ except ValueError:
+ pass
return total
@@ -991,9 +967,7 @@ def extract_mass_from_topology(top_path: Path) -> float:
return _sum_atoms_column(top_path, 7)
-def update_topology_molecules(
- top_path: Path, n_water: int, num_na: int, num_cl: int, water_name: str = "SOL"
-) -> None:
+def update_topology_molecules(top_path: Path, n_water: int, num_na: int, num_cl: int) -> None:
"""Update the [ molecules ] section of a topology file after solvation/ionization.
Appends water and ion molecule counts to the topology.
@@ -1008,17 +982,13 @@ def update_topology_molecules(
Number of Na+ ions.
num_cl : int
Number of Cl- ions.
- water_name : str
- Residue name for water (default "SOL").
"""
+ entries = [("O", n_water), ("[Na+]", num_na), ("[Cl-]", num_cl)]
with open(top_path, "a") as f:
- if n_water > 0:
- f.write(f"{water_name:<20s} {n_water}\n")
- if num_na > 0:
- f.write(f"{'NA':<20s} {num_na}\n")
- if num_cl > 0:
- f.write(f"{'CL':<20s} {num_cl}\n")
+ for smiles, count in entries:
+ if count > 0:
+ f.write(f"{NATIVE_SPECIES[smiles].resname:<20s} {count}\n")
logger.info(f"Topology updated: {n_water} water, {num_na} Na+, {num_cl} Cl-")
@@ -1059,11 +1029,7 @@ def bundle_forcefield_into_topology(top_path: Path) -> str | None:
search_paths: list[Path] | None = None
rewrote_line = False
- for i, line in enumerate(lines):
- stripped = line.strip()
- if not stripped.startswith("#include") or '"' not in stripped:
- continue
- include_path = stripped.split('"')[1]
+ for i, include_path in iter_includes(lines):
if ".ff/" not in include_path:
continue
@@ -1100,7 +1066,7 @@ def bundle_forcefield_into_topology(top_path: Path) -> str | None:
# Relative includes already resolve from the topology's directory; only
# absolute includes need rewriting to the bundled relative form.
if is_absolute:
- lines[i] = line.replace(include_path, f"{ff_dir_name}/{remainder}")
+ lines[i] = lines[i].replace(include_path, f"{ff_dir_name}/{remainder}")
rewrote_line = True
if rewrote_line:
diff --git a/mdfactory/simulation/openmm_utils.py b/mdfactory/simulation/openmm_utils.py
index 3a0be45..4701760 100644
--- a/mdfactory/simulation/openmm_utils.py
+++ b/mdfactory/simulation/openmm_utils.py
@@ -9,10 +9,17 @@
import MDAnalysis as mda
import numpy as np
import openmm as mm
+from loguru import logger
from openmm import app, unit
+# Force constant (kJ/mol/nm^2) of the harmonic protein position restraint.
+RESTRAINT_K = 1000.0
+# Largest relative density overshoot (target_volume / volume - 1) that
+# finish_compression corrects by rescaling; larger overshoots raise.
+OVERSHOOT_TOLERANCE = 0.15
-def _add_protein_rmsd_restraint(system, positions, protein_indices, restraint_k):
+
+def _add_protein_rmsd_restraint(system, positions, protein_indices):
"""Restrain protein shape without opposing barostat translations.
``RMSDForce`` removes the best-fit rigid translation and rotation before
@@ -21,7 +28,7 @@ def _add_protein_rmsd_restraint(system, positions, protein_indices, restraint_k)
rmsd = mm.RMSDForce(positions, list(protein_indices))
restraint = mm.CustomCVForce("k*rmsd^2")
# rmsd^2 is a mean; scaling by the atom count gives k * sum(displacement**2).
- effective_k = restraint_k * len(protein_indices)
+ effective_k = RESTRAINT_K * len(protein_indices)
restraint.addGlobalParameter("k", effective_k * unit.kilojoules_per_mole / unit.nanometer**2)
restraint.addCollectiveVariable("rmsd", rmsd)
system.addForce(restraint)
@@ -41,15 +48,17 @@ def _state_data_reporter(target):
)
-def _universe_from_simulation(simulation, u, stem):
+def _universe_from_simulation(simulation, u, stem, periodic=True):
"""Write ``.pdb``/``.mmcif`` and return a copy of *u* with the final state.
- The returned universe keeps the topology of *u* and takes positions and
- box dimensions from the simulation's current state.
+ The returned universe keeps the topology of *u* and takes positions from the
+ simulation's current state. With *periodic* it also takes the box dimensions
+ and writes the box into the files; otherwise the universe has no box.
"""
state = simulation.context.getState(getPositions=True)
positions = state.getPositions()
- simulation.topology.setPeriodicBoxVectors(state.getPeriodicBoxVectors())
+ if periodic:
+ simulation.topology.setPeriodicBoxVectors(state.getPeriodicBoxVectors())
with open(f"{stem}.pdb", "w") as f:
app.PDBFile.writeFile(simulation.topology, positions, f, keepIds=True)
with open(f"{stem}.mmcif", "w") as f:
@@ -58,7 +67,7 @@ def _universe_from_simulation(simulation, u, stem):
# Use Merge instead of copy() to avoid file reference issues
ret = mda.Merge(u.atoms)
ret.atoms.positions = u_tmp.atoms.positions
- ret.dimensions = u_tmp.dimensions
+ ret.dimensions = u_tmp.dimensions if periodic else None
return ret
@@ -141,7 +150,6 @@ def finish_compression(
converged,
require_convergence,
not_reached_message,
- overshoot_tolerance=0.15,
rigid_indices=None,
):
"""Settle a compressed box exactly on its target volume.
@@ -158,38 +166,29 @@ def finish_compression(
Raise when not converged; otherwise warn and keep the box as is.
not_reached_message : str
Message used when the target was not reached.
- overshoot_tolerance : float, optional
- Largest relative density overshoot (``target_volume / volume - 1``)
- that is corrected by rescaling. Default is 0.15.
rigid_indices : iterable of int or None, optional
Particles moved as one rigid unit during rescaling.
- Returns
- -------
- float or None
- Relative overshoot that was corrected, or None when the target was not
- reached.
-
Raises
------
RuntimeError
If the target was not reached and *require_convergence* is True, or the
- overshoot exceeds *overshoot_tolerance*.
+ overshoot exceeds :data:`OVERSHOOT_TOLERANCE`.
"""
if not converged:
if require_convergence:
raise RuntimeError(not_reached_message)
- print(f"WARNING: {not_reached_message} Continuing with the box as compressed.")
- return None
+ logger.warning(f"{not_reached_message} Continuing with the box as compressed.")
+ return
box = context.getState().getPeriodicBoxVectors()
volume_nm3 = (box[0][0] * box[1][1] * box[2][2]).value_in_unit(unit.nanometer**3)
overshoot = target_volume_nm3 / volume_nm3 - 1.0
- if overshoot > overshoot_tolerance:
+ if overshoot > OVERSHOOT_TOLERANCE:
raise RuntimeError(
f"Box compression overshot the target density by {overshoot:.1%}, exceeding "
- f"the {overshoot_tolerance:.1%} that can be corrected by rescaling."
+ f"the {OVERSHOOT_TOLERANCE:.1%} that can be corrected by rescaling."
)
if volume_nm3 < target_volume_nm3:
scale = rescale_to_volume(context, target_volume_nm3, rigid_indices=rigid_indices)
@@ -197,7 +196,6 @@ def finish_compression(
f"Compression overshot the target density by {overshoot:.2%}; expanded the box "
f"from {volume_nm3:.3f} to {target_volume_nm3:.3f} nm³ (edge x{scale:.4f})."
)
- return max(overshoot, 0.0)
def monitor_volume_convergence(volumes, window_size=10, tolerance=0.02):
@@ -443,17 +441,7 @@ def create_sphere(u, pdb, top, radius, steps: int = 100_000):
simulation.step(steps)
- positions = simulation.context.getState(getPositions=True).getPositions()
- with open("sphere.pdb", "w") as f_sphere_pdb:
- app.PDBFile.writeFile(simulation.topology, positions, f_sphere_pdb, keepIds=True)
- with open("sphere.mmcif", "w") as f_sphere_mmcif:
- app.PDBxFile.writeFile(simulation.topology, positions, f_sphere_mmcif, keepIds=True)
- u_tmp = mda.Universe(app.PDBxFile("sphere.mmcif"))
- # Use Merge instead of copy() to avoid file reference issues
- ret = mda.Merge(u.atoms)
- ret.dimensions = None
- ret.atoms.positions = u_tmp.atoms.positions
- return ret
+ return _universe_from_simulation(simulation, u, "sphere", periodic=False)
def compress_box(
@@ -465,8 +453,6 @@ def compress_box(
steps_compression=200_000,
target_volume=None,
protein_indices=None,
- restraint_k=1000.0,
- overshoot_tolerance=0.15,
require_convergence=True,
):
"""Compress a periodic box to a target density using high-pressure NPT.
@@ -503,15 +489,10 @@ def compress_box(
protein_indices : list of int or None, optional
0-based atom indices to hold with harmonic position restraints. Default
is None (no restraints).
- restraint_k : float, optional
- Restraint force constant in kJ/mol/nm^2. Default is 1000.0.
- overshoot_tolerance : float, optional
- Largest relative density overshoot corrected by rescaling; larger
- overshoots raise. Default is 0.15 (15%).
require_convergence : bool, optional
Raise when the target is not reached within *steps_compression*. When
- False, only warn and return the box as compressed so far (the behavior
- LNP core compression relies on). Default is True.
+ False, only warn and return the box as compressed so far (used for the
+ LNP core, which cannot always reach its nominal density). Default is True.
Returns
-------
@@ -525,7 +506,7 @@ def compress_box(
RuntimeError
If the requested density or volume is not reached (and
*require_convergence* is True), or compression overshoots its target
- beyond *overshoot_tolerance*.
+ beyond :data:`OVERSHOOT_TOLERANCE`.
"""
# TODO: use proper logging
@@ -540,7 +521,7 @@ def compress_box(
)
if protein_indices is not None:
- _add_protein_rmsd_restraint(system, gro.positions, protein_indices, restraint_k)
+ _add_protein_rmsd_restraint(system, gro.positions, protein_indices)
barostat = mm.MonteCarloBarostat(
pressure * unit.bar,
@@ -599,7 +580,6 @@ def compress_box(
target_volume_nm3=target_volume_nm3,
converged=converged,
require_convergence=require_convergence,
- overshoot_tolerance=overshoot_tolerance,
rigid_indices=protein_indices,
not_reached_message=(
f"Box compression did not reach {target_volume_nm3:.3f} nm³ "
@@ -846,7 +826,6 @@ def relax_with_protein_restraints(
pdb,
top,
protein_indices,
- restraint_k=1000.0,
steps=10000,
pressure=1.0,
):
@@ -867,8 +846,6 @@ def relax_with_protein_restraints(
Path to the GROMACS topology file.
protein_indices : list of int
0-based atom indices to restrain (typically protein heavy atoms).
- restraint_k : float, optional
- Restraint force constant in kJ/mol/nm^2. Default is 1000.0.
steps : int, optional
Number of MD steps after minimization. Default is 10000.
pressure : float or None, optional
@@ -893,7 +870,7 @@ def relax_with_protein_restraints(
constraints=app.HBonds,
)
- _add_protein_rmsd_restraint(system, gro.positions, protein_indices, restraint_k)
+ _add_protein_rmsd_restraint(system, gro.positions, protein_indices)
if pressure is not None:
barostat = mm.MonteCarloBarostat(pressure * unit.bar, 300 * unit.kelvin)
diff --git a/mdfactory/tests/conftest.py b/mdfactory/tests/conftest.py
index c809c29..e458412 100644
--- a/mdfactory/tests/conftest.py
+++ b/mdfactory/tests/conftest.py
@@ -11,9 +11,58 @@
import pytest
import yaml
+from mdfactory.models.parametrization import PROTEIN_SIMULATION_TYPES
from mdfactory.prepare import df_to_build_input_models
from mdfactory.settings import Settings
+EXAMPLES_DIR = Path(__file__).resolve().parents[2] / "examples"
+
+PROTEIN_TYPES = sorted(PROTEIN_SIMULATION_TYPES)
+
+#: Smallest valid mixedbox build specification: 100 waters.
+WATER = {
+ "simulation_type": "mixedbox",
+ "parametrization": "charmm",
+ "system": {"species": [{"smiles": "O", "count": 100, "resname": "SOL"}]},
+}
+
+
+def make_forcefield_dir(parent: Path, name: str, files) -> Path:
+ """Create a stub GROMACS force-field directory ``parent/name``.
+
+ Parameters
+ ----------
+ parent : Path
+ Directory that receives the force field.
+ name : str
+ Directory name, e.g. ``charmm36.ff``.
+ files : mapping or iterable of str
+ File name to its text, or just file names (written as ``; ``).
+ """
+ ff_dir = Path(parent) / name
+ ff_dir.mkdir()
+ contents = files if isinstance(files, dict) else {f: f"; {f}\n" for f in files}
+ for filename, text in contents.items():
+ (ff_dir / filename).write_text(text)
+ return ff_dir
+
+
+def default_settings_loader(db_type: str):
+ """Return a ``Settings._load_config`` replacement: defaults only, on one backend.
+
+ Parameters
+ ----------
+ db_type : str
+ Value of ``[database] TYPE`` (``sqlite`` or ``csv``).
+ """
+
+ def load_config(self):
+ for section, options in Settings._get_defaults().items():
+ self.config[section] = dict(options)
+ self.config["database"]["TYPE"] = db_type
+
+ return load_config
+
def make_test_record(hash_value: str, status: str = "build", **overrides: Any) -> dict[str, str]:
"""Create a test database record for RUN_DATABASE.
@@ -135,13 +184,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return original_get_db_path(self, db_name, db_type)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
# Initialize the database via the public init path
from mdfactory.utils.push import init_systems_database
@@ -162,13 +206,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return original_get_db_path(self, db_name, db_type)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
# Initialize the database
from mdfactory.utils.push_analysis import init_analysis_database
diff --git a/mdfactory/tests/test_amber_gmx_build.py b/mdfactory/tests/test_amber_gmx_build.py
index 4135db3..fa2e1d8 100644
--- a/mdfactory/tests/test_amber_gmx_build.py
+++ b/mdfactory/tests/test_amber_gmx_build.py
@@ -128,6 +128,7 @@ def test_amber_sage_composite_topology_passes_grompp_build(tmp_path, monkeypatch
disulfide_bonds=[],
protonation_states={},
output_dir=tmp_path / "pdb2gmx",
+ chains=[],
)
ff_name = protein.bundle_forcefield_into_topology(params.topology_file)
assert ff_name == "amber14sb.ff"
diff --git a/mdfactory/tests/test_csv_backend.py b/mdfactory/tests/test_csv_backend.py
index f719c49..3dad0a3 100644
--- a/mdfactory/tests/test_csv_backend.py
+++ b/mdfactory/tests/test_csv_backend.py
@@ -11,7 +11,7 @@
from mdfactory.utils.data_manager import DataManager
from mdfactory.utils.db_operations import local_upload_with_modes
-from .conftest import make_test_record
+from .conftest import default_settings_loader, make_test_record
# ---------------------------------------------------------------------------
# Fixtures
@@ -30,13 +30,8 @@ def patched_get_csv_path(self, db_name: str) -> str:
return str(csv_path)
return original_get_csv_path(self, db_name)
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "csv"
-
monkeypatch.setattr(Settings, "get_csv_path", patched_get_csv_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("csv"))
from mdfactory.utils.push import init_systems_database
@@ -54,13 +49,8 @@ def patched_get_csv_path(self, db_name: str) -> str:
return str(analysis_dir / f"{db_name}.csv")
return Settings.get_csv_path(self, db_name)
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "csv"
-
monkeypatch.setattr(Settings, "get_csv_path", patched_get_csv_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("csv"))
from mdfactory.utils.push_analysis import init_analysis_database
@@ -115,13 +105,8 @@ def patched_get_csv_path(self, db_name: str) -> str:
return str(nested_path)
return str(tmp_path / f"{db_name}.csv")
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "csv"
-
monkeypatch.setattr(Settings, "get_csv_path", patched_get_csv_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("csv"))
from mdfactory.utils.push import init_systems_database
@@ -370,13 +355,8 @@ def test_csv_database_not_exists(tmp_path, monkeypatch):
def patched_get_csv_path(self, db_name: str) -> str:
return str(csv_path)
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "csv"
-
monkeypatch.setattr(Settings, "get_csv_path", patched_get_csv_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("csv"))
exists, path = DataManager.database_exists("RUN_DATABASE")
assert exists is False
diff --git a/mdfactory/tests/test_e2e_regressions.py b/mdfactory/tests/test_e2e_regressions.py
index b3841f6..9e0780f 100644
--- a/mdfactory/tests/test_e2e_regressions.py
+++ b/mdfactory/tests/test_e2e_regressions.py
@@ -9,6 +9,7 @@
import numpy as np
import openmm as mm
import pytest
+from loguru import logger
from openmm import unit
from mdfactory.models.species import SingleMoleculeSpecies
@@ -98,10 +99,9 @@ def test_refuses_to_shrink(self):
def _finish(context, target, **kwargs):
- """Call finish_compression with a converged, strict, 15%-tolerance default."""
+ """Call finish_compression with a converged, strict default."""
kwargs.setdefault("converged", True)
kwargs.setdefault("require_convergence", True)
- kwargs.setdefault("overshoot_tolerance", 0.15)
kwargs.setdefault("not_reached_message", "did not reach")
return finish_compression(context, target, **kwargs)
@@ -109,21 +109,25 @@ def _finish(context, target, **kwargs):
class TestFinishCompression:
def test_corrects_overshoot_within_tolerance(self):
context = _three_dimer_context() # 27 nm³, i.e. 4.5% below a 28.27 nm³ target
- overshoot = _finish(context, 28.27)
- assert overshoot == pytest.approx(28.27 / 27.0 - 1)
+ assert _finish(context, 28.27) is None
assert _volume(context) == pytest.approx(28.27)
def test_raises_beyond_tolerance(self):
with pytest.raises(RuntimeError, match="overshot the target density"):
_finish(_three_dimer_context(), 40.0)
- def test_unreached_target_raises_or_warns(self, capsys):
+ def test_unreached_target_raises_or_warns(self):
context = _three_dimer_context()
with pytest.raises(RuntimeError, match="did not reach"):
_finish(context, 20.0, converged=False)
- assert _finish(context, 20.0, converged=False, require_convergence=False) is None
+ messages = []
+ sink = logger.add(messages.append, level="WARNING")
+ try:
+ assert _finish(context, 20.0, converged=False, require_convergence=False) is None
+ finally:
+ logger.remove(sink)
assert _volume(context) == pytest.approx(27.0)
- assert "WARNING" in capsys.readouterr().out
+ assert any("did not reach" in str(m) for m in messages)
ETHANOL_ATOM_NAMES = ["C1", "C2", "O3", "H4", "H5", "H6", "H7", "H8", "H9"]
diff --git a/mdfactory/tests/test_models.py b/mdfactory/tests/test_models.py
index 6de1a22..b14aafa 100644
--- a/mdfactory/tests/test_models.py
+++ b/mdfactory/tests/test_models.py
@@ -19,7 +19,13 @@
distribute_counts,
)
from mdfactory.models.input import BuildInput
-from mdfactory.models.species import LipidSpecies, SingleMoleculeSpecies, Species
+from mdfactory.models.species import (
+ LipidSpecies,
+ SingleMoleculeSpecies,
+ SolutionSpecies,
+ Species,
+)
+from mdfactory.utils.setup_utilities import native_atom_names
from mdfactory.utils.utilities import working_directory
@@ -502,3 +508,22 @@ def test_shell_composition_counts_sum_correctly():
assert abs(ratios[0] - 0.33) < 0.05
assert abs(ratios[1] - 0.33) < 0.05
assert abs(ratios[2] - 0.34) < 0.05
+
+
+class TestNativeSpecies:
+ """Water and monatomic ions share one naming table."""
+
+ @pytest.mark.parametrize(
+ ("smiles", "resname", "atom_names"),
+ [("O", "SOL", 3), ("[Na+]", "NA", 1), ("[Cl-]", "CL", 1)],
+ )
+ def test_table_drives_validation_and_atom_names(self, smiles, resname, atom_names):
+ SolutionSpecies(smiles=smiles, resname=resname, count=1)
+ with pytest.raises(ValueError, match="must use resname"):
+ SolutionSpecies(smiles=smiles, resname="XXX", count=1)
+ assert len(native_atom_names(smiles)) == atom_names
+
+ def test_only_ions_are_ions(self):
+ assert SingleMoleculeSpecies(smiles="[Na+]", count=1, resname="NA").is_ion
+ assert not SingleMoleculeSpecies(smiles="O", count=1, resname="SOL").is_ion
+ assert not SingleMoleculeSpecies(smiles="CCO", count=1, resname="ETH").is_ion
diff --git a/mdfactory/tests/test_orchestration_gmx_version.py b/mdfactory/tests/test_orchestration_gmx_version.py
index 5d95c3e..88c85be 100644
--- a/mdfactory/tests/test_orchestration_gmx_version.py
+++ b/mdfactory/tests/test_orchestration_gmx_version.py
@@ -10,13 +10,13 @@
import pytest
from mdfactory.orchestration.apps import _build_grompp_script
-from mdfactory.orchestration.stages import STAGE_BY_NAME, min_gmx_version_for_sim_dir, run_stage
+from mdfactory.orchestration.stages import STAGE_BY_NAME, run_stage
AMBER_TOP = '#include "amber14sb.ff/forcefield.itp"\n#include "amber14sb.ff/tip3p.itp"\n'
CHARMM_TOP = '#include "charmm36-feb2026_cgenff-5.0.ff/forcefield.itp"\n'
-def _script(work_dir, min_gmx_version=""):
+def _script(work_dir, min_gmx_version=None):
return _build_grompp_script(
"em.mdp",
"system.pdb",
@@ -28,15 +28,7 @@ def _script(work_dir, min_gmx_version=""):
class TestRequirementDetection:
- @pytest.mark.parametrize(("top", "expected"), [(AMBER_TOP, "2026"), (CHARMM_TOP, "")])
- def test_from_topology(self, tmp_path, top, expected):
- (tmp_path / "topology.top").write_text(top)
- assert min_gmx_version_for_sim_dir(tmp_path) == expected
-
- def test_missing_topology_has_no_requirement(self, tmp_path):
- assert min_gmx_version_for_sim_dir(tmp_path) == ""
-
- @pytest.mark.parametrize(("top", "expected"), [(AMBER_TOP, "2026"), (CHARMM_TOP, None)])
+ @pytest.mark.parametrize(("top", "expected"), [(AMBER_TOP, (2026, 0)), (CHARMM_TOP, None)])
def test_run_stage_passes_requirement_only_when_needed(self, tmp_path, top, expected):
(tmp_path / "topology.top").write_text(top)
grompp_app = MagicMock(return_value=MagicMock())
@@ -53,16 +45,16 @@ def test_no_check_without_minimum(self, tmp_path):
)
def test_check_runs_after_detection_and_before_grompp(self, tmp_path):
- script = _script(tmp_path, "2026")
+ script = _script(tmp_path, (2026, 0))
assert "GROMACS >= 2026" in script
detect = script.index("GMX_BIN=gmx")
check = script.index("GMX_VERSION=")
assert detect < check < script.index("$GMX_BIN grompp")
- @pytest.mark.parametrize("bad", ["abc", "2026.x", "-1"])
- def test_rejects_invalid_minimum(self, tmp_path, bad):
- with pytest.raises(ValueError, match="Invalid minimum GROMACS version"):
- _script(tmp_path, bad)
+ def test_minor_version_is_part_of_the_threshold(self, tmp_path):
+ script = _script(tmp_path, (2026, 1))
+ assert "GROMACS >= 2026.1" in script
+ assert "GMX_VERSION_NUM < 2026001" in script
def _fake_gmx(bin_dir, name="gmx"):
@@ -109,7 +101,7 @@ def test_fake_gmx(self, tmp_path, version, passes):
_fake_gmx(bin_dir)
log = tmp_path / "grompp.log"
result = subprocess.run(
- ["bash", "-c", _script(sim_dir, "2026")],
+ ["bash", "-c", _script(sim_dir, (2026, 0))],
capture_output=True,
text=True,
env={
diff --git a/mdfactory/tests/test_parametrization.py b/mdfactory/tests/test_parametrization.py
index 315b5aa..22274b6 100644
--- a/mdfactory/tests/test_parametrization.py
+++ b/mdfactory/tests/test_parametrization.py
@@ -36,6 +36,7 @@
from mdfactory.utils.data_manager import DataManager, SQLiteDataSource
from mdfactory.utils.utilities import working_directory
+from .conftest import make_forcefield_dir
from .utils import ProcessWithException
CHARMM_MIXEDBOX = CharmmParametrization().resolve("mixedbox")
@@ -68,7 +69,7 @@ def is_cgenff_available():
@pytest.mark.skipif(not is_cgenff_available(), reason="CGenFF not available")
-def test_parametrization_simple(tmp_path, mocker):
+def test_parametrization_simple(tmp_path, mocker, parameter_store):
mock_data_source = SQLiteDataSource(
db_path=str(tmp_path / "molecule_database.db"), table_name="molecules"
)
@@ -76,8 +77,7 @@ def test_parametrization_simple(tmp_path, mocker):
patcher = patch.object(DataManager, "_initialize_data_source", return_value=mock_data_source)
patcher.start()
- config.parameter_store = tmp_path / "parameters"
- path = config.parameter_store / "gromacs" / "cgenff"
+ path = parameter_store / "gromacs" / "cgenff"
parametrize = partial(parametrize_cgenff_gromacs, spec=CHARMM_MIXEDBOX)
@@ -124,9 +124,8 @@ def test_parametrization_simple(tmp_path, mocker):
@pytest.mark.skipif(not is_cgenff_available(), reason="CGenFF not available")
-def test_parametrization_simultaneously(tmp_path):
- config.parameter_store = tmp_path / "parameters"
- path = config.parameter_store / "gromacs" / "cgenff"
+def test_parametrization_simultaneously(tmp_path, parameter_store):
+ path = parameter_store / "gromacs" / "cgenff"
parametrize = partial(parametrize_cgenff_gromacs, spec=CHARMM_MIXEDBOX)
@@ -288,9 +287,7 @@ def test_cgenff_locks_written_directory_and_regenerates_failed_attempt(tmp_path,
_fake_cgenff(monkeypatch, calls)
locked = _record_locks(monkeypatch)
# Use a temporary force-field directory so the test does not need a CGenFF installation.
- ff_dir = tmp_path / "charmm36.ff"
- ff_dir.mkdir()
- (ff_dir / "forcefield.itp").write_text("; forcefield.itp\n")
+ ff_dir = make_forcefield_dir(tmp_path, "charmm36.ff", ["forcefield.itp"])
native = parametrize_module.NativeForceField("charmm36", ff_dir, "tip3p")
spec = SingleMoleculeSpecies(smiles="CCO", count=1, resname="ETH")
# A failed earlier attempt left a directory without outputs.
@@ -315,10 +312,9 @@ def test_cgenff_locks_written_directory_and_regenerates_failed_attempt(tmp_path,
def test_cgenff_native_water_and_ions_record_forcefield(tmp_path, monkeypatch):
- ff_dir = tmp_path / "charmm36.ff"
- ff_dir.mkdir()
- for name in ("forcefield.itp", "tip3p.itp", "ions.itp"):
- (ff_dir / name).write_text(f"; {name}\n")
+ ff_dir = make_forcefield_dir(
+ tmp_path, "charmm36.ff", ["forcefield.itp", "tip3p.itp", "ions.itp"]
+ )
monkeypatch.setattr(config, "parameter_store", tmp_path / "parameters")
native = parametrize_module.NativeForceField("charmm36m", ff_dir, "tip3p")
diff --git a/mdfactory/tests/test_parametrization_schema.py b/mdfactory/tests/test_parametrization_schema.py
index 1ad43d0..3a78253 100644
--- a/mdfactory/tests/test_parametrization_schema.py
+++ b/mdfactory/tests/test_parametrization_schema.py
@@ -3,30 +3,31 @@
"""Tests for the per-component force-field family specs on BuildInput."""
import json
-from pathlib import Path
import pytest
import yaml
from pydantic import ValidationError
-from mdfactory.models.input import BuildInput
+from mdfactory.models.input import BuildInput, type_mapping
from mdfactory.models.parametrization import (
+ ALL_SIMULATION_TYPES,
DEFAULT_NAGL_MODEL,
DEFAULT_OPENFF_FORCEFIELD,
+ PROTEIN_SIMULATION_TYPES,
+ SMALL_MOLECULE_SIMULATION_TYPES,
AmberParametrization,
CgenffSmallMolecules,
CharmmParametrization,
OpenFFParametrization,
OpenFFSmallMolecules,
)
+from mdfactory.models.simulation_settings import DEFAULT_PRODUCTION_LENGTH_NS
from mdfactory.prepare import load_build_input_csv, resolve_protein_pdb_path
+from mdfactory.run_schedules.gromacs import COUPLING
-from .conftest import protein_input_data
-
-EXAMPLES_DIR = Path(__file__).resolve().parents[2] / "examples"
+from .conftest import EXAMPLES_DIR, PROTEIN_TYPES, protein_input_data
SMALL_MOLECULE_TYPES = ["mixedbox", "bilayer", "lnp"]
-PROTEIN_TYPES = ["proteinbox", "protein_mixedbox"]
VALID_COMBINATIONS = [
("charmm", "mixedbox"),
@@ -219,13 +220,43 @@ def test_charmm_proteinbox_accepts_other_charmm_forcefields_and_waters(protein_p
assert (spec.protein, spec.water_model) == ("charmm27", "spce")
-def test_resolved_spec_is_frozen(protein_pdb):
+def test_resolved_spec_keeps_cgenff_small_molecules(protein_pdb):
spec = _build("mixedbox", "charmm", protein_pdb).parametrization
- with pytest.raises(ValidationError):
- spec.water_model = "spce"
assert isinstance(spec.small_molecules, CgenffSmallMolecules)
+@pytest.mark.parametrize(
+ ("simulation_type", "family", "attribute", "value"),
+ [
+ ("mixedbox", "charmm", "water_model", "spce"),
+ *[(t, f, "protein", "other") for t in PROTEIN_TYPES for f in ("charmm", "amber")],
+ ],
+)
+def test_resolved_spec_is_frozen(simulation_type, family, attribute, value, protein_pdb):
+ spec = _build(simulation_type, family, protein_pdb).parametrization
+ with pytest.raises(ValidationError):
+ setattr(spec, attribute, value)
+
+
+# Unsupported protein force-field and water choices, rejected for both protein types.
+PROTEIN_CHOICE_REJECTIONS = [
+ (sim_type, {"family": family, **overrides}, match)
+ for sim_type in PROTEIN_TYPES
+ for family, overrides, match in [
+ ("charmm", {"protein": "amber99sb-ildn"}, "CHARMM"),
+ ("charmm", {"water_model": "tip4p"}, "3-site"),
+ ("charmm", {"protein": "charmm36m-ljpme"}, "LJ-PME"),
+ ("amber", {"protein": "amber19sb"}, "Only 'amber14sb'"),
+ ("amber", {"protein": "charmm36m"}, "Only 'amber14sb'"),
+ ("amber", {"water_model": "opc"}, "Only 'tip3p'"),
+ ("amber", {"water_model": "spce"}, "Only 'tip3p'"),
+ ]
+] + [
+ ("protein_mixedbox", {"family": "charmm", "protein": "charmm27"}, "only the registered"),
+ ("protein_mixedbox", {"family": "charmm", "water_model": "spce"}, "only 'tip3p'"),
+]
+
+
@pytest.mark.parametrize(
("simulation_type", "parametrization", "match"),
[
@@ -242,10 +273,9 @@ def test_resolved_spec_is_frozen(protein_pdb):
{"family": "charmm", "small_molecules": {"forcefield": "cgenff"}},
"small_molecules does not apply",
),
- ("bilayer", {"family": "openff", "protein": "x"}, "no protein force field"),
+ ("bilayer", {"family": "openff", "protein": "x"}, "protein does not apply"),
# CHARMM restrictions.
("proteinbox", {"family": "charmm", "protein": "amber14sb"}, "only CHARMM force fields"),
- ("proteinbox", {"family": "charmm", "protein": "charmm36m-ljpme"}, "LJ-PME"),
("proteinbox", {"family": "charmm", "water_model": "tip4p"}, "3-site"),
("mixedbox", {"family": "charmm", "water_model": "spce"}, "TIP3P water of its"),
("protein_mixedbox", {"family": "charmm", "protein": "charmm27"}, "'charmm36m' bundle"),
@@ -265,6 +295,7 @@ def test_resolved_spec_is_frozen(protein_pdb):
{"family": "openff", "small_molecules": {"forcefield": "openff-2.2.0"}},
"must be an .offxml file name",
),
+ *PROTEIN_CHOICE_REJECTIONS,
],
)
def test_rejections(simulation_type, parametrization, match, protein_pdb):
@@ -340,3 +371,14 @@ def test_example_csv_parses(path):
assert errors == {}
assert len(models) == len(df) > 0
assert {m.parametrization.family for m in models} <= {"charmm", "amber", "openff"}
+
+
+def test_simulation_type_tables_cover_the_literal():
+ """Per-type tables are keyed by the one SimulationType definition."""
+ assert set(type_mapping) == ALL_SIMULATION_TYPES
+ assert set(DEFAULT_PRODUCTION_LENGTH_NS) == ALL_SIMULATION_TYPES
+ assert set(COUPLING) == ALL_SIMULATION_TYPES
+ assert PROTEIN_SIMULATION_TYPES <= ALL_SIMULATION_TYPES
+ assert SMALL_MOLECULE_SIMULATION_TYPES <= ALL_SIMULATION_TYPES
+ assert CharmmParametrization.simulation_types == ALL_SIMULATION_TYPES
+ assert OpenFFParametrization.simulation_types.isdisjoint(PROTEIN_SIMULATION_TYPES)
diff --git a/mdfactory/tests/test_protein_mixedbox.py b/mdfactory/tests/test_protein_mixedbox.py
index 8cf0525..5f8f9ee 100644
--- a/mdfactory/tests/test_protein_mixedbox.py
+++ b/mdfactory/tests/test_protein_mixedbox.py
@@ -44,7 +44,6 @@
simulate_compression_until_density_reached,
)
from mdfactory.utils.setup_utilities import (
- N_AVOGADRO,
create_protein_mixedbox_universe,
cubic_box_edge_for_density,
native_atom_names,
@@ -57,10 +56,11 @@
unique_atom_names,
write_itp_with_atom_names,
)
+from mdfactory.utils.units import N_AVOGADRO
-from .conftest import protein_input_data
+from .conftest import EXAMPLES_DIR, make_forcefield_dir, protein_input_data
-EXAMPLES_DIR = Path(__file__).resolve().parents[2] / "examples" / "protein_mixedbox"
+MIXEDBOX_EXAMPLES = EXAMPLES_DIR / "protein_mixedbox"
def _write_pdb2gmx_topology(path):
@@ -125,20 +125,6 @@ def test_native_species_require_charmm_resname(self, smiles, wrong_resname):
class TestCharmmSpec:
- @pytest.mark.parametrize(
- ("overrides", "match"),
- [
- ({"protein": "charmm27"}, "only the registered 'charmm36m'"),
- ({"water_model": "spce"}, "only 'tip3p'"),
- ],
- )
- def test_rejects_unpinned_choices(self, protein_pdb, overrides, match):
- data = protein_input_data(
- "protein_mixedbox", protein_pdb, parametrization={"family": "charmm", **overrides}
- )
- with pytest.raises(ValueError, match=match):
- BuildInput(**data)
-
def test_resolves_cgenff_small_molecules(self, protein_pdb):
spec = BuildInput(**protein_input_data("protein_mixedbox", protein_pdb)).parametrization
assert spec.small_molecules.forcefield == "cgenff"
@@ -245,14 +231,14 @@ class TestExampleFiles:
],
)
def test_yaml_loads(self, filename, sizing_cls, total_count):
- inp = BuildInput(**yaml.safe_load((EXAMPLES_DIR / filename).read_text()))
+ inp = BuildInput(**yaml.safe_load((MIXEDBOX_EXAMPLES / filename).read_text()))
assert isinstance(inp.system, ProteinMixedBoxComposition)
assert isinstance(inp.parametrization, CharmmParametrization)
assert isinstance(inp.system.sizing, sizing_cls)
assert inp.system.total_count == total_count
def test_csv_roundtrip(self):
- df = pd.read_csv(EXAMPLES_DIR / "protein_mixedbox.csv")
+ df = pd.read_csv(MIXEDBOX_EXAMPLES / "protein_mixedbox.csv")
models, errors = df_to_build_input_models(df)
assert errors == {}
assert len(models) == 1
@@ -492,10 +478,9 @@ def test_native_charmm_atom_names(self):
def test_native_parameters_come_from_selected_forcefield(self, tmp_path, monkeypatch):
monkeypatch.setattr(settings, "parameter_store", tmp_path / "parameters")
- ff_dir = tmp_path / "selected.ff"
- ff_dir.mkdir()
- for filename in ("forcefield.itp", "tip3p.itp", "ions.itp"):
- (ff_dir / filename).touch()
+ ff_dir = make_forcefield_dir(
+ tmp_path, "selected.ff", ["forcefield.itp", "tip3p.itp", "ions.itp"]
+ )
native = NativeForceField(name="charmm36m", directory=ff_dir, water_model="tip3p")
spec = CharmmParametrization().resolve("protein_mixedbox")
@@ -609,7 +594,7 @@ def test_rmsd_restraint_is_translation_invariant(self):
system = mm.System()
for _ in positions:
system.addParticle(12.0)
- _add_protein_rmsd_restraint(system, positions, [0, 1, 2], 1000.0)
+ _add_protein_rmsd_restraint(system, positions, [0, 1, 2])
integrator = mm.VerletIntegrator(0.001 * unit.picoseconds)
context = mm.Context(system, integrator)
context.setPositions(positions)
@@ -721,10 +706,10 @@ def test_native_only_ordering_and_counts(self, tmp_path):
def test_cgenff_small_molecule_merges_params_and_includes_itp(self, tmp_path):
# A force field that already defines CG331, so it is de-duplicated out of
# the merged small-molecule parameters.
- ff_dir = tmp_path / "charmm36m.ff"
- ff_dir.mkdir()
- (ff_dir / "atomtypes.itp").write_text(
- "[ atomtypes ]\nCG331 6 12.011 0.000 A 0.36 0.28\n"
+ ff_dir = make_forcefield_dir(
+ tmp_path,
+ "charmm36m.ff",
+ {"atomtypes.itp": "[ atomtypes ]\nCG331 6 12.011 0.000 A 0.36 0.28\n"},
)
# The parameter files live in a separate source directory (as in the
@@ -806,17 +791,20 @@ def test_zero_count_species_skipped(self, tmp_path):
def _amber_ff_dir(tmp_path):
"""Write a minimal amber14sb.ff whose forcefield.itp and atom types mirror GROMACS 2026."""
- ff_dir = tmp_path / "amber14sb.ff"
- ff_dir.mkdir()
- (ff_dir / "forcefield.itp").write_text(AMBER_FORCEFIELD_ITP)
- (ff_dir / "ffnonbonded.itp").write_text(
- "[ atomtypes ]\n"
- " CT 6 12.01 0.0000 A 0.339966950842 0.457729600000\n"
- " OW_tip3p 8 16.0 0.0000 A 0.315075240658 0.635968000000\n"
+ return make_forcefield_dir(
+ tmp_path,
+ "amber14sb.ff",
+ {
+ "forcefield.itp": AMBER_FORCEFIELD_ITP,
+ "ffnonbonded.itp": (
+ "[ atomtypes ]\n"
+ " CT 6 12.01 0.0000 A 0.339966950842 0.457729600000\n"
+ " OW_tip3p 8 16.0 0.0000 A 0.315075240658 0.635968000000\n"
+ ),
+ "tip3p.itp": "; tip3p.itp\n",
+ "ions_tip3p.itp": "; ions_tip3p.itp\n",
+ },
)
- for name in ("tip3p.itp", "ions_tip3p.itp"):
- (ff_dir / name).write_text(f"; {name}\n")
- return ff_dir
def _sage_parameters(tmp_path, name="ETH0", **fmt):
@@ -1068,13 +1056,13 @@ class TestAmberExamplesAndSchedule:
],
)
def test_yaml_loads(self, filename, resnames):
- inp = BuildInput(**yaml.safe_load((EXAMPLES_DIR / filename).read_text()))
+ inp = BuildInput(**yaml.safe_load((MIXEDBOX_EXAMPLES / filename).read_text()))
assert inp.parametrization == AMBER_PROTEIN_MIXEDBOX
assert {s.resname for s in inp.system.species} == resnames
def test_csv_loads(self):
models, errors = df_to_build_input_models(
- pd.read_csv(EXAMPLES_DIR / "protein_mixedbox_amber.csv")
+ pd.read_csv(MIXEDBOX_EXAMPLES / "protein_mixedbox_amber.csv")
)
assert errors == {}
assert [m.parametrization for m in models] == [AMBER_PROTEIN_MIXEDBOX]
@@ -1110,3 +1098,19 @@ def test_finalize_writes_scheme_schedule_and_runs_grompp(
# grompp validates the generated em.mdp unconditionally.
assert len(calls) == 1
assert Path(calls[0][0][2]).name == "em.mdp"
+
+
+def test_read_include_layout_finds_anchors():
+ from mdfactory.parametrize import _read_include_layout
+
+ lines = [
+ '#include "amber14sb.ff/forcefield.itp"\n',
+ '#include "topol_Protein.itp"\n',
+ '#include "amber14sb.ff/tip3p.itp"\n',
+ "[ system ]\n",
+ "Protein\n",
+ "[ molecules ]\n",
+ "Protein 1\n",
+ ]
+ assert _read_include_layout(lines) == (0, 2, 3, 5)
+ assert _read_include_layout(lines[1:2]) == (None, None, None, None)
diff --git a/mdfactory/tests/test_proteinbox.py b/mdfactory/tests/test_proteinbox.py
index 090741a..2f2cd6c 100644
--- a/mdfactory/tests/test_proteinbox.py
+++ b/mdfactory/tests/test_proteinbox.py
@@ -15,8 +15,14 @@
import pytest
import yaml
-from mdfactory import workflows
-from mdfactory.build import _center_in_cubic_box, _ion_counts, _prepare_protein
+from mdfactory import build, cli
+from mdfactory.build import (
+ _bounding_box,
+ _center_in_cubic_box,
+ _ion_counts,
+ _prepare_protein,
+ _solvate_with_water,
+)
from mdfactory.models.composition import ProteinBoxComposition, ProteinMixedBoxComposition
from mdfactory.models.input import BuildInput
from mdfactory.models.parametrization import (
@@ -30,6 +36,7 @@
from mdfactory.setup.protein import (
_apply_protonation_states,
_build_disulfide_prompt_input,
+ _parse_protein_molecule_names,
_sum_atoms_column,
bundle_forcefield_into_topology,
check_forcefield_available,
@@ -40,9 +47,8 @@
update_topology_molecules,
)
-from .conftest import protein_input_data
+from .conftest import EXAMPLES_DIR, PROTEIN_TYPES, protein_input_data
-PROTEIN_SIM_TYPES = ["proteinbox", "protein_mixedbox"]
COMPOSITION_CLASSES = {
"proteinbox": ProteinBoxComposition,
"protein_mixedbox": ProteinMixedBoxComposition,
@@ -98,6 +104,7 @@ def _run(tmp_path, pdb, **kwargs):
kwargs.setdefault("merge_all", False)
kwargs.setdefault("disulfide_bonds", [])
kwargs.setdefault("protonation_states", {})
+ kwargs.setdefault("chains", [])
kwargs.setdefault("output_dir", tmp_path / "out")
return run_pdb2gmx(pdb_path=pdb, **kwargs)
@@ -172,31 +179,6 @@ def test_creation_with_custom_ionization(self, protein_pdb):
assert comp.ionization.concentration == 0.10
-@pytest.mark.parametrize("sim_type", PROTEIN_SIM_TYPES)
-class TestCharmmProteinSpec:
- """CHARMM protein force-field and water constraints for both protein types."""
-
- def test_frozen(self, sim_type, protein_pdb):
- inp = BuildInput(**protein_input_data(sim_type, protein_pdb))
- with pytest.raises(Exception):
- inp.parametrization.protein = "other"
-
- @pytest.mark.parametrize(
- ("overrides", "match"),
- [
- ({"protein": "amber99sb-ildn"}, "CHARMM"),
- ({"water_model": "tip4p"}, "3-site"),
- ({"protein": "charmm36m-ljpme"}, "LJ-PME"),
- ],
- )
- def test_rejects_unsupported_choices(self, sim_type, protein_pdb, overrides, match):
- data = protein_input_data(
- sim_type, protein_pdb, parametrization={"family": "charmm", **overrides}
- )
- with pytest.raises(ValueError, match=match):
- BuildInput(**data)
-
-
class TestForceFieldCheck:
def test_available_forcefield_passes(self, tmp_path, gmx_search_paths):
(tmp_path / "charmm27.ff").mkdir()
@@ -258,7 +240,7 @@ def test_rejects_mismatched_sha256(self, tmp_path, monkeypatch):
assert not (tmp_path / "test.ff").exists()
-@pytest.mark.parametrize("sim_type", PROTEIN_SIM_TYPES)
+@pytest.mark.parametrize("sim_type", PROTEIN_TYPES)
class TestProteinInput:
"""BuildInput validation shared by proteinbox and protein_mixedbox."""
@@ -389,7 +371,7 @@ def test_update_topology_molecules(self, tmp_path):
(
[("HIS", 15), ("HIS", 15), ("ALA", 16), ("GLU", 35)],
{"HIS15": "HIE", "GLU35": "GLH"},
- None,
+ "amber14sb",
["HIE", "HIE", "ALA", "GLH"],
),
# CHARMM histidine aliases.
@@ -397,45 +379,48 @@ def test_update_topology_molecules(self, tmp_path):
# An existing tautomer (HID) still matches an unqualified HIS key.
([("HID", 15)], {"HIS15": "HIE"}, "charmm36m", ["HSE"]),
# A chain-qualified key renames only the residue in that chain.
- ([("HIS", 15, "A"), ("HIS", 15, "B")], {"A:HIS15": "HIE"}, None, ["HIE", "HIS"]),
+ ([("HIS", 15, "A"), ("HIS", 15, "B")], {"A:HIS15": "HIE"}, "amber14sb", ["HIE", "HIS"]),
# Insertion codes distinguish otherwise identical residue numbers.
(
[("HIS", 100, "H"), ("HIS", 100, "H", 0.0, "A")],
{"H:HIS100A": "HIE"},
- None,
+ "amber14sb",
["HIS", "HIE"],
),
(
[("HIS", 100, "H"), ("HIS", 100, "H", 0.0, "A")],
{"H:HIS100": "HID"},
- None,
+ "amber14sb",
["HID", "HIS"],
),
],
)
def test_apply_protonation_states(self, tmp_path, atoms, states, forcefield, expected):
pdb = _write_pdb(tmp_path, *atoms)
- kwargs = {"forcefield": forcefield} if forcefield else {}
- result = _apply_protonation_states(pdb, states, tmp_path, **kwargs)
+ result = _apply_protonation_states(pdb, states, tmp_path, forcefield)
lines = [ln for ln in result.read_text().splitlines() if ln.startswith("ATOM")]
assert [ln[17:20].strip() for ln in lines] == expected
@pytest.mark.parametrize(
("atoms", "states", "forcefield", "match"),
[
- ([("GLU", 35)], {"GLU35": "GLUP"}, None, "cannot be written to a PDB"),
+ ([("GLU", 35)], {"GLU35": "GLUP"}, "amber14sb", "cannot be written to a PDB"),
# GLH aliases to CHARMM GLUP (4 chars), which cannot fit a PDB column.
([("GLU", 35)], {"GLU35": "GLH"}, "charmm36m", "cannot be written to a PDB"),
- ([("HIS", 15)], {"HIS99": "HIE"}, None, "matched no residue"),
+ ([("HIS", 15)], {"HIS99": "HIE"}, "amber14sb", "matched no residue"),
# HIS15 exists in both chains; an unqualified key cannot say which.
- ([("HIS", 15, "A"), ("HIS", 15, "B")], {"HIS15": "HIE"}, None, "ambiguous across"),
+ (
+ [("HIS", 15, "A"), ("HIS", 15, "B")],
+ {"HIS15": "HIE"},
+ "amber14sb",
+ "ambiguous across",
+ ),
],
)
def test_apply_protonation_states_rejects(self, tmp_path, atoms, states, forcefield, match):
pdb = _write_pdb(tmp_path, *atoms)
- kwargs = {"forcefield": forcefield} if forcefield else {}
with pytest.raises(ValueError, match=match):
- _apply_protonation_states(pdb, states, tmp_path, **kwargs)
+ _apply_protonation_states(pdb, states, tmp_path, forcefield)
@pytest.mark.parametrize(
("atoms", "bonds", "expected"),
@@ -531,7 +516,6 @@ def test_run_pdb2gmx_multichain_collects_chain_includes(self, tmp_path, monkeypa
params = _run(tmp_path, pdb, chains=chains)
- assert params.chains == chains
assert sorted(p.name for p in params.topology_include_files) == sorted(extra)
@pytest.mark.parametrize(
@@ -696,7 +680,7 @@ def test_models_internal_gaps_but_not_terminal_tails(self, tmp_path):
assert "internal" in text
-@pytest.mark.parametrize("sim_type", PROTEIN_SIM_TYPES)
+@pytest.mark.parametrize("sim_type", PROTEIN_TYPES)
class TestPdbPathResolution:
def test_relative_pdb_path_resolved_against_yaml_dir(self, tmp_path, sim_type):
base = tmp_path / "specs"
@@ -757,15 +741,13 @@ def test_csv_push_and_prepare_compute_identical_hashes(
assert prepared[0].system.protein.pdb_path == pdb.resolve()
assert [m.hash for m in pushed] == [m.hash for m in prepared]
- @pytest.mark.parametrize("sim_type", PROTEIN_SIM_TYPES)
- def test_run_from_file_resolves_relative_path(
- self, tmp_path, protein_pdb, monkeypatch, sim_type
- ):
+ @pytest.mark.parametrize("sim_type", PROTEIN_TYPES)
+ def test_yaml_build_resolves_relative_path(self, tmp_path, protein_pdb, monkeypatch, sim_type):
_, pdb, yaml_file, _ = _write_spec_dir(tmp_path, protein_pdb, sim_type)
- captured = {}
- monkeypatch.setattr(workflows, "run_build_from_dict", captured.update)
- workflows.run_build_from_file(yaml_file)
- assert captured["system"]["protein"]["pdb_path"] == str(pdb.resolve())
+ captured = []
+ monkeypatch.setattr(cli, "run_build_from_dict", captured.append)
+ cli._build_from_yaml(yaml_file, tmp_path / "out")
+ assert captured[0].system.protein.pdb_path == pdb.resolve()
def test_yaml_dry_run_resolves_before_dispatch(self, tmp_path, protein_pdb, monkeypatch):
import mdfactory.orchestration
@@ -834,6 +816,45 @@ def test_padding_guaranteed_on_every_face(self):
assert abs((box_size - bbox_max[long_axis]) - padding) < 1e-3
+def test_bounding_box_returns_corners_and_center():
+ positions = np.array([[0.0, 2.0, -4.0], [10.0, 6.0, 0.0], [3.0, 4.0, 2.0]])
+ lower, upper, center = _bounding_box(positions)
+ assert lower.tolist() == [0.0, 2.0, -4.0]
+ assert upper.tolist() == [10.0, 6.0, 2.0]
+ assert center.tolist() == [5.0, 4.0, -1.0]
+
+
+class TestSolvateWithWater:
+ @staticmethod
+ def _solvated(n_water):
+ u = mda.Universe.empty(
+ n_water, n_residues=n_water, atom_resindex=np.arange(n_water), trajectory=True
+ )
+ u.add_TopologyAttr("resnames", ["SOL"] * n_water)
+ u.dimensions = [50.0, 50.0, 50.0, 90, 90, 90]
+ return u
+
+ def test_returns_water_species_and_keeps_the_box(self, monkeypatch):
+ solvated = self._solvated(4)
+ solvated.dimensions = None
+ monkeypatch.setattr(build, "solvate", lambda u, **kwargs: solvated)
+ reference = self._solvated(1)
+ reference.dimensions = [20.0, 20.0, 20.0, 90, 90, 90]
+
+ result, water = _solvate_with_water(reference, prune_in_z=False)
+
+ assert result is solvated
+ assert list(result.dimensions) == [20.0, 20.0, 20.0, 90, 90, 90]
+ assert (water.smiles, water.resname, water.count) == ("O", "SOL", 4)
+
+ def test_no_water_raises(self, monkeypatch):
+ dry = self._solvated(1)
+ dry.residues.resnames = ["LIP"]
+ monkeypatch.setattr(build, "solvate", lambda u, **kwargs: dry)
+ with pytest.raises(ValueError, match="did not add any water"):
+ _solvate_with_water(self._solvated(1))
+
+
class TestIonCounts:
"""ionize_solvated_system returns only the ions it added; counts map by resname."""
@@ -855,11 +876,11 @@ def test_ion_counts(self, ions, expected):
assert _ion_counts(species) == expected
-EXAMPLES_DIR = Path(__file__).resolve().parents[2] / "examples" / "proteinbox"
+PROTEINBOX_EXAMPLES = EXAMPLES_DIR / "proteinbox"
class TestAmberSpec:
- @pytest.mark.parametrize("sim_type", PROTEIN_SIM_TYPES)
+ @pytest.mark.parametrize("sim_type", PROTEIN_TYPES)
def test_family_defaults(self, sim_type, protein_pdb):
inp = BuildInput(**protein_input_data(sim_type, protein_pdb, parametrization="amber"))
spec = inp.parametrization
@@ -872,28 +893,6 @@ def test_family_defaults(self, sim_type, protein_pdb):
# A proteinbox has no small molecules, so the component is not shown.
assert spec.small_molecules is None
- def test_frozen(self, protein_pdb):
- inp = BuildInput(**protein_input_data("proteinbox", protein_pdb, parametrization="amber"))
- with pytest.raises(Exception):
- inp.parametrization.protein = "amber19sb"
-
- @pytest.mark.parametrize(
- ("overrides", "match"),
- [
- ({"protein": "amber19sb"}, "Only 'amber14sb'"),
- ({"protein": "charmm36m"}, "Only 'amber14sb'"),
- ({"water_model": "opc"}, "Only 'tip3p'"),
- ({"water_model": "spce"}, "Only 'tip3p'"),
- ],
- )
- @pytest.mark.parametrize("sim_type", PROTEIN_SIM_TYPES)
- def test_rejects_unsupported_choices(self, sim_type, protein_pdb, overrides, match):
- data = protein_input_data(
- sim_type, protein_pdb, parametrization={"family": "amber", **overrides}
- )
- with pytest.raises(ValueError, match=match):
- BuildInput(**data)
-
class TestAmberInput:
def test_small_molecule_override_for_protein_mixedbox(self, protein_pdb):
@@ -954,9 +953,13 @@ def test_example_yaml_and_csv(self):
from mdfactory.prepare import df_to_build_input_models
- inp = BuildInput(**yaml.safe_load((EXAMPLES_DIR / "lysozyme_amber14sb.yaml").read_text()))
+ inp = BuildInput(
+ **yaml.safe_load((PROTEINBOX_EXAMPLES / "lysozyme_amber14sb.yaml").read_text())
+ )
assert isinstance(inp.parametrization, AmberParametrization)
- models, errors = df_to_build_input_models(pd.read_csv(EXAMPLES_DIR / "proteins_amber.csv"))
+ models, errors = df_to_build_input_models(
+ pd.read_csv(PROTEINBOX_EXAMPLES / "proteins_amber.csv")
+ )
assert errors == {}
assert [m.parametrization.family for m in models] == ["amber"]
assert models[0].parametrization == inp.parametrization
@@ -1064,14 +1067,13 @@ def test_prepare_protein_takes_forcefield_from_spec_and_options_from_protein(
resname="LYZ", pdb_path=protein_pdb, ignore_hydrogens=False, merge_all=True
)
- params = _prepare_protein(species, spec)
+ _prepare_protein(species, spec)
cmd = calls["cmd"]
assert cmd[cmd.index("-ff") + 1] == "charmm27"
assert cmd[cmd.index("-water") + 1] == "spce"
assert "-ignh" not in cmd
assert cmd[cmd.index("-merge") + 1] == "all"
- assert (params.forcefield, params.water_model) == ("charmm27", "spce")
def test_run_pdb2gmx_amber_passes_forcefield_and_amber_protonation(tmp_path, monkeypatch):
@@ -1079,7 +1081,7 @@ def test_run_pdb2gmx_amber_passes_forcefield_and_amber_protonation(tmp_path, mon
calls = {}
_patch_pdb2gmx(monkeypatch, ["Protein_chain_A"], ["posre.itp"], calls)
- params = _run(
+ _run(
tmp_path,
pdb,
forcefield="amber14sb",
@@ -1095,4 +1097,13 @@ def test_run_pdb2gmx_amber_passes_forcefield_and_amber_protonation(tmp_path, mon
resnames = [ln[17:20] for ln in adjusted.read_text().splitlines() if ln.startswith("ATOM")]
# Amber names pass through unchanged (no CHARMM HSE alias).
assert resnames == ["HIE", "ALA"]
- assert params.forcefield == "amber14sb"
+
+
+def test_parse_protein_molecule_names_lists_only_protein_entries(tmp_path):
+ top = tmp_path / "topol.top"
+ top.write_text(
+ "[ system ]\nProtein_chain_X in water\n\n"
+ "[ molecules ]\n; Compound #mols\n"
+ "Protein_chain_A 1\nProtein_chain_B 1 ; second\nSOL 100\nNA 2\n"
+ )
+ assert _parse_protein_molecule_names(top) == ["Protein_chain_A", "Protein_chain_B"]
diff --git a/mdfactory/tests/test_push.py b/mdfactory/tests/test_push.py
index 7eaaf60..68f7cf8 100644
--- a/mdfactory/tests/test_push.py
+++ b/mdfactory/tests/test_push.py
@@ -28,7 +28,12 @@
)
# Import shared helpers from conftest
-from .conftest import _create_simulation_folder, _write_csv_with_models, make_test_record
+from .conftest import (
+ _create_simulation_folder,
+ _write_csv_with_models,
+ default_settings_loader,
+ make_test_record,
+)
def _record(hash_value: str, status: str = "build") -> dict[str, str]:
@@ -811,12 +816,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
results = init_systems_database()
@@ -847,12 +847,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
results_first = init_systems_database()
results_second = init_systems_database()
@@ -871,13 +866,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return str(db_path)
return original_get_db_path(self, db_name, db_type)
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
# Init and insert a row
init_systems_database()
@@ -907,12 +897,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
monkeypatch.chdir(tmp_path)
sync_init_systems()
diff --git a/mdfactory/tests/test_push_analysis.py b/mdfactory/tests/test_push_analysis.py
index 21b0a25..d2791f1 100644
--- a/mdfactory/tests/test_push_analysis.py
+++ b/mdfactory/tests/test_push_analysis.py
@@ -33,7 +33,7 @@
)
# Import shared helpers from conftest (temp_analysis_db fixture is auto-discovered)
-from .conftest import _write_csv_with_models
+from .conftest import _write_csv_with_models, default_settings_loader
def _create_simulation_folder_analysis(tmp_path: Path, hash_value: str, build_input) -> Path:
@@ -260,12 +260,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
results = init_analysis_database()
@@ -551,13 +546,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return str(db_path)
return original_get_db_path(self, db_name, db_type)
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
init_analysis_database()
# Create simulation folder with analysis data
@@ -594,13 +584,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return original_get_db_path(self, db_name, db_type)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
init_analysis_database()
# Create simulation folder with analysis data
@@ -629,13 +614,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return original_get_db_path(self, db_name, db_type)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
init_analysis_database()
# Create simulation folder with analysis
@@ -663,13 +643,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return str(db_path)
return original_get_db_path(self, db_name, db_type)
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
init_analysis_database()
_, models = _write_csv_with_models(tmp_path)
@@ -706,12 +681,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
monkeypatch.chdir(tmp_path)
sync_init_analysis()
@@ -784,13 +754,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return original_get_db_path(self, db_name, db_type)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
init_analysis_database()
_, models = _write_csv_with_models(tmp_path)
@@ -839,13 +804,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return str(db_path)
return original_get_db_path(self, db_name, db_type)
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
init_analysis_database()
_, models = _write_csv_with_models(tmp_path)
@@ -884,13 +844,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return str(db_path)
return original_get_db_path(self, db_name, db_type)
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
results = init_analysis_database()
@@ -928,13 +883,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return str(db_path)
return original_get_db_path(self, db_name, db_type)
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
# Init and insert data
init_analysis_database()
@@ -970,13 +920,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return str(db_path)
return original_get_db_path(self, db_name, db_type)
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
results_first = init_analysis_database()
results_second = init_analysis_database()
diff --git a/mdfactory/tests/test_push_artifacts.py b/mdfactory/tests/test_push_artifacts.py
index eee33d7..a6d4eb8 100644
--- a/mdfactory/tests/test_push_artifacts.py
+++ b/mdfactory/tests/test_push_artifacts.py
@@ -28,7 +28,7 @@
)
# Import shared helpers from conftest (temp_analysis_db-like fixture will be temp_artifact_db)
-from .conftest import _write_csv_with_models
+from .conftest import _write_csv_with_models, default_settings_loader
def _create_simulation_folder_artifact(tmp_path: Path, hash_value: str, build_input) -> Path:
@@ -138,12 +138,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
# Initialize both analysis and artifact database
init_analysis_database()
@@ -167,12 +162,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
results = init_artifact_database()
@@ -204,12 +194,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
init_analysis_database()
init_artifact_database()
@@ -338,12 +323,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
init_artifact_database()
# Create simulation folder with artifact data
@@ -386,12 +366,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
init_analysis_database()
init_artifact_database()
@@ -427,12 +402,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
init_analysis_database()
init_artifact_database()
@@ -510,12 +480,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
monkeypatch.chdir(tmp_path)
sync_init_artifacts()
@@ -590,12 +555,7 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
# Force sqlite backend
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
init_analysis_database()
init_artifact_database()
@@ -649,13 +609,8 @@ def patched_get_db_path(self, db_name: str, db_type: str = "sqlite") -> str:
return str(db_path)
return original_get_db_path(self, db_name, db_type)
- def patched_init(self):
- for section, options in Settings.DEFAULT_CONFIG.items():
- self.config[section] = dict(options)
- self.config["database"]["TYPE"] = "sqlite"
-
monkeypatch.setattr(Settings, "get_db_path", patched_get_db_path)
- monkeypatch.setattr(Settings, "_load_config", patched_init)
+ monkeypatch.setattr(Settings, "_load_config", default_settings_loader("sqlite"))
init_analysis_database()
init_artifact_database()
diff --git a/mdfactory/tests/test_run_schedules.py b/mdfactory/tests/test_run_schedules.py
index 1c096ca..4fa9903 100644
--- a/mdfactory/tests/test_run_schedules.py
+++ b/mdfactory/tests/test_run_schedules.py
@@ -27,6 +27,8 @@
write_mdp_files,
)
+from .conftest import PROTEIN_TYPES, WATER
+
# Every simulation type with each family that supports it (both nonbonded schemes).
COMBINATIONS = [
("mixedbox", CharmmParametrization),
@@ -40,7 +42,6 @@
("protein_mixedbox", CharmmParametrization),
("protein_mixedbox", AmberParametrization),
]
-PROTEIN_TYPES = ("proteinbox", "protein_mixedbox")
def _inp(simulation_type, family=CharmmParametrization, **settings):
@@ -348,12 +349,6 @@ def test_sage_scheme_matches_the_openff_force_fields(tmp_path):
# --- re-render existing directories -----------------------------------------------------
-WATER = {
- "simulation_type": "mixedbox",
- "parametrization": "charmm",
- "system": {"species": [{"smiles": "O", "count": 100, "resname": "SOL"}]},
-}
-
@pytest.fixture()
def prepared(tmp_path):
diff --git a/mdfactory/tests/test_simulation_settings.py b/mdfactory/tests/test_simulation_settings.py
index 184a799..1258cc7 100644
--- a/mdfactory/tests/test_simulation_settings.py
+++ b/mdfactory/tests/test_simulation_settings.py
@@ -2,7 +2,6 @@
# ABOUTME: Covers defaults, validation, hashing, YAML/CSV/file/object inputs, precedence, duplicates
"""Tests for :mod:`mdfactory.models.simulation_settings` and its BuildInput integration."""
-from pathlib import Path
from unittest.mock import patch
import pandas as pd
@@ -23,13 +22,7 @@
load_build_input_csv,
)
-WATER = {
- "simulation_type": "mixedbox",
- "parametrization": "charmm",
- "system": {"species": [{"smiles": "O", "count": 100, "resname": "SOL"}]},
-}
-
-EXAMPLES = Path(__file__).resolve().parents[2] / "examples"
+from .conftest import EXAMPLES_DIR, WATER
CSV_HEADER = "simulation_type,parametrization,system.species.SOL.smiles,system.species.SOL.count"
@@ -121,7 +114,7 @@ def test_metadata_and_record_carry_settings(tmp_path):
def test_record_round_trip_keeps_hash_of_dict_fields(tmp_path):
"""Dict fields (LNP input_fractions) keep their order, so the record keeps the hash."""
- data = yaml.safe_load((EXAMPLES / "lnp" / "lnp_alc0315_smirnoff.yaml").read_text())
+ data = yaml.safe_load((EXAMPLES_DIR / "lnp" / "lnp_alc0315_smirnoff.yaml").read_text())
model = BuildInput(**data)
fractions = model.system.shell.input_fractions
assert list(fractions) != sorted(fractions) # the case that a sorted dump breaks
diff --git a/mdfactory/tests/test_sync_config_local_paths.py b/mdfactory/tests/test_sync_config_local_paths.py
index e4a3127..ff51921 100644
--- a/mdfactory/tests/test_sync_config_local_paths.py
+++ b/mdfactory/tests/test_sync_config_local_paths.py
@@ -50,7 +50,7 @@ def test_config_resolves_relative_local_paths_against_config_directory(monkeypat
def patched_load_config(self):
self.config = configparser.ConfigParser()
- for section, options in self.DEFAULT_CONFIG.items():
+ for section, options in self._get_defaults().items():
self.config[section] = {}
for key, value in options.items():
self.config[section][key] = str(value)
diff --git a/mdfactory/tests/test_topology_utilities.py b/mdfactory/tests/test_topology_utilities.py
index 76624aa..f03987d 100644
--- a/mdfactory/tests/test_topology_utilities.py
+++ b/mdfactory/tests/test_topology_utilities.py
@@ -13,9 +13,13 @@
from mdfactory.utils.topology_utilities import (
_atomtype_signature,
_get_sybyl_atom_type,
+ _iter_prm_atomtypes,
collect_forcefield_atomtype_definitions,
extract_reusable_parts_from_cgenff_gmx_top,
+ iter_includes,
+ iter_topology_entries,
merge_extra_parameter_itps,
+ read_topology_defaults,
rename_atomtypes,
validate_cgenff_atomtype_compatibility,
write_mol2_from_rdkit,
@@ -139,7 +143,9 @@ def test_cgenff_atomtype_compatibility_ignores_inline_comment(tmp_path):
generated.write_text("[ atomtypes ]\nCG331 6 12.011 0.0 A 0.36 0.28 ; sp3 carbon\n")
# Identical physical parameters, only a comment differs -> no conflict raised.
- validate_cgenff_atomtype_compatibility([generated], ff_dir)
+ validate_cgenff_atomtype_compatibility(
+ [generated], collect_forcefield_atomtype_definitions(ff_dir)
+ )
def test_cgenff_atomtype_compatibility_rejects_conflict(tmp_path):
@@ -150,7 +156,9 @@ def test_cgenff_atomtype_compatibility_rejects_conflict(tmp_path):
generated.write_text("[ atomtypes ]\nCG331 6 12.011 0.0 A 0.40 0.28\n")
with pytest.raises(ValueError, match="conflicts with the selected CHARMM"):
- validate_cgenff_atomtype_compatibility([generated], ff_dir)
+ validate_cgenff_atomtype_compatibility(
+ [generated], collect_forcefield_atomtype_definitions(ff_dir)
+ )
def test_collect_atomtypes_honors_preprocessor_conditionals(tmp_path):
@@ -435,3 +443,89 @@ def test_oco2_carboxylate_partner():
o_neutral = mol.GetAtomWithIdx(2)
assert o_neutral.GetFormalCharge() == 0
assert _get_sybyl_atom_type(o_neutral) == "O.co2"
+
+
+ENTRY_TEXT = """\
+; title
+#define FLEXIBLE
+[ Atoms ]
+; nr type
+1 CT ; first
+#ifdef X
+2 CT
+#endif
+
+[ bonds ] ; trailing comment
+1 2
+"""
+
+
+def test_iter_topology_entries_tracks_sections_and_keeps_comments():
+ entries = list(iter_topology_entries(ENTRY_TEXT.splitlines(keepends=True)))
+ # Blank, comment-only and preprocessor lines are not data; sections are lower-cased.
+ assert [(e.section, e.code) for e in entries] == [
+ ("atoms", "1 CT"),
+ ("atoms", "2 CT"),
+ ("bonds", "1 2"),
+ ]
+ assert entries[0].comment == "; first"
+ assert entries[1].comment == ""
+ # The index addresses the original line, so a caller can rewrite it.
+ assert ENTRY_TEXT.splitlines()[entries[2].index] == "1 2"
+
+
+def test_iter_topology_entries_headers_are_opt_in():
+ lines = ENTRY_TEXT.splitlines()
+ headers = [e for e in iter_topology_entries(lines, headers=True) if e.is_header]
+ assert [(e.section, e.index) for e in headers] == [("atoms", 2), ("bonds", 9)]
+ assert not any(e.is_header for e in iter_topology_entries(lines))
+
+
+def test_iter_topology_entries_supports_early_stop():
+ seen = []
+ for entry in iter_topology_entries(ENTRY_TEXT.splitlines()):
+ if entry.section != "atoms":
+ break
+ seen.append(entry.code)
+ assert seen == ["1 CT", "2 CT"]
+
+
+def test_iter_includes_reports_paths_and_line_indices():
+ lines = [
+ '#include "charmm36.ff/forcefield.itp"\n',
+ '; #include "commented.itp"\n',
+ ' #include "local.itp" ; trailing\n',
+ "#include \n",
+ "[ system ]\n",
+ ]
+ assert list(iter_includes(lines)) == [(0, "charmm36.ff/forcefield.itp"), (2, "local.itp")]
+
+
+def test_read_topology_defaults(tmp_path):
+ top = tmp_path / "forcefield.itp"
+ top.write_text(
+ "[ defaults ]\n; nbfunc comb\n1 2 YES 0.5 0.8333 ; amber\n[ atomtypes ]\nX 1\n"
+ )
+ assert read_topology_defaults(top) == (1, 2, "yes", 0.5, 0.8333)
+
+ top.write_text("[ atomtypes ]\nX 1\n")
+ assert read_topology_defaults(top) is None
+
+ top.write_text("[ defaults ]\n[ atomtypes ]\nX 1\n")
+ with pytest.raises(ValueError, match="has no entry"):
+ read_topology_defaults(top)
+
+ top.write_text("[ defaults ]\n1 2 yes\n")
+ with pytest.raises(ValueError, match="Malformed"):
+ read_topology_defaults(top)
+
+
+def test_iter_prm_atomtypes_reads_every_file(tmp_path):
+ first = tmp_path / "a.itp"
+ first.write_text("[ atomtypes ]\n; name\nCG331 6 12.011 0.0 A 0.36 0.28\n")
+ second = tmp_path / "b.itp"
+ second.write_text("[ atomtypes ]\nNA1 11 22.99 1.0 A 0.2 0.1 ; sodium\n")
+ assert [(prm, name) for prm, name, _ in _iter_prm_atomtypes([first, second])] == [
+ (first, "CG331"),
+ (second, "NA1"),
+ ]
diff --git a/mdfactory/utils/setup_utilities.py b/mdfactory/utils/setup_utilities.py
index f9b4548..5194e79 100644
--- a/mdfactory/utils/setup_utilities.py
+++ b/mdfactory/utils/setup_utilities.py
@@ -20,16 +20,14 @@
from mdfactory.models.composition import BilayerComposition
from mdfactory.models.input import MixedBoxComposition
-from mdfactory.models.species import LipidSpecies
+from mdfactory.models.species import NATIVE_SPECIES, LipidSpecies
+from mdfactory.utils.units import A3_TO_L, CM3_TO_A3, N_AVOGADRO
from mdfactory.utils.utilities import temporary_working_directory, working_directory
-# Avogadro's number, matching the g<->Dalton conventions used in composition.py.
-N_AVOGADRO = 6.022e23
-
-# A protein_mixedbox is packed at 0.1x the final density so water and solutes have
-# room to be placed around the fixed protein before compression squeezes the box
-# down to its requested size (mirrors the mixedbox 0.1x loose-packing idiom).
-PROTEIN_MIXEDBOX_DILUTION_FACTOR = 0.1
+# Mixed boxes (and protein_mixedbox) are packed at 0.1x the final density so the
+# molecules have room to be placed before compression squeezes the box down to
+# its requested size.
+LOOSE_PACKING_FACTOR = 0.1
def generate_lipid_structure_legacy(
@@ -742,18 +740,18 @@ def protein_displaced_volume_a3(
``partial_specific_volume`` is in mL/g. The mass converts to grams via
Avogadro's number, multiplies by the specific volume to give cm³, and cm³
- converts to ų (1 cm³ = 1e24 ų).
+ converts to ų.
"""
mass_g = protein_mass_dalton / N_AVOGADRO
volume_cm3 = mass_g * partial_specific_volume
- return volume_cm3 * 1e24
+ return volume_cm3 * CM3_TO_A3
def cubic_box_edge_for_density(total_mass_dalton: float, target_density: float) -> float:
"""Return box edge (Å) holding ``total_mass_dalton`` at the density."""
mass_g = total_mass_dalton / N_AVOGADRO
volume_cm3 = mass_g / target_density
- volume_a3 = volume_cm3 * 1e24
+ volume_a3 = volume_cm3 * CM3_TO_A3
return float(volume_a3 ** (1.0 / 3.0))
@@ -815,8 +813,8 @@ def resolve_solution_counts(
v_basis_a3 = concentration_reference_volume_a3(
box_volume_a3, protein_mass_dalton, partial_specific_volume, basis
)
- # 1 ų = 1e-27 L, so count = concentration[mol/L] * V[L] * N_A.
- v_basis_liters = v_basis_a3 * 1e-27
+ # count = concentration[mol/L] * V[L] * N_A.
+ v_basis_liters = v_basis_a3 * A3_TO_L
counts = []
for spec in species:
if spec.count is not None:
@@ -988,11 +986,8 @@ def native_atom_names(smiles: str) -> list[str] | None:
Atom names in packing order, or None if the species is not native.
"""
- return {
- "O": ["OW", "HW1", "HW2"],
- "[Na+]": ["NA"],
- "[Cl-]": ["CL"],
- }.get(smiles)
+ native = NATIVE_SPECIES.get(smiles)
+ return list(native.atom_names) if native is not None else None
def topology_atom_names(species, parameter_set) -> list[str] | None:
diff --git a/mdfactory/utils/topology_utilities.py b/mdfactory/utils/topology_utilities.py
index 9ff5f22..8555ebf 100644
--- a/mdfactory/utils/topology_utilities.py
+++ b/mdfactory/utils/topology_utilities.py
@@ -6,9 +6,10 @@
import subprocess
import warnings
from collections import defaultdict
-from collections.abc import Collection
+from collections.abc import Collection, Iterable, Iterator
from itertools import groupby
from pathlib import Path
+from typing import NamedTuple
from rdkit import Chem
@@ -348,6 +349,63 @@ def run_cgenff_to_gmx(mol2_file):
raise RuntimeError("cgenff failed.") from e
+class TopologyEntry(NamedTuple):
+ """One line of a GROMACS topology as seen by :func:`iter_topology_entries`."""
+
+ index: int
+ #: Lower-cased name of the section the line belongs to ("" before the first header).
+ section: str
+ #: The line without its comment, stripped.
+ code: str
+ #: The ``;`` comment including the semicolon ("" when there is none).
+ comment: str
+ is_header: bool = False
+
+
+def iter_topology_entries(
+ lines: Iterable[str], *, headers: bool = False
+) -> Iterator[TopologyEntry]:
+ """Walk the data lines of GROMACS topology text, tracking the current section.
+
+ Blank lines, comment-only lines and preprocessor (``#``) lines are skipped.
+ Section names are lower-cased.
+
+ Parameters
+ ----------
+ lines : iterable of str
+ Topology lines, with or without line endings.
+ headers : bool, optional
+ Also yield the ``[ section ]`` header lines (``is_header=True``, with
+ the new section as ``section``). Default is False.
+
+ Yields
+ ------
+ TopologyEntry
+ Line index, section, comment-free code and comment of each data line.
+
+ """
+ section = ""
+ for index, raw in enumerate(lines):
+ code, sep, comment = raw.rstrip("\n").partition(";")
+ code = code.strip()
+ if not code or code.startswith("#"):
+ continue
+ is_header = code.startswith("[") and code.endswith("]")
+ if is_header:
+ section = code[1:-1].strip().lower()
+ if not headers:
+ continue
+ yield TopologyEntry(index, section, code, sep + comment if sep else "", is_header)
+
+
+def iter_includes(lines: Iterable[str]) -> Iterator[tuple[int, str]]:
+ """Yield ``(line index, path)`` for each ``#include "path"`` line."""
+ for index, raw in enumerate(lines):
+ stripped = raw.strip()
+ if stripped.startswith("#include") and '"' in stripped:
+ yield index, stripped.split('"')[1]
+
+
def count_contiguous_strings(str_list: list[str]) -> list[tuple[str, int]]:
"""Generate from a list of strings a list of tuples with the string
as the first element and the number of continuous repetitions as
@@ -555,21 +613,15 @@ def read_itp_atoms(itp_path: Path) -> list[tuple[str, float]]:
"""
atoms: list[tuple[str, float]] = []
- in_atoms = False
- for raw in Path(itp_path).read_text().splitlines():
- line = raw.split(";", 1)[0].strip()
- if not line or line.startswith("#"):
- continue
- if line.startswith("[") and line.endswith("]"):
- if in_atoms:
+ for entry in iter_topology_entries(Path(itp_path).read_text().splitlines()):
+ if entry.section != "atoms":
+ if atoms:
break
- in_atoms = line[1:-1].strip().lower() == "atoms"
continue
- if in_atoms:
- fields = line.split()
- if len(fields) < 8:
- raise ValueError(f"Atom line without a mass in {itp_path}: {raw!r}")
- atoms.append((fields[4], float(fields[7])))
+ fields = entry.code.split()
+ if len(fields) < 8:
+ raise ValueError(f"Atom line without a mass in {itp_path}: {entry.code!r}")
+ atoms.append((fields[4], float(fields[7])))
if not atoms:
raise ValueError(f"No [ atoms ] section found in {itp_path}.")
return atoms
@@ -631,36 +683,26 @@ def write_itp_with_atom_names(src: str | Path, dest: str | Path, names: list[str
"""
lines = Path(src).read_text().splitlines(keepends=True)
- in_atoms = False
- done = False
index = 0
- for i, raw in enumerate(lines):
- line = raw.split(";", 1)[0].strip()
- if done or not line or line.startswith("#"):
- continue
- if line.startswith("[") and line.endswith("]"):
- if in_atoms:
- done = True
- continue
- in_atoms = line[1:-1].strip().lower() == "atoms"
+ for entry in iter_topology_entries(lines):
+ if entry.section != "atoms":
+ if index:
+ break
continue
- if in_atoms:
- if index >= len(names):
- raise ValueError(f"{src} has more atoms than the {len(names)} names given.")
- fields = line.split()
- fields[4] = names[index]
- comment = raw.split(";", 1)[1].rstrip("\n") if ";" in raw else ""
- lines[i] = " ".join(fields) + (f" ;{comment}" if comment else "") + "\n"
- index += 1
+ if index >= len(names):
+ raise ValueError(f"{src} has more atoms than the {len(names)} names given.")
+ fields = entry.code.split()
+ fields[4] = names[index]
+ lines[entry.index] = _join_fields(fields, entry.comment)
+ index += 1
if index != len(names):
raise ValueError(f"{src} has {index} atoms but {len(names)} names were given.")
Path(dest).write_text("".join(lines))
-def _split_code_comment(raw: str) -> tuple[str, str]:
- """Split a topology line into its code and its ``;`` comment (without newline)."""
- code, sep, comment = raw.rstrip("\n").partition(";")
- return code, (sep + comment if sep else "")
+def _join_fields(fields: list[str], comment: str) -> str:
+ """Rebuild a topology line from its fields and trailing comment."""
+ return " ".join(fields) + (f" {comment}" if comment else "") + "\n"
def rename_atomtypes(itp_path: str | Path, atomtypes_path: str | Path, prefix: str) -> dict:
@@ -699,18 +741,10 @@ def rename_atomtypes(itp_path: str | Path, atomtypes_path: str | Path, prefix: s
atomtypes_lines = Path(atomtypes_path).read_text().splitlines(keepends=True)
mapping: dict[str, str] = {}
- section = None
- for i, raw in enumerate(atomtypes_lines):
- code, comment = _split_code_comment(raw)
- stripped = code.strip()
- if not stripped or stripped.startswith("#"):
+ for entry in iter_topology_entries(atomtypes_lines):
+ if entry.section != "atomtypes":
continue
- if stripped.startswith("[") and stripped.endswith("]"):
- section = stripped[1:-1].strip().lower()
- continue
- if section != "atomtypes":
- continue
- fields = stripped.split()
+ fields = entry.code.split()
old = fields[0]
if old in mapping:
raise ValueError(f"Atom type '{old}' is defined twice in {atomtypes_path}.")
@@ -719,24 +753,16 @@ def rename_atomtypes(itp_path: str | Path, atomtypes_path: str | Path, prefix: s
# An optional bond_type column repeats the name.
if len(fields) > 1 and fields[1] == old:
fields[1] = mapping[old]
- atomtypes_lines[i] = " ".join(fields) + (f" {comment}" if comment else "") + "\n"
+ atomtypes_lines[entry.index] = _join_fields(fields, entry.comment)
if not mapping:
raise ValueError(f"No [ atomtypes ] entries found in {atomtypes_path}.")
itp_lines = Path(itp_path).read_text().splitlines(keepends=True)
- section = None
n_atoms = 0
- for i, raw in enumerate(itp_lines):
- code, comment = _split_code_comment(raw)
- stripped = code.strip()
- if not stripped or stripped.startswith("#"):
- continue
- if stripped.startswith("[") and stripped.endswith("]"):
- section = stripped[1:-1].strip().lower()
- continue
- if section != "atoms":
+ for entry in iter_topology_entries(itp_lines):
+ if entry.section != "atoms":
continue
- fields = stripped.split()
+ fields = entry.code.split()
# nr type resnr residue atom cgnr charge mass [typeB chargeB massB]
for column in (1, 8):
if column >= len(fields):
@@ -747,7 +773,7 @@ def rename_atomtypes(itp_path: str | Path, atomtypes_path: str | Path, prefix: s
f"{atomtypes_path}."
)
fields[column] = mapping[fields[column]]
- itp_lines[i] = " ".join(fields) + (f" {comment}" if comment else "") + "\n"
+ itp_lines[entry.index] = _join_fields(fields, entry.comment)
n_atoms += 1
if n_atoms == 0:
raise ValueError(f"No [ atoms ] entries found in {itp_path}.")
@@ -757,18 +783,8 @@ def rename_atomtypes(itp_path: str | Path, atomtypes_path: str | Path, prefix: s
return mapping
-def validate_cgenff_atomtype_compatibility(
- prm_files: list[str | Path],
- ff_dir: str | Path,
- forcefield_types: dict[str, tuple[str, ...]] | None = None,
-) -> None:
- """Reject CGenFF atom types that conflict with the selected CHARMM bundle.
-
- *forcefield_types* may pass the already-parsed result of
- :func:`collect_forcefield_atomtype_definitions` to avoid re-reading *ff_dir*.
- """
- if forcefield_types is None:
- forcefield_types = collect_forcefield_atomtype_definitions(ff_dir)
+def _iter_prm_atomtypes(prm_files: list[str | Path]) -> Iterator[tuple[str | Path, str, str]]:
+ """Yield ``(file, atom-type name, definition line)`` for the ``[ atomtypes ]`` of each file."""
for prm in prm_files:
*_, sections = extract_reusable_parts_from_cgenff_gmx_top(prm, must_have_moleculetype=False)
for section_name, section_text in sections:
@@ -778,17 +794,29 @@ def validate_cgenff_atomtype_compatibility(
stripped = line.strip()
if not stripped or stripped.startswith((";", "[")):
continue
- name = stripped.split()[0]
- existing = forcefield_types.get(name)
- if existing is None:
- continue
- generated = _atomtype_signature(stripped)
- if generated != existing:
- raise ValueError(
- f"CGenFF atom type '{name}' conflicts with the selected CHARMM "
- f"force-field bundle ({generated} != {existing}). Use a CGenFF "
- "installation matching the registered charmm36m bundle."
- )
+ yield prm, stripped.split()[0], stripped
+
+
+def validate_cgenff_atomtype_compatibility(
+ prm_files: list[str | Path],
+ forcefield_types: dict[str, tuple[str, ...]],
+) -> None:
+ """Reject CGenFF atom types that conflict with the selected CHARMM bundle.
+
+ *forcefield_types* is the result of :func:`collect_forcefield_atomtype_definitions`
+ for the bundle.
+ """
+ for _, name, line in _iter_prm_atomtypes(prm_files):
+ existing = forcefield_types.get(name)
+ if existing is None:
+ continue
+ generated = _atomtype_signature(line)
+ if generated != existing:
+ raise ValueError(
+ f"CGenFF atom type '{name}' conflicts with the selected CHARMM "
+ f"force-field bundle ({generated} != {existing}). Use a CGenFF "
+ "installation matching the registered charmm36m bundle."
+ )
DEFAULTS_FUDGE_TOLERANCE = 1e-6
@@ -816,33 +844,27 @@ def read_topology_defaults(path: str | Path) -> tuple[int, int, str, float, floa
"""
in_defaults = False
- with open(path) as fb:
- for raw in fb:
- line = raw.split(";", 1)[0].strip()
- if not line or line.startswith("#"):
- continue
- if line.startswith("[") and line.endswith("]"):
- if in_defaults:
- break
- in_defaults = line[1:-1].strip().lower() == "defaults"
- continue
- if not in_defaults:
- continue
- fields = line.split()
- if len(fields) < 5:
- raise ValueError(f"Malformed [ defaults ] entry in {path}: {raw.strip()!r}")
- try:
- return (
- int(fields[0]),
- int(fields[1]),
- fields[2].lower(),
- float(fields[3]),
- float(fields[4]),
- )
- except ValueError as exc:
- raise ValueError(
- f"Malformed [ defaults ] entry in {path}: {raw.strip()!r}"
- ) from exc
+ for entry in iter_topology_entries(Path(path).read_text().splitlines(), headers=True):
+ if entry.is_header:
+ if in_defaults:
+ break
+ in_defaults = entry.section == "defaults"
+ continue
+ if not in_defaults:
+ continue
+ fields = entry.code.split()
+ if len(fields) < 5:
+ raise ValueError(f"Malformed [ defaults ] entry in {path}: {entry.code!r}")
+ try:
+ return (
+ int(fields[0]),
+ int(fields[1]),
+ fields[2].lower(),
+ float(fields[3]),
+ float(fields[4]),
+ )
+ except ValueError as exc:
+ raise ValueError(f"Malformed [ defaults ] entry in {path}: {entry.code!r}") from exc
if in_defaults:
raise ValueError(f"[ defaults ] section in {path} has no entry.")
return None
@@ -851,12 +873,11 @@ def read_topology_defaults(path: str | Path) -> tuple[int, int, str, float, floa
def check_defaults_compatible(
reference_path: str | Path,
other_path: str | Path,
- tolerance: float = DEFAULTS_FUDGE_TOLERANCE,
) -> None:
"""Raise unless two files declare equivalent ``[ defaults ]`` entries.
``nbfunc``, ``comb-rule`` and ``gen-pairs`` must match exactly; ``fudgeLJ``
- and ``fudgeQQ`` must agree within *tolerance*, because force fields write
+ and ``fudgeQQ`` must agree within :data:`DEFAULTS_FUDGE_TOLERANCE`, because force fields write
the 1-4 scaling with different precision (``0.833333`` versus
``0.83333333333333333``).
@@ -866,8 +887,6 @@ def check_defaults_compatible(
File whose ``[ defaults ]`` is kept, e.g. the force field's ``forcefield.itp``.
other_path : str or Path
File whose ``[ defaults ]`` would be dropped, e.g. a SMIRNOFF parameter file.
- tolerance : float, optional
- Absolute tolerance for the fudge factors. Default is 1e-6.
Raises
------
@@ -881,7 +900,9 @@ def check_defaults_compatible(
missing = reference_path if reference is None else other_path
raise ValueError(f"No [ defaults ] section found in {missing}.")
exact_match = reference[:3] == other[:3]
- fudge_match = all(abs(a - b) <= tolerance for a, b in zip(reference[3:], other[3:]))
+ fudge_match = all(
+ abs(a - b) <= DEFAULTS_FUDGE_TOLERANCE for a, b in zip(reference[3:], other[3:])
+ )
if not (exact_match and fudge_match):
raise ValueError(
"Incompatible [ defaults ] (nbfunc, comb-rule, gen-pairs, fudgeLJ, fudgeQQ): "
@@ -916,29 +937,20 @@ def check_atomtype_collisions(
"""
seen: dict[str, tuple[tuple[str, ...], str]] = {}
- for prm in prm_files:
- *_, sections = extract_reusable_parts_from_cgenff_gmx_top(prm, must_have_moleculetype=False)
- for section_name, section_text in sections:
- if section_name != "atomtypes":
- continue
- for line in section_text.splitlines():
- stripped = line.strip()
- if not stripped or stripped.startswith((";", "[")):
- continue
- name = stripped.split()[0]
- if name in forcefield_types:
- raise ValueError(
- f"Atom type '{name}' from {prm} is already defined by the force "
- "field; the parameter sets cannot be combined."
- )
- signature = _atomtype_signature(stripped)
- previous = seen.get(name)
- if previous is not None and previous[0] != signature:
- raise ValueError(
- f"Atom type '{name}' is defined with different parameters in "
- f"{previous[1]} and {prm}; the parameter sets cannot be combined."
- )
- seen.setdefault(name, (signature, str(prm)))
+ for prm, name, line in _iter_prm_atomtypes(prm_files):
+ if name in forcefield_types:
+ raise ValueError(
+ f"Atom type '{name}' from {prm} is already defined by the force "
+ "field; the parameter sets cannot be combined."
+ )
+ signature = _atomtype_signature(line)
+ previous = seen.get(name)
+ if previous is not None and previous[0] != signature:
+ raise ValueError(
+ f"Atom type '{name}' is defined with different parameters in "
+ f"{previous[1]} and {prm}; the parameter sets cannot be combined."
+ )
+ seen.setdefault(name, (signature, str(prm)))
def merge_extra_parameter_itps(
diff --git a/mdfactory/utils/units.py b/mdfactory/utils/units.py
new file mode 100644
index 0000000..9db8236
--- /dev/null
+++ b/mdfactory/utils/units.py
@@ -0,0 +1,12 @@
+# ABOUTME: Physical constants and unit conversions shared by the build code
+# ABOUTME: Leaf module (imports nothing from mdfactory) so models and utils can both use it
+"""Physical constants and unit conversions shared by the build code."""
+
+#: Avogadro's number, matching the g <-> Dalton convention used throughout.
+N_AVOGADRO = 6.022e23
+#: ų -> liter (1 ų = 1e-27 L), for molar concentrations.
+A3_TO_L = 1e-27
+#: ų -> cm³ (1 ų = 1e-24 cm³), for densities in g/cm³.
+A3_TO_CM3 = 1e-24
+#: cm³ -> ų.
+CM3_TO_A3 = 1e24
diff --git a/mdfactory/workflows.py b/mdfactory/workflows.py
index 4935af2..e6839f6 100644
--- a/mdfactory/workflows.py
+++ b/mdfactory/workflows.py
@@ -13,8 +13,6 @@
build_proteinbox,
)
from .models.input import BuildInput, write_build_record
-from .prepare import resolve_protein_pdb_path
-from .utils.utilities import load_yaml_file
DISPATCH_BUILD = {
"mixedbox": build_mixedbox,
@@ -51,20 +49,3 @@ def run_build_from_dict(inp_dict: dict[Any, Any] | BuildInput):
)
write_build_record(build_model, Path.cwd())
build_function(build_model)
-
-
-def run_build_from_file(fname: Path):
- """Load a YAML build specification and dispatch the build.
-
- Relative protein PDB paths are resolved against the YAML file's directory so
- the build finds them regardless of the process working directory.
-
- Parameters
- ----------
- fname : Path
- Path to the YAML build specification file
-
- """
- dct = load_yaml_file(fname)
- resolve_protein_pdb_path(dct, Path(fname).parent)
- run_build_from_dict(dct)