Allow the user to supply starting conformers from an SDF

docs/changelog.md

@@ -4,11 +4,17 @@
 
 ### Features
 
+- Add an optional `starting_conformers` setting to each sampling stage (`training_sampling_settings`, `testing_sampling_settings`) and to `msm_settings`. When set to an SDF path, that stage starts from the supplied conformers (matched to each molecule by graph isomorphism and realigned automatically) instead of generating them with ETKDG; `n_conformers` is ignored for that stage. The default remains ETKDG.In [#78](https://github.com/cole-group/presto/pull/78).
 - Add `presto.create_types.add_library_charges_to_forcefield` to write custom partial charges from OpenFF `Molecule` objects into a force field as library charges, addressing [#64](https://github.com/cole-group/presto/issues/64).
 
+### Fixes
+
+- Fix a latent `IndexError` in the MSM step when fewer conformers were available than `n_conformers`; the conformer loop now iterates the conformers actually present.
+
 ### Documentation
 
 - Add [Use custom charges](how-to/use-custom-charges.md) how-to guide.
+- Add [Use your own starting conformers](how-to/use-starting-conformers.md) how-to guide.
 - Add a `CITATION.cff` file and cite the presto preprint in the README and docs, and point users to the OpenFF publications to cite, in [#76](https://github.com/cole-group/presto/pull/76).
 - Document installing `presto` from `conda-forge` as `presto-fit` (note this comes without the MLP dependencies, which must be installed separately) in [#70](https://github.com/cole-group/presto/pull/70).
 

docs/how-to/index.md

@@ -7,6 +7,7 @@ Short, task-oriented recipes. Each page assumes you've finished **[Get started](
 - **[Fit a single molecule](fit-single-molecule.md)** — defaults you may want to change.
 - **[Fit a congeneric series](fit-congeneric-series.md)** — share parameters across related molecules using `max_extend_distance`.
 - **[Use SDF inputs](use-sdf-inputs.md)** — switch from SMILES to one or more `.sdf` files.
+- **[Use your own starting conformers](use-starting-conformers.md)** — seed a sampling stage (or MSM) from an SDF instead of ETKDG.
 - **[Use custom charges](use-custom-charges.md)** — bake your own partial charges into the force field as library charges.
 - **[Wipe output and rerun](clean-rerun.md)**
 

docs/how-to/use-starting-conformers.md

@@ -0,0 +1,69 @@
+# Use your own starting conformers
+
+By default each sampling stage generates its starting conformers with RDKit ETKDG. You can
+instead start from your own conformers (for example a curated ensemble, docked or crystal
+poses, or a torsion scan) by pointing a stage at an SDF file with `starting_conformers`.
+
+This is available **per stage**, so you can mix supplied conformers and ETKDG freely:
+
+- `training_sampling_settings.starting_conformers` — training-data MD
+- `testing_sampling_settings.starting_conformers` — test-data MD
+- `param_settings.msm_settings.starting_conformers` — the MSM Hessian step that seeds the
+  initial bond/angle parameters
+
+Any stage left unset keeps generating conformers with ETKDG.
+
+## Behaviour
+
+- **The SDF takes precedence.** When `starting_conformers` is set for a stage, that stage
+  starts from **every** conformer in the file, and the stage's `n_conformers` is **ignored**
+  (a log line makes this explicit).
+- **Conformers are matched by connectivity.** Each record is matched to the molecule being
+  fitted by graph isomorphism, and its atom ordering is **remapped automatically** to match.
+  You do not need to worry about atom ordering, and a single SDF can hold conformers for
+  several molecules in a multi-molecule fit — each molecule picks up only its own records.
+- **Missing molecules fail fast.** If a configured SDF contains no conformer for a molecule
+  being fitted, the run stops during settings validation, before any expensive work.
+
+## CLI form
+
+```bash
+presto train \
+    --param-settings.molecules "CCO" \
+    --training-sampling-settings.starting-conformers my_conformers.sdf
+```
+
+## YAML form
+
+```yaml
+param_settings:
+    molecule_input_type: smiles
+    molecules:
+        - CCO
+    msm_settings:
+        # optional: also seed the MSM Hessian step from the same (or a different) SDF
+        starting_conformers: my_conformers.sdf
+
+training_sampling_settings:
+    sampling_protocol: mm_md_metadynamics_torsion_minimisation
+    starting_conformers: my_conformers.sdf
+
+testing_sampling_settings:
+    sampling_protocol: ml_md
+    # left unset -> ETKDG
+```
+
+## Preparing the SDF
+
+The SDF should contain one record per conformer, all of the same molecule (multiple
+molecules are fine — group them however you like, they are matched by connectivity). For
+example, from an existing OpenFF molecule:
+
+```python
+from rdkit import Chem
+
+rdkit_molecule = molecule.to_rdkit()  # molecule with N conformers
+with Chem.SDWriter("my_conformers.sdf") as writer:
+    for conformer_id in range(rdkit_molecule.GetNumConformers()):
+        writer.write(rdkit_molecule, confId=conformer_id)
+```

presto/load_molecules.py

@@ -5,8 +5,11 @@
 from typing import Literal
 
 from openff.toolkit import Molecule
+from openff.units import Quantity
 from rdkit import Chem
 
+from .utils.typing import PathLike
+
 MoleculeInputType = Literal["smiles", "sdf"]
 
 MoleculeLoader = Callable[[str], list[Molecule]]
@@ -72,3 +75,81 @@ def load_sdf_molecules(input_value: str) -> list[Molecule]:
     "smiles": load_smiles_molecules,
     "sdf": load_sdf_molecules,
 }
+
+
+def load_conformers_for_molecule(
+    molecule: Molecule, sdf_path: PathLike
+) -> list[Quantity]:
+    """Load the conformers of ``molecule`` from an SDF, aligned to its atom ordering.
+
+    Every record in the SDF that is graph-isomorphic to ``molecule`` is treated as a
+    conformer of it. Because the atom ordering in the SDF need not match ``molecule``,
+    each matching record is remapped onto ``molecule``'s atom ordering before its
+    coordinates are extracted, guaranteeing the returned conformers are valid starting
+    positions for a topology built from ``molecule``. Records that do not match are
+    ignored (they may belong to another molecule in a multi-molecule SDF).
+
+    Parameters
+    ----------
+    molecule : openff.toolkit.Molecule
+        The molecule whose conformers should be loaded. Defines the canonical atom
+        ordering the returned conformers are aligned to.
+    sdf_path : PathLike
+        Path to an SDF file containing one or more conformers of ``molecule`` (and,
+        optionally, of other molecules).
+
+    Returns:
+    -------
+    list[openff.units.Quantity]
+        The matching conformers, each aligned to ``molecule``'s atom ordering.
+
+    Raises:
+    ------
+    ValueError
+        If the path does not exist, does not end in ``.sdf``, or contains no record
+        matching ``molecule``.
+    """
+    path = Path(sdf_path)
+
+    if not path.exists():
+        raise ValueError(f"SDF file does not exist: {path}")
+
+    if path.suffix.lower() != ".sdf":
+        raise ValueError(f"Expected an SDF file path ending in .sdf: {path}")
+
+    try:
+        supplier = Chem.SDMolSupplier(str(path), removeHs=False)
+    except Exception as exc:
+        raise ValueError(f"Failed to read SDF file: {path}") from exc
+
+    conformers: list[Quantity] = []
+
+    for rdkit_molecule in supplier:
+        if rdkit_molecule is None:
+            continue
+
+        try:
+            record = Molecule.from_rdkit(rdkit_molecule, allow_undefined_stereo=True)
+        except Exception:
+            # Records that cannot be interpreted as molecules cannot be a conformer of
+            # ``molecule``, so skip them rather than failing the whole load.
+            continue
+
+        matched, atom_map = Molecule.are_isomorphic(
+            record, molecule, return_atom_map=True
+        )
+        if not matched or atom_map is None:
+            continue
+
+        # ``atom_map`` maps record atom indices -> ``molecule`` atom indices, so remapping
+        # the record with ``current_to_new`` yields a molecule in ``molecule``'s ordering.
+        aligned = record.remap(atom_map, current_to_new=True)
+        conformers.extend(aligned.conformers)
+
+    if not conformers:
+        raise ValueError(
+            f"SDF file {path} contains no conformers matching the molecule "
+            f"{_molecule_identity(molecule.to_rdkit())}."
+        )
+
+    return conformers

presto/msm.py

@@ -698,9 +698,11 @@ def apply_msm_to_molecule(
         Tuple of (bond_params_dict, angle_params_dict). If multiple conformers are used,
         the parameters are averaged over all conformers.
     """
-    # Generate conformers
+    # Generate conformers (or load the supplied starting conformers)
     mol_with_conformers = _copy_mol_and_add_conformers(
-        mol, n_conformers=settings.n_conformers
+        mol,
+        n_conformers=settings.n_conformers,
+        starting_conformers=settings.starting_conformers,
     )
     simulation, integrator = _build_ml_simulation(
         mol_with_conformers,
@@ -717,14 +719,12 @@ def apply_msm_to_molecule(
         defaultdict(list)
     )
 
-    for conf_idx in track(
-        range(settings.n_conformers),
+    for conformer in track(
+        mol_with_conformers.conformers,
         description="Finding MSM parameters for conformers",
     ):
         # Set positions for this conformer
-        simulation.context.setPositions(
-            mol_with_conformers.conformers[conf_idx].to_openmm()
-        )
+        simulation.context.setPositions(conformer.to_openmm())
 
         # Minimize to local energy minimum
         simulation.minimizeEnergy(maxIterations=0, tolerance=settings.tolerance)

presto/sample.py

@@ -33,6 +33,7 @@
     DEFAULT_TORSIONS_TO_INCLUDE_SMARTS,
     get_rot_torsions_by_rot_bond,
 )
+from .load_molecules import load_conformers_for_molecule
 from .metadynamics import Metadynamics
 from .outputs import OutputType, get_mol_path
 from .utils.gpu import cleanup_simulation
@@ -77,16 +78,36 @@ def __call__(self, **kwargs: Unpack[SampleFnArgs]) -> list[datasets.Dataset]:
 def _copy_mol_and_add_conformers(
     mol: openff.toolkit.Molecule,
     n_conformers: int,
+    starting_conformers: pathlib.Path | None = None,
 ) -> openff.toolkit.Molecule:
-    """Copy a molecule and add conformers to it."""
+    """Copy a molecule and add starting conformers to it.
+
+    If ``starting_conformers`` is None (default), ``n_conformers`` conformers are
+    generated with ETKDG. Otherwise the conformers are loaded from the given SDF (matched
+    to ``mol`` by graph and aligned to its atom ordering) and ``n_conformers`` is ignored.
+    """
     mol = copy.deepcopy(mol)
-    mol.generate_conformers(n_conformers=n_conformers, rms_cutoff=0.0 * _ANGSTROM)
-    n_gen_conformers = len(mol.conformers)
-    if n_gen_conformers < n_conformers:
-        logger.warning(
-            f"Only {n_gen_conformers} conformers were generated, which is less than the requested {n_conformers}."
-            f" As a result, {n_gen_conformers / n_conformers * 100:.1f}% of the requested samples will be generated."
-        )
+
+    if starting_conformers is None:
+        mol.generate_conformers(n_conformers=n_conformers, rms_cutoff=0.0 * _ANGSTROM)
+        n_gen_conformers = len(mol.conformers)
+        if n_gen_conformers < n_conformers:
+            logger.warning(
+                f"Only {n_gen_conformers} conformers were generated, which is less than the requested {n_conformers}."
+                f" As a result, {n_gen_conformers / n_conformers * 100:.1f}% of the requested samples will be generated."
+            )
+        return mol
+
+    # Drop any incidental conformer (e.g. one carried by an SDF-loaded molecule) before
+    # attaching the supplied starting conformers.
+    mol._conformers = []
+    conformers = load_conformers_for_molecule(mol, starting_conformers)
+    logger.info(
+        f"Starting from {len(conformers)} supplied conformers in {starting_conformers} "
+        f"(n_conformers={n_conformers} ignored)."
+    )
+    for conformer in conformers:
+        mol.add_conformer(conformer)
     return mol
 
 
@@ -338,7 +359,9 @@ def sample_mmmd(
     all_datasets = []
 
     for mol_idx, mol in enumerate(mols):
-        mol_with_conformers = _copy_mol_and_add_conformers(mol, settings.n_conformers)
+        mol_with_conformers = _copy_mol_and_add_conformers(
+            mol, settings.n_conformers, settings.starting_conformers
+        )
         interchange = openff.interchange.Interchange.from_smirnoff(
             off_ff, openff.toolkit.Topology.from_molecules(mol_with_conformers)
         )
@@ -433,7 +456,9 @@ def sample_mlmd(
     all_datasets = []
 
     for mol_idx, mol in enumerate(mols):
-        mol_with_conformers = _copy_mol_and_add_conformers(mol, settings.n_conformers)
+        mol_with_conformers = _copy_mol_and_add_conformers(
+            mol, settings.n_conformers, settings.starting_conformers
+        )
         ml_simulation, integrator = _build_ml_simulation(
             mol_with_conformers,
             mol_with_conformers.to_topology().to_openmm(),
@@ -559,7 +584,9 @@ def sample_mmmd_metadynamics(
     all_datasets = []
 
     for mol_idx, mol in enumerate(mols):
-        mol_with_conformers = _copy_mol_and_add_conformers(mol, settings.n_conformers)
+        mol_with_conformers = _copy_mol_and_add_conformers(
+            mol, settings.n_conformers, settings.starting_conformers
+        )
         interchange = openff.interchange.Interchange.from_smirnoff(
             off_ff, openff.toolkit.Topology.from_molecules(mol_with_conformers)
         )
@@ -1175,7 +1202,9 @@ def sample_mmmd_metadynamics_with_torsion_minimisation(
     all_datasets = []
 
     for mol_idx, mol in enumerate(mols):
-        mol_with_conformers = _copy_mol_and_add_conformers(mol, settings.n_conformers)
+        mol_with_conformers = _copy_mol_and_add_conformers(
+            mol, settings.n_conformers, settings.starting_conformers
+        )
         interchange = openff.interchange.Interchange.from_smirnoff(
             off_ff, openff.toolkit.Topology.from_molecules(mol_with_conformers)
         )