From bb5a58b9b52bd6be67c2d7d6c355deaace50c34b Mon Sep 17 00:00:00 2001 From: Corin Wagen Date: Fri, 3 Apr 2026 14:59:07 -0400 Subject: [PATCH 1/2] add constraints + analogue mode --- README.md | 95 ++++++++++-- SCIENCE.md | 32 ++++ openconf/__init__.py | 8 +- openconf/api.py | 91 ++++++++++- openconf/config.py | 55 ++++++- openconf/perceive.py | 89 +++++++++++ openconf/propose/hybrid.py | 243 +++++++++++++++++++++++++++-- pyproject.toml | 2 +- tests/test_constrained.py | 307 +++++++++++++++++++++++++++++++++++++ 9 files changed, 890 insertions(+), 32 deletions(-) create mode 100644 tests/test_constrained.py diff --git a/README.md b/README.md index 440ef5f..49bbf3e 100644 --- a/README.md +++ b/README.md @@ -40,7 +40,7 @@ ensemble.to_xyz("output.xyz") ### Named Presets -Four use-case presets are available out of the box: +Five use-case presets are available out of the box: ```python from openconf import generate_conformers @@ -51,6 +51,8 @@ ensemble = generate_conformers(mol, preset="spectroscopic") # NMR / IR / VCD ensemble = generate_conformers(mol, preset="docking") # docking pose recovery ``` +For FEP-style analogue generation from a fixed pose, see [`generate_conformers_from_pose`](#5-analogue--fep-style-r-group-exploration) below. + ### Custom Configuration For full control, pass a `ConformerConfig` directly. You can also use @@ -106,7 +108,7 @@ config.max_out = 200 # override a single field ensemble = generate_conformers(mol, config=config) ``` -Available presets: `"rapid"`, `"ensemble"`, `"spectroscopic"`, `"docking"`. +Available presets: `"rapid"`, `"ensemble"`, `"spectroscopic"`, `"docking"`, `"analogue"`. Below are representative wall-clock timings measured on a single CPU core (Apple M2 Pro), mean over 3 runs. @@ -325,18 +327,85 @@ ensemble.to_sdf("docking_input.sdf") --- +### 5. Analogue / FEP-style R-group exploration + +Generate conformers for an MCS-aligned analogue while keeping the core scaffold +exactly fixed at the input pose. The correct entry point here is +`generate_conformers_from_pose` rather than `generate_conformers`. + +- Starts from the **supplied conformer** — no ETKDG seeding +- Only **free terminal rotors** are explored (those whose moving fragment is + entirely outside the constrained core) +- MMFF minimization uses stiff **position restraints** on all constrained atoms, + then snaps them to exact starting coordinates so there is zero drift +- **Global shake** is suppressed to avoid thrashing the starting pose + +```python +from rdkit import Chem +from rdkit.Chem import AllChem +from openconf import generate_conformers_from_pose + +# Suppose we have an MCS-aligned analogue with a propyl substituent +# replacing the butyl chain on a benzene scaffold. +mol = Chem.MolFromSmiles("CCCc1ccccc1") +mol = Chem.AddHs(mol) +AllChem.EmbedMolecule(mol, AllChem.ETKDGv3()) + +# Ring heavy-atom indices — these must not move +ring_atoms = [3, 4, 5, 6, 7, 8] + +ensemble = generate_conformers_from_pose(mol, constrained_atoms=ring_atoms) +ensemble.to_sdf("analogues.sdf") +``` + +The default preset (`"analogue"`) returns up to 50 conformers. Pass `preset=` or +`config=` to override: + +```python +from openconf import ConformerConfig, generate_conformers_from_pose + +# Fewer conformers, faster turnaround +config = ConformerConfig(max_out=10, n_steps=60, pool_max=200) +ensemble = generate_conformers_from_pose(mol, constrained_atoms=ring_atoms, config=config) +``` + +
Full analogue preset equivalent + +```python +from openconf import ConformerConfig, PrismConfig, generate_conformers_from_pose + +config = ConformerConfig( + max_out=50, + pool_max=500, + n_steps=150, + energy_window_kcal=10.0, + do_final_refine=True, + minimize_batch_size=8, + parent_strategy="softmax", + final_select="diverse", + prism_config=PrismConfig(energy_window_kcal=10.0), +) +ensemble = generate_conformers_from_pose(mol, constrained_atoms=ring_atoms, config=config) +``` + +
+ +--- + ### Configuration Quick Reference -| Parameter | Rapid | Ensemble | Spectroscopic | Docking | -|---|---|---|---|---| -| `max_out` | 5 | 50 | 100 | 250 | -| `n_steps` | 30 | 200 | 400 | 500 | -| `energy_window_kcal` | 20 | 10 | 5 | 18 | -| `seed_n_per_rotor` | 2 | 3 | 5 | 4 | -| `seed_prune_rms_thresh` | 1.5 | 1.0 | 0.5 | 0.8 | -| `do_final_refine` | False | True | True | False | -| `parent_strategy` | softmax | softmax | softmax | uniform | -| `final_select` | diverse | diverse | energy | diverse | +| Parameter | Rapid | Ensemble | Spectroscopic | Docking | Analogue | +|---|---|---|---|---|---| +| `max_out` | 5 | 50 | 100 | 250 | 50 | +| `n_steps` | 30 | 200 | 400 | 500 | 150 | +| `energy_window_kcal` | 20 | 10 | 5 | 18 | 10 | +| `seed_n_per_rotor` | 2 | 3 | 5 | 4 | — | +| `seed_prune_rms_thresh` | 1.5 | 1.0 | 0.5 | 0.8 | — | +| `do_final_refine` | False | True | True | False | True | +| `parent_strategy` | softmax | softmax | softmax | uniform | softmax | +| `final_select` | diverse | diverse | energy | diverse | diverse | + +*Analogue mode uses `generate_conformers_from_pose`; seeding parameters are unused because ETKDG is skipped.* ## How It Works @@ -380,11 +449,13 @@ openconf is not recommended for macrocycles (ring size ≥ 12). On macrocyclic r - `generate_conformers(mol, method="hybrid", config=None)` - Main entry point - `generate_conformers_from_smiles(smiles, ...)` - Convenience wrapper +- `generate_conformers_from_pose(mol, constrained_atoms, config=None)` - FEP-style analogue generation from an aligned pose ### Configuration Classes - `ConformerConfig` - Main configuration - `PrismConfig` - PRISM Pruner settings +- `ConstraintSpec` - Positional constraints for pose-locked generation ### Data Classes diff --git a/SCIENCE.md b/SCIENCE.md index f7608cf..dbcd9f8 100644 --- a/SCIENCE.md +++ b/SCIENCE.md @@ -12,6 +12,38 @@ Seed count is computed automatically from the molecular topology rather than set The key to making this efficient is aggressive deduplication. Without it, the search quickly fills with near-identical structures that differ only in insignificant ways. openconf uses PRISM Pruner, which implements a cached divide-and-conquer algorithm for comparing conformers by RMSD and moment of inertia. Rather than doing O(N²) all-to-all comparisons, PRISM sorts conformers by energy and recursively partitions them, exploiting the fact that similar structures tend to cluster. This lets openconf maintain large internal pools (thousands of conformers) while keeping only the truly unique ones. The final selection step returns the lowest-energy conformers after PRISM deduplication, ensuring the output ensemble contains distinct, low-strain geometries. +## Pose-Constrained Generation (FEP / Analogue Mode) + +A common need in lead optimization is generating conformers for an analogue where the core scaffold is already placed — for example, the result of an MCS alignment from a co-crystal structure. Standard ETKDG seeding is not appropriate here: it randomizes the entire molecule and cannot be guaranteed to recover the aligned pose. Instead, openconf supports constrained generation via `generate_conformers_from_pose`. + +The idea is to identify a set of *constrained atoms* (the MCS core) whose Cartesian positions must remain fixed at the input geometry, and to explore only the remaining degrees of freedom — the *free rotors*, i.e., bonds whose entire moving fragment lies outside the constrained set. The algorithm adapts in three ways: + +**Seeding.** ETKDG is skipped entirely. The single input conformer is used directly as the seed, immediately fast-minimized with position restraints to relax any bond-length or angle strain in the free fragment before exploration begins. + +**Move set.** Only free rotors (and ring flips of entirely free rings) are sampled. The global shake move is suppressed: it would randomize 50–80% of all rotors, some of which pass through the core, defeating the constraint. Move probability that would have gone to global shake is redistributed to `single_rotor` moves. + +**Minimization.** Every MMFF minimization call applies `MMFFAddPositionConstraint` to each constrained atom with a stiff harmonic force constant (default 1000 kcal/mol/Ų). After minimization converges, constrained atom coordinates are snapped back to the exact reference values, eliminating any residual drift. Final refinement applies the same treatment. + +The rotor filtering logic uses a BFS traversal: for each bond (i, j) in the rotor list, we compute the set of atoms reachable from j without crossing (i, j). If this moving fragment contains no constrained atoms, the rotor is free. If the moving fragment *does* contain constrained atoms but the opposite side (reachable from i) does not, the rotor is flipped so that the free atoms become the moving side. Rotors where constrained atoms appear on both sides are excluded entirely — they cannot be rotated without moving the core. + +```python +from rdkit import Chem +from rdkit.Chem import AllChem +from openconf import generate_conformers_from_pose + +# MCS-aligned analogue: ethyl group added to a benzene scaffold +mol = Chem.MolFromSmiles("CCc1ccccc1") +mol = Chem.AddHs(mol) +AllChem.EmbedMolecule(mol, AllChem.ETKDGv3()) # or use your aligned pose + +# Constrain the benzene ring (heavy-atom indices 2–7) +ensemble = generate_conformers_from_pose(mol, constrained_atoms=list(range(2, 8))) +``` + +The `"analogue"` preset (50 conformers, 150 steps, softmax parent strategy, full refinement) is the default. Because the free rotor space is much smaller than the full conformational space, 150 steps is usually sufficient for thorough coverage of simple R-groups. For larger substituents with many free rotors, increase `n_steps` accordingly. + +**Atom index convention.** `Chem.AddHs` appends new H atoms after all existing atoms, preserving all prior indices. So whether you pass a heavy-atom-only mol (indices 0…N−1 for heavy atoms) or an H-added mol (same heavy-atom indices plus new H indices), the constrained atom indices remain valid after the internal `AddHs` call. + ## Tuning Guide For the most common workflows, named presets are the easiest starting point: diff --git a/openconf/__init__.py b/openconf/__init__.py index d892bf8..fa9ac79 100644 --- a/openconf/__init__.py +++ b/openconf/__init__.py @@ -18,12 +18,13 @@ ConformerEnsemble, ConformerRecord, generate_conformers, + generate_conformers_from_pose, generate_conformers_from_smiles, ) -from .config import ConformerConfig, ConformerPreset, PrismConfig, preset_config +from .config import ConformerConfig, ConformerPreset, ConstraintSpec, PrismConfig, preset_config from .dedupe import prism_dedupe from .io import mol_to_smiles, read_sdf, smiles_to_mol, write_sdf, write_xyz -from .perceive import Rotor, RotorModel, build_rotor_model, prepare_molecule +from .perceive import Rotor, RotorModel, build_rotor_model, filter_constrained_rotors, prepare_molecule from .relax import RDKitMMFFMinimizer, get_minimizer from .torsionlib import TorsionLibrary, TorsionRule @@ -32,6 +33,7 @@ "ConformerEnsemble", "ConformerPreset", "ConformerRecord", + "ConstraintSpec", "PrismConfig", "RDKitMMFFMinimizer", "Rotor", @@ -39,7 +41,9 @@ "TorsionLibrary", "TorsionRule", "build_rotor_model", + "filter_constrained_rotors", "generate_conformers", + "generate_conformers_from_pose", "generate_conformers_from_smiles", "get_minimizer", "mol_to_smiles", diff --git a/openconf/api.py b/openconf/api.py index deb395a..55f831f 100644 --- a/openconf/api.py +++ b/openconf/api.py @@ -4,7 +4,7 @@ from rdkit import Chem -from .config import ConformerConfig, ConformerPreset, preset_config +from .config import ConformerConfig, ConformerPreset, ConstraintSpec, preset_config from .perceive import build_rotor_model, prepare_molecule from .pool import ConformerRecord from .propose.hybrid import run_hybrid_generation @@ -189,6 +189,95 @@ def generate_conformers( return ConformerEnsemble(mol=mol, records=records) +def generate_conformers_from_pose( + mol: Chem.Mol, + constrained_atoms: list[int] | frozenset[int], + config: ConformerConfig | None = None, + preset: ConformerPreset | None = None, +) -> ConformerEnsemble: + """Generate conformers for an MCS-aligned pose, keeping core atoms fixed. + + Designed for FEP-style analogue generation: you supply a molecule with an + existing conformer (e.g. the result of MCS alignment) and the indices of + the atoms that must remain fixed (the MCS core / scaffold). Only the free + terminal rotors are explored, so the search is faster and stays consistent + with the bound pose. + + **Atom index convention:** pass indices as they appear in *mol*. If *mol* + already has explicit hydrogens the indices are used as-is. If *mol* has only + heavy atoms, ``Chem.AddHs`` is called internally — it appends new H atoms at + the end and leaves all existing indices unchanged, so heavy-atom indices + remain valid after H addition. + + Args: + mol: RDKit molecule with at least one conformer (the MCS-aligned pose). + If multiple conformers are present, the first is used as the seed. + constrained_atoms: Atom indices of the core scaffold that must not move. + These are indices into *mol* as supplied (see note above). + config: Configuration options. Defaults to the ``"analogue"`` preset. + preset: Named preset. Defaults to ``"analogue"`` when neither *config* + nor *preset* is given. Mutually exclusive with *config*. + + Returns: + ConformerEnsemble with terminal-group conformational diversity while + preserving the input core geometry. + + Raises: + ValueError: If mol has no conformers, if both *config* and *preset* are + supplied, or if the method is unknown. + + Examples: + >>> from rdkit import Chem + >>> from openconf import generate_conformers_from_pose + >>> mol = Chem.MolFromSmiles("CCCc1ccccc1") + >>> # (assume mol already has an aligned conformer) + >>> ensemble = generate_conformers_from_pose( # doctest: +SKIP + ... mol, constrained_atoms=[4, 5, 6, 7, 8, 9] + ... ) + """ + if config is not None and preset is not None: + raise ValueError("Specify at most one of 'config' or 'preset', not both.") + + if mol.GetNumConformers() == 0: + raise ValueError( + "mol must have at least one conformer for pose-constrained generation. " + "Supply the MCS-aligned pose as conformer 0." + ) + + # Resolve config — default to "analogue" preset for this entry point. + if config is not None: + resolved_config = ConformerConfig(**{ + k: v for k, v in config.__dict__.items() if k != "constraint_spec" + }) + elif preset is not None: + resolved_config = preset_config(preset) + else: + resolved_config = preset_config("analogue") + + # Attach the constraint spec (overrides any constraint_spec already in config). + resolved_config.constraint_spec = ConstraintSpec( + constrained_atoms=frozenset(constrained_atoms) + ) + + # Prepare molecule — AddHs is a no-op if Hs are already present. + prepped_mol = prepare_molecule(mol, add_hs=True) + + rotor_model = build_rotor_model(prepped_mol) + + prepped_mol, conf_ids, energies = run_hybrid_generation(prepped_mol, rotor_model, resolved_config) + + records = [ + ConformerRecord( + conf_id=cid, + energy_kcal=energy, + source="analogue", + ) + for cid, energy in zip(conf_ids, energies, strict=True) + ] + + return ConformerEnsemble(mol=prepped_mol, records=records) + + def generate_conformers_from_smiles( smiles: str, method: str = "hybrid", diff --git a/openconf/config.py b/openconf/config.py index 4a4fa43..ac69e66 100644 --- a/openconf/config.py +++ b/openconf/config.py @@ -3,7 +3,29 @@ from dataclasses import dataclass, field from typing import Literal -ConformerPreset = Literal["rapid", "ensemble", "spectroscopic", "docking"] +ConformerPreset = Literal["rapid", "ensemble", "spectroscopic", "docking", "analogue"] + + +@dataclass +class ConstraintSpec: + """Specification for positional constraints during conformer generation. + + Used for FEP-style analogue generation where an MCS-aligned pose is provided + and a subset of atoms (the core scaffold) must remain fixed while terminal + groups are explored. + + Attributes: + constrained_atoms: Atom indices that must not move. These are indices + into the molecule after hydrogen addition — since ``Chem.AddHs`` + preserves all existing atom indices, you can pass indices from either + the heavy-atom mol or the H-added mol interchangeably. + position_force_constant: MMFF force constant (kcal/mol/Ų) for the + harmonic position restraints applied to constrained atoms. + The default (1000.0) is very stiff and effectively freezes the core. + """ + + constrained_atoms: frozenset[int] + position_force_constant: float = 1000.0 @dataclass @@ -89,6 +111,13 @@ class ConformerConfig: refinement pass. Lower than fast_dielectric (default 4.0) to give more physically meaningful energies representative of a condensed-phase environment (protein interior / organic solvent). + constraint_spec: Optional positional constraints for FEP-style analogue + generation. When set, the specified atoms are pinned to their starting + coordinates via MMFF position restraints, ETKDG seeding is skipped + (the input conformer is used as the sole seed), global shake moves are + suppressed, and only rotors whose moving fragment is entirely outside + the constrained set are explored. Normally set via + ``generate_conformers_from_pose`` rather than directly. """ max_out: int = 200 @@ -113,6 +142,7 @@ class ConformerConfig: random_seed: int | None = None num_threads: int = 0 use_heavy_atoms_only: bool = True + constraint_spec: "ConstraintSpec | None" = None clash_threshold: float = 1.5 fast_minimization_iters: int = 20 max_minimization_iters: int = 200 @@ -147,9 +177,14 @@ def preset_config(preset: ConformerPreset) -> "ConformerConfig": workflows. Wide energy window, uniform parent sampling, no final refinement (docking programs minimize in-situ). + - ``"analogue"`` — FEP-style analogue / R-group enumeration. Intended for + use with ``generate_conformers_from_pose``, which supplies the + ``constraint_spec`` automatically. 50 conformers, full refinement, + softmax parent strategy to stay near the constrained energy basin. + Args: preset: One of ``"rapid"``, ``"ensemble"``, ``"spectroscopic"``, - ``"docking"``. + ``"docking"``, ``"analogue"``. Returns: ConformerConfig configured for the requested use case. @@ -219,7 +254,21 @@ def preset_config(preset: ConformerPreset) -> "ConformerConfig": final_select="diverse", prism_config=PrismConfig(energy_window_kcal=18.0), ) + case "analogue": + return ConformerConfig( + max_out=50, + pool_max=500, + n_steps=150, + energy_window_kcal=10.0, + seed_n_per_rotor=1, + seed_prune_rms_thresh=1.0, + do_final_refine=True, + minimize_batch_size=8, + parent_strategy="softmax", + final_select="diverse", + prism_config=PrismConfig(energy_window_kcal=10.0), + ) case _: raise ValueError( - f"Unknown preset {preset!r}. Choose from: 'rapid', 'ensemble', 'spectroscopic', 'docking'." + f"Unknown preset {preset!r}. Choose from: 'rapid', 'ensemble', 'spectroscopic', 'docking', 'analogue'." ) diff --git a/openconf/perceive.py b/openconf/perceive.py index 25c6fc0..12ae5c6 100644 --- a/openconf/perceive.py +++ b/openconf/perceive.py @@ -324,6 +324,95 @@ def _find_ring_flips(mol: Chem.Mol, atom_rings: list[tuple[int, ...]]) -> list[R return flips +def _atoms_on_side(mol: Chem.Mol, start_atom: int, excluded_atom: int) -> frozenset[int]: + """BFS from start_atom without crossing the bond to excluded_atom. + + Returns the set of atom indices reachable from start_atom when the bond + between excluded_atom and start_atom is severed. For acyclic bonds this + gives the atoms in one fragment; ring bonds are already excluded from the + rotor list so the result is always well-defined. + + Args: + mol: RDKit molecule. + start_atom: Index to start BFS from (one end of the bond). + excluded_atom: The other end of the bond — never visited. + + Returns: + frozenset of atom indices on start_atom's side of the bond. + """ + visited: set[int] = {start_atom} + queue = [start_atom] + while queue: + current = queue.pop() + for neighbor in mol.GetAtomWithIdx(current).GetNeighbors(): + nb_idx = neighbor.GetIdx() + if nb_idx not in visited and nb_idx != excluded_atom: + visited.add(nb_idx) + queue.append(nb_idx) + return frozenset(visited) + + +def filter_constrained_rotors(rotor_model: "RotorModel", constrained_atoms: frozenset[int]) -> "RotorModel": + """Return a new RotorModel containing only rotors that move no constrained atoms. + + A rotor around bond (i, j) is included if the fragment on one side of the + bond contains no constrained atoms (meaning we can rotate that fragment + freely). When the free fragment happens to be on the i-side rather than the + j-side as stored in dihedral_atoms, the rotor is flipped so the moving + fragment is always the free one. + + Ring flips are kept only when the entire ring is free of constrained atoms. + + Args: + rotor_model: Full rotor model built by build_rotor_model. + constrained_atoms: Atom indices that must not move. + + Returns: + New RotorModel with only free rotors and free ring flips. + """ + mol = rotor_model.mol + free_rotors: list[Rotor] = [] + + for rotor in rotor_model.rotors: + atom_i, atom_j = rotor.atom_idxs + moving_j = _atoms_on_side(mol, atom_j, atom_i) + + if not constrained_atoms & moving_j: + # j-side is free — use rotor as-is (SetDihedralDeg moves j-side) + free_rotors.append(rotor) + else: + moving_i = _atoms_on_side(mol, atom_i, atom_j) + if not constrained_atoms & moving_i: + # i-side is free — flip so the free side is the moving side + a, i, j, b = rotor.dihedral_atoms + free_rotors.append( + Rotor( + bond_idx=rotor.bond_idx, + atom_idxs=(atom_j, atom_i), + dihedral_atoms=(b, j, i, a), + rotor_type=rotor.rotor_type, + ) + ) + # else: constrained atoms on both sides — exclude this rotor + + # Ring flips: keep only rings with no constrained atoms + free_ring_flips = [ + rf for rf in rotor_model.ring_flips if not constrained_atoms & frozenset(rf.ring_atoms) + ] + + adj = _build_rotor_adjacency(free_rotors, mol) + + return RotorModel( + mol=mol, + rotors=free_rotors, + adj=adj, + ring_info=rotor_model.ring_info, + ring_flips=free_ring_flips, + heavy_atom_indices=rotor_model.heavy_atom_indices, + n_rotatable=len(free_rotors), + ) + + def build_rotor_model(mol: Chem.Mol) -> RotorModel: """Build a rotor model for a molecule. diff --git a/openconf/propose/hybrid.py b/openconf/propose/hybrid.py index 276a63b..77d7eeb 100644 --- a/openconf/propose/hybrid.py +++ b/openconf/propose/hybrid.py @@ -9,8 +9,8 @@ from rdkit import Chem from rdkit.Chem import AllChem, rdMolTransforms -from ..config import ConformerConfig -from ..perceive import RotorModel +from ..config import ConformerConfig, ConstraintSpec +from ..perceive import RotorModel, filter_constrained_rotors from ..pool import ConformerPool from ..relax import get_minimizer, minimize_confs_mmff from ..torsionlib import TorsionLibrary @@ -148,6 +148,7 @@ def __init__( rotor_model: RotorModel, torsion_lib: TorsionLibrary, config: ConformerConfig, + constraint_spec: ConstraintSpec | None = None, ): """Initialize the hybrid proposer. @@ -156,7 +157,12 @@ def __init__( rotor_model: Rotor model for the molecule. torsion_lib: Torsion library for biased sampling. config: Generation configuration. + constraint_spec: Optional positional constraints. When provided, + ETKDG seeding is replaced by seeding from an existing conformer, + global shake moves are suppressed, and minimization applies MMFF + position restraints to keep constrained atoms fixed. """ + self.constraint_spec = constraint_spec self.mol = mol self.rotor_model = rotor_model self.torsion_lib = torsion_lib @@ -194,6 +200,15 @@ def __init__( if self._staging_mmff_props is not None: self._staging_mmff_props.SetMMFFDielectricConstant(config.fast_dielectric) + # Store reference positions for constrained atoms so we can snap them + # back to exactly the starting pose after each minimization. The stiff + # position restraints keep drift tiny, but this eliminates it entirely. + self._constrained_ref_pos: dict[int, np.ndarray] = {} + if constraint_spec is not None and mol.GetNumConformers() > 0: + ref_conf = mol.GetConformer(mol.GetConformers()[0].GetId()) + for idx in constraint_spec.constrained_atoms: + self._constrained_ref_pos[idx] = np.array(ref_conf.GetAtomPosition(idx)) + # Set random seed if config.random_seed is not None: random.seed(config.random_seed) @@ -254,6 +269,105 @@ def generate_seeds(self, n_seeds: int) -> list[tuple[int, float]]: for cid in conf_ids ] + def _reset_constrained_positions(self, mol: Chem.Mol, conf_id: int) -> None: + """Snap constrained atoms back to their exact reference coordinates. + + Called after every constrained minimization to eliminate any residual + drift that the position restraints did not fully suppress. + + Args: + mol: RDKit molecule containing the conformer. + conf_id: Conformer ID to update in place. + """ + if not self._constrained_ref_pos: + return + conf = mol.GetConformer(conf_id) + for idx, pos in self._constrained_ref_pos.items(): + conf.SetAtomPosition(idx, pos.tolist()) + + def seed_from_conformer(self, conf_id: int) -> list[tuple[int, float]]: + """Use an existing conformer as the sole seed (constrained mode). + + Fast-minimizes the conformer with position restraints applied to + constrained atoms so the core stays at the MCS-aligned pose. + + Args: + conf_id: ID of the starting conformer already present in self.mol. + + Returns: + List of one (conf_id, energy_kcal) tuple, or empty list on failure. + """ + energy = self._minimize_constrained(self.mol, conf_id, use_fast=True) + if not np.isfinite(energy): + return [] + return [(conf_id, energy)] + + def _minimize_constrained(self, mol: Chem.Mol, conf_id: int, use_fast: bool = True) -> float: + """Minimize a conformer with MMFF position restraints on constrained atoms. + + Args: + mol: RDKit molecule containing the conformer. + conf_id: Conformer ID to minimize in place. + use_fast: If True, use fast_minimizer iterations; otherwise full. + + Returns: + Energy in kcal/mol after minimization, or inf on failure. + """ + assert self.constraint_spec is not None + minimizer = self.fast_minimizer if use_fast else self.full_minimizer + props = getattr(minimizer, "_mmff_props", None) + + try: + if props is not None: + ff = AllChem.MMFFGetMoleculeForceField(mol, props, confId=int(conf_id)) + if ff is None: + return float("inf") + for idx in self.constraint_spec.constrained_atoms: + ff.MMFFAddPositionConstraint(idx, 0.0, self.constraint_spec.position_force_constant) + ff.Minimize(maxIts=int(minimizer.max_iters)) + energy = float(ff.CalcEnergy()) + else: + # UFF fallback — no position constraints available, minimize freely + ff = AllChem.UFFGetMoleculeForceField(mol, confId=int(conf_id)) + if ff is None: + return float("inf") + ff.Minimize(maxIts=int(minimizer.max_iters)) + energy = float(ff.CalcEnergy()) + except (ValueError, RuntimeError): + return float("inf") + + # Snap constrained atoms back to exact reference coordinates, eliminating + # any residual drift that the position restraints did not fully suppress. + self._reset_constrained_positions(mol, conf_id) + return energy + + def _propose_constrained(self, pool: "ConformerPool", step: int) -> tuple[int, float, str] | None: + """Propose a single conformer using constrained minimization. + + Args: + pool: Conformer pool for parent selection. + step: Current step number. + + Returns: + Tuple of (conf_id, energy, source) or None if failed. + """ + result = self._generate_candidate(pool, step) + if result is None: + return None + new_conf_id, move_type = result + + try: + energy = self._minimize_constrained(self.mol, new_conf_id, use_fast=True) + except Exception: + self.mol.RemoveConformer(new_conf_id) + return None + + if not np.isfinite(energy): + self.mol.RemoveConformer(new_conf_id) + return None + + return (new_conf_id, energy, f"hybrid_{move_type}") + def _sample_angle(self, rotor_idx: int) -> float: """Sample an angle from the torsion library. @@ -400,12 +514,19 @@ def _select_move_type(self, step: int) -> str: Returns: Move type string. """ - # Periodic global shake - if step > 0 and step % self.config.shake_period == 0: + # Periodic global shake — suppressed in constrained mode (would thrash + # the carefully placed starting pose). + if self.constraint_spec is None and step > 0 and step % self.config.shake_period == 0: return "global_shake" - # Weighted random selection, suppressing ring_flip when no rings exist. + # Weighted random selection, suppressing unavailable move types. probs = dict(self.config.move_probs) + + # Remove global_shake from weighted pool in constrained mode. + if self.constraint_spec is not None: + extra = probs.pop("global_shake", 0.0) + probs["single_rotor"] = probs.get("single_rotor", 0.0) + extra + if not self.rotor_model.ring_flips and "ring_flip" in probs: extra = probs.pop("ring_flip") # redistribute to single_rotor @@ -497,6 +618,10 @@ def propose_batch(self, pool: ConformerPool, step: int) -> list[tuple[int, float staging props (fast_dielectric applied), then transfers accepted (finite-energy) conformers back to self.mol. + When constraint_spec is set, falls back to sequential per-conformer + minimization with MMFF position restraints (MMFFOptimizeMoleculeConfs + does not support custom force field terms). + Args: pool: Conformer pool for parent selection. step: Current step number (used for move-type selection of first item). @@ -504,6 +629,15 @@ def propose_batch(self, pool: ConformerPool, step: int) -> list[tuple[int, float Returns: List of (conf_id, energy, source) tuples for accepted conformers. """ + # Constrained mode: per-conformer MMFF with position restraints. + if self.constraint_spec is not None: + results: list[tuple[int, float, str]] = [] + for i in range(self.config.minimize_batch_size): + result = self._propose_constrained(pool, step + i) + if result is not None: + results.append(result) + return results + batch_size = self.config.minimize_batch_size # 1. Generate candidates on self.mol (move applied, clash-checked). @@ -599,6 +733,65 @@ def full_refine_final( return energies + def full_refine_final_constrained( + self, + mol: Chem.Mol, + final_ids: list[int], + max_iters: int = 200, + dielectric: float = 4.0, + ) -> list[float]: + """Run full MMFF minimization on final conformers with position restraints. + + Used in constrained mode so the core remains pinned to the MCS pose + during final refinement. + + Args: + mol: RDKit molecule containing the conformers to refine. + final_ids: Conformer IDs to refine. + max_iters: Maximum MMFF iterations. + dielectric: Dielectric constant for the refinement pass. + + Returns: + List of refined energies in kcal/mol, aligned to final_ids. + """ + assert self.constraint_spec is not None + + # Keep only finals before refining + final_set = set(final_ids) + for conf in list(mol.GetConformers()): + if conf.GetId() not in final_set: + mol.RemoveConformer(conf.GetId()) + + props = AllChem.MMFFGetMoleculeProperties(mol, mmffVariant="MMFF94s") + energies: list[float] = [] + + for cid in final_ids: + try: + if props is not None: + props.SetMMFFDielectricConstant(dielectric) + ff = AllChem.MMFFGetMoleculeForceField(mol, props, confId=int(cid)) + if ff is None: + energies.append(float("inf")) + continue + for idx in self.constraint_spec.constrained_atoms: + ff.MMFFAddPositionConstraint(idx, 0.0, self.constraint_spec.position_force_constant) + ff.Minimize(maxIts=int(max_iters)) + energies.append(float(ff.CalcEnergy())) + else: + ff = AllChem.UFFGetMoleculeForceField(mol, confId=int(cid)) + if ff is None: + energies.append(float("inf")) + continue + ff.Minimize(maxIts=int(max_iters)) + energies.append(float(ff.CalcEnergy())) + except (ValueError, RuntimeError): + energies.append(float("inf")) + continue + self._reset_constrained_positions(mol, cid) + + return energies + + def run_hybrid_generation( mol: Chem.Mol, rotor_model: RotorModel, @@ -614,16 +807,35 @@ def run_hybrid_generation( Returns: Tuple of (mol, conf_ids, energies). """ + constraint_spec = config.constraint_spec + + # Filter rotors before building the proposer so _rotor_angles is computed + # only for free rotors. + if constraint_spec is not None: + rotor_model = filter_constrained_rotors(rotor_model, constraint_spec.constrained_atoms) + torsion_lib = TorsionLibrary() - proposer = HybridProposer(mol, rotor_model, torsion_lib, config) + proposer = HybridProposer(mol, rotor_model, torsion_lib, config, constraint_spec=constraint_spec) pool = ConformerPool(mol, config) - # Resolve seed count: explicit config value or auto-computed from topology. - n_seeds = config.n_seeds if config.n_seeds is not None else _compute_n_seeds(rotor_model, config.seed_n_per_rotor) - seeds = proposer.generate_seeds(n_seeds) + if constraint_spec is not None: + # Constrained mode: seed from the single starting conformer already in mol. + existing_ids = [c.GetId() for c in mol.GetConformers()] + if not existing_ids: + raise ValueError( + "Constrained conformer generation requires a starting conformer. " + "Use generate_conformers_from_pose to supply one." + ) + seeds = proposer.seed_from_conformer(existing_ids[0]) + seed_source = "seed_pose" + else: + # Standard mode: ETKDG seeding. + n_seeds = config.n_seeds if config.n_seeds is not None else _compute_n_seeds(rotor_model, config.seed_n_per_rotor) + seeds = proposer.generate_seeds(n_seeds) + seed_source = "seed_etkdg" for conf_id, energy in seeds: - pool.insert(conf_id, energy, source="seed_etkdg") + pool.insert(conf_id, energy, source=seed_source) # Run exploration. # Batch mode: accumulate minimize_batch_size proposals, minimize in parallel. @@ -653,9 +865,14 @@ def run_hybrid_generation( # Full refinement on the final set (optional — skip for docking-prep workflows). if config.do_final_refine: - final_energies = proposer.full_refine_final( - mol, final_ids, config.num_threads, config.max_minimization_iters, dielectric=config.final_dielectric - ) + if constraint_spec is not None: + final_energies = proposer.full_refine_final_constrained( + mol, final_ids, config.max_minimization_iters, dielectric=config.final_dielectric + ) + else: + final_energies = proposer.full_refine_final( + mol, final_ids, config.num_threads, config.max_minimization_iters, dielectric=config.final_dielectric + ) else: # Return the fast-minimized energies already stored in the pool. energy_map = {cid: (rec.energy_kcal or float("inf")) for cid, rec in pool.records.items()} diff --git a/pyproject.toml b/pyproject.toml index 923fd72..9a5e562 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,7 +1,7 @@ [project] name = "openconf" description = "Conformer generation for drug-like molecules" -version = "0.0.1" +version = "0.0.2" readme = "README.md" license = "MIT" authors = [{name = "Corin Wagen"}] diff --git a/tests/test_constrained.py b/tests/test_constrained.py new file mode 100644 index 0000000..9045374 --- /dev/null +++ b/tests/test_constrained.py @@ -0,0 +1,307 @@ +"""Tests for pose-constrained (FEP-style analogue) conformer generation.""" + +import numpy as np +import pytest +from rdkit import Chem +from rdkit.Chem import AllChem + + +def _make_butylbenzene_pose() -> tuple[Chem.Mol, list[int]]: + """Return an H-added butylbenzene with one embedded conformer and ring atom indices.""" + mol = Chem.MolFromSmiles("CCCCc1ccccc1") + mol = Chem.AddHs(mol) + AllChem.EmbedMolecule(mol, AllChem.ETKDGv3()) + ring_atoms = [4, 5, 6, 7, 8, 9] # benzene ring heavy-atom indices + return mol, ring_atoms + + +# --------------------------------------------------------------------------- +# Imports / preset +# --------------------------------------------------------------------------- + + +def test_imports(): + """New public symbols are importable from the top-level package.""" + from openconf import ConstraintSpec, filter_constrained_rotors, generate_conformers_from_pose # noqa: F401 + + +def test_analogue_preset_values(): + """'analogue' preset returns the expected config.""" + from openconf import preset_config + + cfg = preset_config("analogue") + assert cfg.max_out == 50 + assert cfg.n_steps == 150 + assert cfg.do_final_refine is True + assert cfg.parent_strategy == "softmax" + assert cfg.final_select == "diverse" + + +# --------------------------------------------------------------------------- +# Core correctness +# --------------------------------------------------------------------------- + + +def test_constrained_atoms_exactly_fixed(): + """Constrained atom positions must be bit-for-bit identical to the input pose.""" + from openconf import generate_conformers_from_pose + + mol, ring_atoms = _make_butylbenzene_pose() + + ref_pos = {i: np.array(mol.GetConformer(0).GetAtomPosition(i)) for i in ring_atoms} + ensemble = generate_conformers_from_pose(mol, constrained_atoms=ring_atoms) + + for rec in ensemble.records: + conf = ensemble.mol.GetConformer(rec.conf_id) + for i in ring_atoms: + pos = np.array(conf.GetAtomPosition(i)) + assert np.allclose(pos, ref_pos[i], atol=1e-9), ( + f"Atom {i} drifted in conformer {rec.conf_id}: " + f"ref={ref_pos[i]}, got={pos}" + ) + + +def test_generate_conformers_from_pose_returns_conformers(): + """Basic smoke test: function runs and produces at least one conformer.""" + from openconf import generate_conformers_from_pose + + mol, ring_atoms = _make_butylbenzene_pose() + ensemble = generate_conformers_from_pose(mol, constrained_atoms=ring_atoms) + + assert ensemble.n_conformers > 0 + assert ensemble.n_conformers <= 50 # analogue preset max_out + assert all(np.isfinite(e) for e in ensemble.energies) + + +def test_free_atoms_do_move(): + """At least one free (non-constrained) atom should differ across conformers.""" + from openconf import generate_conformers_from_pose + + mol, ring_atoms = _make_butylbenzene_pose() + ensemble = generate_conformers_from_pose(mol, constrained_atoms=ring_atoms) + + # Butyl chain atom indices (heavy atoms 0-3) + free_atoms = [0, 1, 2, 3] + + # Collect positions of free atoms across all conformers + all_positions = [] + for rec in ensemble.records: + conf = ensemble.mol.GetConformer(rec.conf_id) + pos = np.array([conf.GetAtomPosition(i) for i in free_atoms]) + all_positions.append(pos) + + # At least two conformers should differ in at least one free atom + assert len(all_positions) >= 2 + assert not all( + np.allclose(all_positions[0], p, atol=0.01) for p in all_positions[1:] + ), "All conformers have identical free-atom coordinates — no exploration occurred" + + +# --------------------------------------------------------------------------- +# Input variants +# --------------------------------------------------------------------------- + + +def test_mol_without_hs_input(): + """Works when the input mol has only heavy atoms (AddHs applied internally).""" + from openconf import generate_conformers_from_pose + + mol = Chem.MolFromSmiles("CCCCc1ccccc1") + # Embed without adding Hs — only heavy-atom conformer + AllChem.EmbedMolecule(mol, AllChem.ETKDGv3()) + ring_atoms = [4, 5, 6, 7, 8, 9] + + ensemble = generate_conformers_from_pose(mol, constrained_atoms=ring_atoms) + assert ensemble.n_conformers > 0 + + +def test_first_conformer_used_as_seed(): + """When mol has multiple conformers, only the first is used as the seed.""" + from openconf import ConformerConfig, generate_conformers_from_pose + + mol = Chem.MolFromSmiles("CCCCc1ccccc1") + mol = Chem.AddHs(mol) + AllChem.EmbedMultipleConfs(mol, numConfs=3, randomSeed=0) + ring_atoms = [4, 5, 6, 7, 8, 9] + + ref_pos = {i: np.array(mol.GetConformer(0).GetAtomPosition(i)) for i in ring_atoms} + config = ConformerConfig(max_out=5, n_steps=20, pool_max=50) + ensemble = generate_conformers_from_pose(mol, constrained_atoms=ring_atoms, config=config) + + for rec in ensemble.records: + conf = ensemble.mol.GetConformer(rec.conf_id) + for i in ring_atoms: + pos = np.array(conf.GetAtomPosition(i)) + assert np.allclose(pos, ref_pos[i], atol=1e-9) + + +def test_no_conformer_raises(): + """Passing a mol with no conformer raises ValueError.""" + from openconf import generate_conformers_from_pose + + mol = Chem.MolFromSmiles("CCCCc1ccccc1") + mol = Chem.AddHs(mol) + # deliberately do NOT embed + + with pytest.raises(ValueError, match="conformer"): + generate_conformers_from_pose(mol, constrained_atoms=[4, 5, 6, 7, 8, 9]) + + +def test_config_and_preset_both_raises(): + """Passing both config and preset raises ValueError.""" + from openconf import ConformerConfig, generate_conformers_from_pose + + mol, ring_atoms = _make_butylbenzene_pose() + with pytest.raises(ValueError): + generate_conformers_from_pose( + mol, + constrained_atoms=ring_atoms, + config=ConformerConfig(), + preset="analogue", + ) + + +# --------------------------------------------------------------------------- +# Rotor filtering +# --------------------------------------------------------------------------- + + +def test_filter_constrained_rotors_eliminates_double_sided(): + """A rotor whose both fragments contain constrained atoms is excluded. + + In butane (CCCC) with heavy atoms {1, 2} constrained, the central C1–C2 + bond has constrained atoms on both sides and must be eliminated. + The terminal bonds C0–C1 and C2–C3 each have one free side and are kept. + """ + from openconf import build_rotor_model, filter_constrained_rotors, prepare_molecule + + mol = prepare_molecule(Chem.MolFromSmiles("CCCC")) + rm = build_rotor_model(mol) + full_count = rm.n_rotatable # 3 bonds: C0-C1, C1-C2, C2-C3 + + # Constrain the two middle carbons; the central bond has constrained atoms + # on both sides and should be removed. + constrained = frozenset([1, 2]) + filtered = filter_constrained_rotors(rm, constrained) + + assert filtered.n_rotatable < full_count + assert filtered.n_rotatable > 0 # terminal bonds should remain + + +def test_filter_constrained_rotors_free_side_reoriented(): + """When the free fragment is on the 'wrong' side of a rotor it is flipped. + + Butylbenzene with only ring atoms constrained: the ring-to-chain bond's + free side (butyl chain) becomes the moving side. All 4 rotors are kept. + """ + from openconf import build_rotor_model, filter_constrained_rotors, prepare_molecule + from openconf.perceive import _atoms_on_side + + mol = prepare_molecule(Chem.MolFromSmiles("CCCCc1ccccc1")) + rm = build_rotor_model(mol) + + constrained = frozenset(range(4, 10)) # ring atoms + filtered = filter_constrained_rotors(rm, constrained) + + # All 4 chain rotors should be preserved (just possibly reoriented). + assert filtered.n_rotatable == rm.n_rotatable + + # For every remaining rotor the moving fragment must be entirely free. + for rotor in filtered.rotors: + atom_i, atom_j = rotor.atom_idxs + moving = _atoms_on_side(mol, atom_j, atom_i) + assert not constrained & moving, ( + f"Rotor {rotor.atom_idxs} has constrained atoms in moving fragment: " + f"{constrained & moving}" + ) + + +def test_filter_constrained_rotors_no_constrained_atoms(): + """Empty constraint set leaves the rotor model unchanged.""" + from openconf import build_rotor_model, filter_constrained_rotors, prepare_molecule + + mol = prepare_molecule(Chem.MolFromSmiles("CCCCc1ccccc1")) + rm = build_rotor_model(mol) + filtered = filter_constrained_rotors(rm, frozenset()) + + assert filtered.n_rotatable == rm.n_rotatable + + +def test_filter_constrained_ring_flips(): + """Ring flips involving constrained ring atoms are removed.""" + from openconf import build_rotor_model, filter_constrained_rotors, prepare_molecule + + # Cyclohexyl-methyl: the ring should flip; constrain the ring + mol = prepare_molecule(Chem.MolFromSmiles("CC1CCCCC1")) + rm = build_rotor_model(mol) + assert len(rm.ring_flips) == 1 + + # Constrain all ring atoms (indices vary with H-addition; use atom symbols) + ring_indices = frozenset( + i for i, a in enumerate(mol.GetAtoms()) + if a.IsInRing() and a.GetAtomicNum() == 6 + ) + filtered = filter_constrained_rotors(rm, ring_indices) + assert len(filtered.ring_flips) == 0, "Ring flip should be removed when ring is constrained" + + +def test_filter_free_ring_flips_preserved(): + """Ring flips whose atoms are all free are kept after filtering.""" + from openconf import build_rotor_model, filter_constrained_rotors, prepare_molecule + + # Butyl-cyclohexane: constrain only the butyl chain, ring stays free + mol = prepare_molecule(Chem.MolFromSmiles("CCCCC1CCCCC1")) + rm = build_rotor_model(mol) + assert len(rm.ring_flips) == 1 + + # Constrain first 4 heavy atoms (butyl chain, not the ring) + constrained = frozenset([0, 1, 2, 3]) + filtered = filter_constrained_rotors(rm, constrained) + assert len(filtered.ring_flips) == 1, "Free ring flip should be preserved" + + +# --------------------------------------------------------------------------- +# Move suppression +# --------------------------------------------------------------------------- + + +def test_global_shake_suppressed_in_constrained_mode(): + """_select_move_type never returns 'global_shake' when constraint_spec is set.""" + from openconf.config import ConformerConfig, ConstraintSpec + from openconf.perceive import build_rotor_model, prepare_molecule + from openconf.propose.hybrid import HybridProposer + from openconf.torsionlib import TorsionLibrary + + mol = prepare_molecule(Chem.MolFromSmiles("CCCCc1ccccc1")) + AllChem.EmbedMolecule(mol, randomSeed=0) + rm = build_rotor_model(mol) + + spec = ConstraintSpec(constrained_atoms=frozenset(range(4, 10))) + config = ConformerConfig(random_seed=0, shake_period=1) # shake_period=1 forces periodic shakes + proposer = HybridProposer(mol, rm, TorsionLibrary(), config, constraint_spec=spec) + + # Run enough steps to trigger shake_period multiple times + move_types = {proposer._select_move_type(step) for step in range(200)} + assert "global_shake" not in move_types, f"global_shake appeared in constrained mode: {move_types}" + + +# --------------------------------------------------------------------------- +# SDF output +# --------------------------------------------------------------------------- + + +def test_constrained_ensemble_to_sdf(tmp_path): + """Constrained ensemble can be written to and read back from SDF.""" + from openconf import ConformerConfig, generate_conformers_from_pose + + mol, ring_atoms = _make_butylbenzene_pose() + config = ConformerConfig(max_out=5, n_steps=20, pool_max=50) + ensemble = generate_conformers_from_pose(mol, constrained_atoms=ring_atoms, config=config) + + out = tmp_path / "analogue.sdf" + ensemble.to_sdf(str(out)) + assert out.exists() + + supplier = Chem.SDMolSupplier(str(out), removeHs=False) + mols_read = [m for m in supplier if m is not None] + assert len(mols_read) == ensemble.n_conformers From 1b65ab4615e8f57d037e5d377c0e2628a9f6d201 Mon Sep 17 00:00:00 2001 From: Corin Wagen Date: Fri, 3 Apr 2026 15:04:44 -0400 Subject: [PATCH 2/2] lints etc --- openconf/api.py | 8 ++------ openconf/perceive.py | 4 +--- openconf/propose/hybrid.py | 5 +++-- pyproject.toml | 2 ++ tests/test_constrained.py | 17 ++++++----------- uv.lock | 2 +- 6 files changed, 15 insertions(+), 23 deletions(-) diff --git a/openconf/api.py b/openconf/api.py index 55f831f..e3e662d 100644 --- a/openconf/api.py +++ b/openconf/api.py @@ -246,18 +246,14 @@ def generate_conformers_from_pose( # Resolve config — default to "analogue" preset for this entry point. if config is not None: - resolved_config = ConformerConfig(**{ - k: v for k, v in config.__dict__.items() if k != "constraint_spec" - }) + resolved_config = ConformerConfig(**{k: v for k, v in config.__dict__.items() if k != "constraint_spec"}) elif preset is not None: resolved_config = preset_config(preset) else: resolved_config = preset_config("analogue") # Attach the constraint spec (overrides any constraint_spec already in config). - resolved_config.constraint_spec = ConstraintSpec( - constrained_atoms=frozenset(constrained_atoms) - ) + resolved_config.constraint_spec = ConstraintSpec(constrained_atoms=frozenset(constrained_atoms)) # Prepare molecule — AddHs is a no-op if Hs are already present. prepped_mol = prepare_molecule(mol, add_hs=True) diff --git a/openconf/perceive.py b/openconf/perceive.py index 12ae5c6..37559e1 100644 --- a/openconf/perceive.py +++ b/openconf/perceive.py @@ -396,9 +396,7 @@ def filter_constrained_rotors(rotor_model: "RotorModel", constrained_atoms: froz # else: constrained atoms on both sides — exclude this rotor # Ring flips: keep only rings with no constrained atoms - free_ring_flips = [ - rf for rf in rotor_model.ring_flips if not constrained_atoms & frozenset(rf.ring_atoms) - ] + free_ring_flips = [rf for rf in rotor_model.ring_flips if not constrained_atoms & frozenset(rf.ring_atoms)] adj = _build_rotor_adjacency(free_rotors, mol) diff --git a/openconf/propose/hybrid.py b/openconf/propose/hybrid.py index 77d7eeb..aff9cec 100644 --- a/openconf/propose/hybrid.py +++ b/openconf/propose/hybrid.py @@ -732,7 +732,6 @@ def full_refine_final( energies.append(ff.CalcEnergy() if ff else float("inf")) return energies - def full_refine_final_constrained( self, mol: Chem.Mol, @@ -830,7 +829,9 @@ def run_hybrid_generation( seed_source = "seed_pose" else: # Standard mode: ETKDG seeding. - n_seeds = config.n_seeds if config.n_seeds is not None else _compute_n_seeds(rotor_model, config.seed_n_per_rotor) + n_seeds = ( + config.n_seeds if config.n_seeds is not None else _compute_n_seeds(rotor_model, config.seed_n_per_rotor) + ) seeds = proposer.generate_seeds(n_seeds) seed_source = "seed_etkdg" diff --git a/pyproject.toml b/pyproject.toml index 9a5e562..7ac6507 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -65,6 +65,8 @@ ignore = [ "PLR0915", # Too many statements "PLR1702", # Too many nested-blocks "PLR2004", # Magic value in comparison + "PT011", # `pytest.raises(ValueError)` is too broad + "RUF002", # Docstring contains ambiguous `–` (EN DASH). "TID252", # Prefer absolute over relative imports ] diff --git a/tests/test_constrained.py b/tests/test_constrained.py index 9045374..f4e7b1f 100644 --- a/tests/test_constrained.py +++ b/tests/test_constrained.py @@ -56,8 +56,7 @@ def test_constrained_atoms_exactly_fixed(): for i in ring_atoms: pos = np.array(conf.GetAtomPosition(i)) assert np.allclose(pos, ref_pos[i], atol=1e-9), ( - f"Atom {i} drifted in conformer {rec.conf_id}: " - f"ref={ref_pos[i]}, got={pos}" + f"Atom {i} drifted in conformer {rec.conf_id}: ref={ref_pos[i]}, got={pos}" ) @@ -92,9 +91,9 @@ def test_free_atoms_do_move(): # At least two conformers should differ in at least one free atom assert len(all_positions) >= 2 - assert not all( - np.allclose(all_positions[0], p, atol=0.01) for p in all_positions[1:] - ), "All conformers have identical free-atom coordinates — no exploration occurred" + assert not all(np.allclose(all_positions[0], p, atol=0.01) for p in all_positions[1:]), ( + "All conformers have identical free-atom coordinates — no exploration occurred" + ) # --------------------------------------------------------------------------- @@ -211,8 +210,7 @@ def test_filter_constrained_rotors_free_side_reoriented(): atom_i, atom_j = rotor.atom_idxs moving = _atoms_on_side(mol, atom_j, atom_i) assert not constrained & moving, ( - f"Rotor {rotor.atom_idxs} has constrained atoms in moving fragment: " - f"{constrained & moving}" + f"Rotor {rotor.atom_idxs} has constrained atoms in moving fragment: {constrained & moving}" ) @@ -237,10 +235,7 @@ def test_filter_constrained_ring_flips(): assert len(rm.ring_flips) == 1 # Constrain all ring atoms (indices vary with H-addition; use atom symbols) - ring_indices = frozenset( - i for i, a in enumerate(mol.GetAtoms()) - if a.IsInRing() and a.GetAtomicNum() == 6 - ) + ring_indices = frozenset(i for i, a in enumerate(mol.GetAtoms()) if a.IsInRing() and a.GetAtomicNum() == 6) filtered = filter_constrained_rotors(rm, ring_indices) assert len(filtered.ring_flips) == 0, "Ring flip should be removed when ring is constrained" diff --git a/uv.lock b/uv.lock index aa2023b..df19d5e 100644 --- a/uv.lock +++ b/uv.lock @@ -269,7 +269,7 @@ wheels = [ [[package]] name = "openconf" -version = "0.0.1" +version = "0.0.2" source = { editable = "." } dependencies = [ { name = "numpy" },