From 4da95ff5d290d4eb7643d594b041547d26185637 Mon Sep 17 00:00:00 2001 From: ncasetti7 Date: Mon, 15 Jun 2026 14:09:49 -0400 Subject: [PATCH 1/4] Add macrocycle preset, low mode following, macrocyclic amide flip --- README.md | 46 +-- SCIENCE.md | 43 +-- openconf/config.py | 80 +++++- openconf/data/runtime_tuning.json | 12 +- openconf/dedupe.py | 5 +- openconf/pool.py | 1 + openconf/propose/__init__.py | 3 +- openconf/propose/amide_seeds.py | 66 +++++ openconf/propose/hybrid.py | 37 +++ openconf/propose/low_mode.py | 296 +++++++++++++++++++ openconf/propose/moves.py | 64 +++++ openconf/propose/stats.py | 2 + openconf/tuning.py | 2 + tests/test_low_mode.py | 457 ++++++++++++++++++++++++++++++ tests/test_macrocycles.py | 113 ++++++++ uv.lock | 2 +- 16 files changed, 1182 insertions(+), 47 deletions(-) create mode 100644 openconf/propose/low_mode.py create mode 100644 tests/test_low_mode.py diff --git a/README.md b/README.md index 94e79da..ae40ceb 100644 --- a/README.md +++ b/README.md @@ -39,7 +39,7 @@ ensemble.to_xyz("output.xyz") ### Named Presets -Five use-case presets are available out of the box: +Six use-case presets are available out of the box: ```python from openconf import generate_conformers @@ -48,10 +48,13 @@ ensemble = generate_conformers(mol, preset="rapid") # fast virtual scree ensemble = generate_conformers(mol, preset="ensemble") # property prediction ensemble = generate_conformers(mol, preset="spectroscopic") # NMR / IR / VCD ensemble = generate_conformers(mol, preset="docking") # docking pose recovery +ensemble = generate_conformers(mol, preset="macrocycle") # macrocyclic ring systems ``` For FEP-style analogue generation from a fixed pose, see [`generate_conformers_from_pose`](#5-analogue--fep-style-r-group-exploration) below. +For macrocyclic ring systems, use the `"macrocycle"` preset, which enables low-mode following and a wide energy window to capture the full range of ring-pucker conformations. + ### Custom Configuration For full control, pass a `ConformerConfig` directly. You can also use @@ -94,7 +97,7 @@ config.max_out = 200 # override a single field ensemble = generate_conformers(mol, config=config) ``` -Available presets: `"rapid"`, `"ensemble"`, `"spectroscopic"`, `"docking"`, `"analogue"`. +Available presets: `"rapid"`, `"ensemble"`, `"spectroscopic"`, `"docking"`, `"analogue"`, `"macrocycle"`. Below are representative wall-clock timings measured on a single CPU core (Apple M2 Pro), mean over 3 runs. @@ -279,7 +282,7 @@ ensemble.to_sdf("output.sdf")
Full config equivalent ```python -from openconf import ConformerConfig, PrismConfig, generate_conformers +from openconf import ConformerConfig, generate_conformers config = ConformerConfig( max_out=250, @@ -292,7 +295,6 @@ config = ConformerConfig( minimize_batch_size=8, parent_strategy="uniform", final_select="diverse", - prism_config=PrismConfig(energy_window_kcal=18.0), ) ensemble = generate_conformers("CC(C)Cc1ccc(cc1)C(C)C(=O)O", config=config) ensemble.to_sdf("docking_input.sdf") @@ -354,7 +356,7 @@ ensemble = generate_conformers_from_pose(mol, constrained_atoms=ring_atoms, conf
Full analogue preset equivalent ```python -from openconf import ConformerConfig, PrismConfig, generate_conformers_from_pose +from openconf import ConformerConfig, generate_conformers_from_pose config = ConformerConfig( max_out=50, @@ -365,7 +367,6 @@ config = ConformerConfig( 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) ``` @@ -376,19 +377,22 @@ ensemble = generate_conformers_from_pose(mol, constrained_atoms=ring_atoms, conf ### Configuration Quick Reference -| 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 | +| Parameter | Rapid | Ensemble | Spectroscopic | Docking | Analogue | Macrocycle | +|---|---|---|---|---|---|---| +| `max_out` | 5 | 50 | 100 | 250 | 50 | 200 | +| `n_steps` | 30 | 200 | 400 | 500 | 150 | 500 | +| `energy_window_kcal` | 20 | 10 | 5 | 18 | 10 | 100 | +| `seed_n_per_rotor` | 2 | 3 | 5 | 4 | — | 3 | +| `seed_prune_rms_thresh` | 1.5 | 1.0 | 0.5 | 0.8 | — | 1.0 | +| `do_final_refine` | False | True | True | False | True | True | +| `parent_strategy` | softmax | softmax | softmax | uniform | softmax | softmax | +| `final_select` | diverse | diverse | energy | diverse | diverse | diverse | +| `use_low_mode_following` | False | False | False | False | False | True | *Analogue mode uses `generate_conformers_from_pose`; seeding parameters are unused because ETKDG is skipped.* +`use_low_mode_following` is enabled automatically by the `"macrocycle"` preset. Enable it manually via `ConformerConfig(use_low_mode_following=True)` for other highly coupled systems. See [SCIENCE.md](SCIENCE.md) for the full parameter reference. + ## How It Works ### 1. Seeding @@ -397,6 +401,8 @@ Generates initial conformers using RDKit's ETKDGv3 algorithm. The seed count is When `n_seeds=None`, openconf resolves a seed plan before embedding. Explicit `n_seeds` values are always honored. Macrocycles keep dense topology-derived seed budgets for ring-pucker discovery, while simple low-flexibility acyclic molecules and large flexible hydrocarbons use data-backed reduced budgets to avoid redundant RDKit embeddings. +When `use_low_mode_following=True`, an additional seeding step runs after ETKDG: openconf numerically evaluates the Hessian at each source seed, identifies soft eigenvectors (eigenvalue < `low_mode_eigenvalue_threshold`), and scans displaced geometries along each mode in both directions. Each scan endpoint is minimized to a new local minimum, providing seeds that capture collective ring-puckering and correlated torsion motions that independent torsion moves miss. This mirrors the LMOD procedure (Kolossváry & Guida, *JACS*, 1996) and is most effective for macrocycles and other highly coupled systems. Disabled by default; Hessian evaluation costs 6N MMFF gradient calls per seed. The `"macrocycle"` preset enables low-mode following automatically. See [SCIENCE.md](SCIENCE.md) for full parameter details. + ### 2. Hybrid Exploration The default "hybrid" strategy combines: @@ -405,6 +411,8 @@ The default "hybrid" strategy combines: - MCMM moves: random torsion perturbations biased by the library - Correlated moves: simultaneous changes to adjacent rotors - Ring flip moves: SVD plane-reflection of non-aromatic 5–7-membered rings to sample chair/envelope inversions +- Macrocycle ring moves: closure-preserving crankshaft and kinematic ring-closure (KIC) moves that rotate ring arcs without breaking bond lengths +- Amide-flip moves: 180° rotation about an in-ring amide C–N bond to invert its cis/trans configuration in macrocycles - Global shakes: periodic large perturbations to escape local basins ### 3. Minimization @@ -425,20 +433,18 @@ For the full algorithm description and parameter tuning guide, see [SCIENCE.md]( Validated on the [Iridium benchmark](docs/benchmark_report.md) (120 drug-like molecules, bioactive conformer recovery from crystal structures). At N=200, openconf achieves a median best-RMSD of 0.58 Å vs. 0.63 Å for ETKDG+MMFF94s, at 10–15× lower wall time. The advantage is concentrated in flexible molecules (7–9 rotatable bonds), where openconf's torsion-library biasing and ring flip moves outperform pure distance-geometry seeding. -openconf is not recommended for macrocycles (ring size ≥ 12). On macrocyclic ring systems, ETKDG+MMFF94s with a large conformer budget outperforms openconf in both RMSD and ensemble coverage metrics. ETKDGv3 has dedicated macrocycle distance-geometry bounds that openconf does not replicate. +For macrocyclic systems (ring size ≥ 12) the default configuration trails ETKDG+MMFF94s in both RMSD and ensemble coverage metrics; ETKDGv3 has dedicated macrocycle distance-geometry bounds that the default MCMM exploration does not replicate. Setting `use_low_mode_following=True` in `ConformerConfig` adds Hessian-guided seeds targeting collective ring-puckering modes and is the recommended path for challenging macrocyclic systems where the default configuration is insufficient. ## API Reference ### Main Functions -- `generate_conformers(mol, method="hybrid", config=None)` - Main entry point -- `generate_conformers_from_smiles(smiles, ...)` - Convenience wrapper +- `generate_conformers(mol, method="hybrid", config=None)` - Main entry point; accepts a `Chem.Mol` or SMILES string - `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 1a4b249..e338789 100644 --- a/SCIENCE.md +++ b/SCIENCE.md @@ -10,10 +10,34 @@ Ring conformations require special treatment beyond simple torsion moves. For no For non-aromatic rings of size ≥ 6, openconf also applies a *crankshaft move*: two non-adjacent ring atoms are chosen as anchors, and the arc of ring atoms between them is rotated rigidly about the axis passing through the anchors. Because both anchors lie on the rotation axis, every bond length in the ring is preserved exactly — only the bond angles at the anchors deform, and those relax in the subsequent MMFF minimization. This gives a closure-preserving macrocycle move that single-torsion rotations cannot match. Rotation angles are drawn from [30°, 120°] with random sign; the arc length is drawn uniformly from 1 to n−3 ring atoms. Substituent subtrees rotate rigidly with their host ring atom. The crankshaft is selected with 12% probability per step when a crankable ring exists, and the move type is exempt from the usual static clash filter — post-rotation geometries are often momentarily strained but MMFF relaxes them reliably. Seed pruning (`pruneRmsThresh`) is also disabled automatically when a macrocycle is present, because the 1.0 Å default collapses distinct puckers before they reach MMFF. +Macrocyclic backbones often hinge on the cis/trans configuration of their amide bonds, and the cis↔trans barrier is far too high for ordinary torsion jitter to cross during exploration. openconf therefore adds an *amide-flip move* targeting in-ring amide bonds (a non-aromatic carbonyl carbon single-bonded to nitrogen, both in a ring of ≥ 8 atoms; both secondary and tertiary amides qualify). The move rotates the rest of the ring 180° (± the usual torsion jitter) about the carbonyl-carbon→nitrogen axis, inverting the amide while — exactly as in the crankshaft — preserving every ring bond length, since both the carbonyl carbon and nitrogen lie on the rotation axis. Only the bond angles at those two atoms strain, and MMFF relaxes them. Each flippable amide bond is associated with the largest ring containing it (the macrocycle in fused systems) so the most flexible arc moves. Like the other ring moves it is selected with a small per-step probability when an in-ring amide exists and is exempt from the static clash filter. This complements the cis/trans seed enumeration: the seeds plant both configurations up front, while the move lets the walk discover a flip after other torsional rearrangements have accumulated. + Seed count is computed automatically from the molecular topology rather than set to a fixed value. The base formula is `max(20, n_rotatable × 3)` (controlled by `seed_n_per_rotor`, default 3), plus 5 seeds per flippable ring and `ring_size²` seeds per macrocycle ring, capped at 500. A simple drug-like molecule with 8 rotatable bonds gets ~24 seeds; a steroid with three non-aromatic rings gets ~35; a 12-membered macrocycle gets ~164; a 16-membered macrocycle gets ~276. The super-linear macrocycle term reflects the fact that the low-energy pucker fraction drops rapidly with ring size, so dense seeding is the cheapest way to ensure the global basin is sampled. For simple low-flexibility acyclic molecules and large flexible hydrocarbons, openconf applies data-backed seed-plan reductions so redundant RDKit embeddings are skipped and MC torsion moves do more of the exploration. You can override this by setting `n_seeds` explicitly in `ConformerConfig`. 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. +## Low-Mode Following (optional) + +When `use_low_mode_following=True`, openconf supplements the ETKDG seeds with a Hessian-guided step inspired by the LMOD procedure of Kolossváry & Guida (*JACS*, 1996). After each source seed is minimized, openconf numerically evaluates the 3N×3N Hessian matrix of the MMFF energy via central finite differences of the gradient: each Cartesian coordinate is displaced by ±`fd_step` (default 0.005 Å) and the gradient at the two displaced geometries is used to estimate one row of the Hessian. The resulting matrix is symmetrized. + +Diagonalizing the Hessian yields a spectrum of normal modes. The first several eigenvalues are near zero (rigid-body translations and rotations); openconf detects the rigid-body/conformational boundary via the largest eigenvalue gap in the first eight entries rather than assuming exactly six zero modes, which correctly handles linear molecules and cases where minimization is not fully converged. The next eigenvectors — the **soft modes** — correspond to collective conformational motions: ring puckerings, coupled backbone torsions, and other deformations that cannot be decomposed into independent single-bond rotations. + +For each soft mode whose eigenvalue is below `low_mode_eigenvalue_threshold` (default 100 kcal/mol/Ų, capturing torsional and ring-puckering modes while excluding bond stretches), openconf scans the geometry in both the positive and negative directions along the eigenvector. At each step the geometry is displaced by `low_mode_scan_step_size` Å (3N Euclidean norm), the MMFF energy is evaluated, and scanning continues until the per-step energy increase exceeds `low_mode_scan_energy_threshold` (default ~2390 kcal/mol, equivalent to the paper's 10 000 kJ/mol) or `low_mode_scan_max_steps` is reached. The large default energy threshold means scanning passes through conformational barriers and stops only at severe steric clashes, placing the displaced endpoint on the far side of a barrier where a new local minimum can be found. Each scan endpoint that differs from the starting geometry is minimized; endpoints that minimize back to the starting geometry are discarded. + +The dominant cost is Hessian evaluation: 6N MMFF force-field constructions for an N-atom molecule. `low_mode_n_source_seeds` (default 1) caps how many seeds receive this treatment — only the lowest-energy seeds are selected. With the default settings of `low_mode_n_source_seeds=1` and `low_mode_max_modes=5`, at most 10 new seeds (2 directions × 5 modes) are added per Hessian evaluation. + +**When to enable:** macrocycles and other highly coupled flexible systems where independent single-rotor moves miss collective soft motions. The `"macrocycle"` preset enables low-mode following automatically alongside a wide 100 kcal/mol energy window. For simple acyclic drug-like molecules the MCMM exploration already covers these degrees of freedom and the Hessian cost is not worth paying. + +| Parameter | Default | Meaning | +|---|---|---| +| `use_low_mode_following` | `False` | Enable low-mode seeding | +| `low_mode_eigenvalue_threshold` | 100.0 | Eigenvalue cutoff (kcal/mol/Ų) for soft-mode selection | +| `low_mode_max_modes` | 5 | Maximum soft modes to scan per source seed | +| `low_mode_scan_step_size` | 0.25 | Displacement per scan step in Å (3N Euclidean norm) | +| `low_mode_scan_energy_threshold` | 2390.0 | Per-step ΔE limit before stopping scan (kcal/mol) | +| `low_mode_scan_max_steps` | 10 | Maximum scan steps per direction | +| `low_mode_n_source_seeds` | 1 | Number of minimized seeds from which Hessians are computed | + ## 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`. @@ -57,6 +81,7 @@ ensemble = generate_conformers(mol, preset="rapid") # fast virtual scree ensemble = generate_conformers(mol, preset="ensemble") # property prediction (50 conformers) ensemble = generate_conformers(mol, preset="spectroscopic") # Boltzmann ensemble for NMR/IR ensemble = generate_conformers(mol, preset="docking") # maximize bioactive recall +ensemble = generate_conformers(mol, preset="macrocycle") # macrocyclic ring systems ``` Each preset is a fully specified `ConformerConfig`; you can inspect and override individual @@ -112,23 +137,7 @@ config.energy_window_kcal = 20.0 ensemble = generate_conformers(mol, config=config) ``` -Macrocycles (ring size ≥ 10): openconf uses ETKDGv3 macrocycle distance bounds at seed time, disables seed RMSD pruning (which otherwise collapses distinct puckers), and applies crankshaft moves during exploration. On cycloalkanes C6–C14 the default config recovers the ETKDG+MMFF94s reference minimum within 1 kcal/mol, and on C16 within ~0.2 kcal/mol. Cyclododecane is an edge case — its global [3333] basin occupies <1% of seeds and needs `n_seeds ≥ 300` for reliable recovery (default is ~164); override `n_seeds` if you need the exact global minimum for this specific ring size. For very flexible acyclic molecules (>10 rotatable bonds), increasing `n_steps` and `pool_max` helps ensure thorough exploration. - -### PRISM Pruner Settings - -`PrismConfig.energy_window_kcal` (default: 15.0) — Only compare conformers within this energy window during pruning. Should be ≥ your main `energy_window_kcal` setting to avoid pruning conformers that are still within your desired energy range. - -```python -from openconf import PrismConfig - -prism_config = PrismConfig(energy_window_kcal=20.0) - -config = ConformerConfig( - max_out=200, - energy_window_kcal=15.0, - prism_config=prism_config, -) -``` +Macrocycles (ring size ≥ 10): openconf uses ETKDGv3 macrocycle distance bounds at seed time, disables seed RMSD pruning (which otherwise collapses distinct puckers), applies crankshaft and kinematic ring-closure moves during exploration, and supports an amide-flip move for macrocycles containing in-ring amide bonds (both secondary and tertiary). On cycloalkanes C6–C14 the default config recovers the ETKDG+MMFF94s reference minimum within 1 kcal/mol. Cyclododecane is an edge case — its global [3333] basin occupies <1% of seeds and needs `n_seeds ≥ 300` for reliable recovery (default is ~164); override `n_seeds` if you need the exact global minimum for this specific ring size. The `"macrocycle"` preset is the recommended starting point: it enables low-mode following for Hessian-guided seeds and uses a 100 kcal/mol energy window to capture the full range of ring-pucker minima. For very flexible acyclic molecules (>10 rotatable bonds), increasing `n_steps` and `pool_max` helps ensure thorough exploration. ### Performance Tips diff --git a/openconf/config.py b/openconf/config.py index ff14ddb..7b0bcfe 100644 --- a/openconf/config.py +++ b/openconf/config.py @@ -5,7 +5,7 @@ from .tuning import get_default_move_probs -ConformerPreset = Literal["rapid", "ensemble", "spectroscopic", "docking", "analogue"] +ConformerPreset = Literal["rapid", "ensemble", "spectroscopic", "docking", "analogue", "macrocycle"] _SUPPORTED_MOVE_TYPES = frozenset( { "single_rotor", @@ -15,6 +15,7 @@ "ring_flip", "crankshaft", "ring_kic", + "amide_flip", } ) _SUPPORTED_PARENT_STRATEGIES = frozenset({"softmax", "uniform", "best"}) @@ -104,6 +105,9 @@ class ConformerConfig: random_seed: Random seed for reproducibility. num_threads: thread count for parallel operations. use_heavy_atoms_only: Use only heavy atoms for RMSD calculations. + prism_max_deviation: MoI similarity threshold for PRISM deduplication. Two conformers are + considered duplicates if all three principal moments differ by less than this fraction. + Smaller values mean less aggressive pruning (more conformers survive). Default 0.01 (1%). clash_threshold: Distance threshold for clash detection (Angstroms). fast_minimization_iters: Iterations for quick minimization. max_minimization_iters: Maximum iterations for minimization. @@ -192,6 +196,40 @@ class ConformerConfig: to try for clash-filtered torsion moves. The candidate with the lowest clash score is kept, reducing wasted minimizations on crowded intermediates. Set to 1 to recover single-proposal behavior. + use_low_mode_following: when True, supplement ETKDG seeds with + conformers generated by displacing each minimized seed along its + low-frequency Hessian eigenvectors and re-minimizing. Adds coverage + for collective ring-puckering and correlated torsion motions that + independent single-rotor moves miss. Disabled by default because + Hessian evaluation costs 6N MMFF gradient calls per seed (N = atom + count); enable explicitly for macrocycles or other challenging + systems where the extra sampling is warranted. + low_mode_eigenvalue_threshold: Hessian eigenvalue cutoff (kcal/mol/Ų) + below which eigenvectors are treated as conformationally soft and + followed. Lower values select only the very softest modes; higher + values include stiffer bending modes. Default 100.0 captures + torsional and ring-puckering modes while excluding bond stretches. + low_mode_max_modes: maximum number of low-frequency eigenvectors to + follow per seed conformer. Caps the number of additional + minimizations regardless of how many modes satisfy + low_mode_eigenvalue_threshold. + low_mode_scan_step_size: distance advanced per scan step in Å + (3N Euclidean norm of the displacement vector). Smaller values + give finer resolution of the stopping point at the cost of more + force-field evaluations per scan. Default 0.25 Å. + low_mode_scan_energy_threshold: maximum per-step energy increase + (kcal/mol) before scanning along a mode stops. The default + (~2390 kcal/mol, equivalent to the paper's 10 000 kJ/mol) allows + the scan to pass through conformational barriers and terminates + only at severe steric clashes. + low_mode_scan_max_steps: upper bound on scan steps per direction; + acts as a safety cap on total displacement regardless of the + energy criterion. + low_mode_n_source_seeds: number of minimized ETKDG seeds to compute + Hessians for, selected by ascending energy. Remaining seeds are + skipped for low mode following but still enter the MCMM pool + normally. Default 1 (only the lowest-energy seed) keeps the + Hessian cost fixed regardless of total seed count. """ max_out: int = 200 @@ -207,6 +245,7 @@ class ConformerConfig: random_seed: int | None = None num_threads: int = 0 use_heavy_atoms_only: bool = True + prism_max_deviation: float = 0.01 constraint_spec: ConstraintSpec | None = None clash_threshold: float = 1.5 fast_minimization_iters: int = 20 @@ -232,6 +271,13 @@ class ConformerConfig: auto_tune_large_flexible: bool = True collect_stats: bool = False torsion_multitry_attempts: int = 4 + use_low_mode_following: bool = False + low_mode_eigenvalue_threshold: float = 100.0 + low_mode_max_modes: int = 5 + low_mode_scan_step_size: float = 0.25 + low_mode_scan_energy_threshold: float = 2390.0 + low_mode_scan_max_steps: int = 10 + low_mode_n_source_seeds: int = 1 def __post_init__(self) -> None: _require_at_least("max_out", self.max_out, 1) @@ -257,6 +303,12 @@ def __post_init__(self) -> None: _require_fraction("adapt_decay", self.adapt_decay) _require_at_least("patience", self.patience, 0) _require_at_least("torsion_multitry_attempts", self.torsion_multitry_attempts, 1) + _require_greater_than("low_mode_eigenvalue_threshold", self.low_mode_eigenvalue_threshold, 0.0) + _require_at_least("low_mode_max_modes", self.low_mode_max_modes, 1) + _require_greater_than("low_mode_scan_step_size", self.low_mode_scan_step_size, 0.0) + _require_greater_than("low_mode_scan_energy_threshold", self.low_mode_scan_energy_threshold, 0.0) + _require_at_least("low_mode_scan_max_steps", self.low_mode_scan_max_steps, 1) + _require_at_least("low_mode_n_source_seeds", self.low_mode_n_source_seeds, 1) _validate_move_probs(self.move_probs) if self.minimizer != "rdkit_mmff": @@ -305,9 +357,14 @@ def preset_config(preset: ConformerPreset) -> ConformerConfig: `constraint_spec` automatically. 50 conformers, full refinement, softmax parent strategy to stay near the constrained energy basin. + - `"macrocycle"` — Macrocyclic ring systems (ring size ≥ 10). Wide 100 + kcal/mol energy window to capture the full range of ring-pucker minima, + and low-mode following enabled to supplement ETKDG seeds with + Hessian-guided scan points along collective ring-puckering modes. + Args: preset: one of `"rapid"`, `"ensemble"`, `"spectroscopic"`, - `"docking"`, `"analogue"`. + `"docking"`, `"analogue"`, `"macrocycle"`. Returns: Configuration for requested use case @@ -322,6 +379,9 @@ def preset_config(preset: ConformerPreset) -> ConformerConfig: >>> config = preset_config("rapid") >>> config.do_final_refine False + >>> config = preset_config("macrocycle") + >>> config.energy_window_kcal + 100.0 """ match preset: case "rapid": @@ -393,7 +453,21 @@ def preset_config(preset: ConformerPreset) -> ConformerConfig: final_select="diverse", adaptive_moves=True, ) + case "macrocycle": + return ConformerConfig( + max_out=200, + n_steps=500, + energy_window_kcal=100.0, + use_low_mode_following=True, + seed_n_per_rotor=3, + seed_prune_rms_thresh=1.0, + do_final_refine=True, + minimize_batch_size=8, + parent_strategy="softmax", + final_select="diverse", + ) case _: raise ValueError( - f"Unknown preset {preset!r}. Choose from: 'rapid', 'ensemble', 'spectroscopic', 'docking', 'analogue'." + f"Unknown preset {preset!r}. " + "Choose from: 'rapid', 'ensemble', 'spectroscopic', 'docking', 'analogue', 'macrocycle'." ) diff --git a/openconf/data/runtime_tuning.json b/openconf/data/runtime_tuning.json index f213053..b56c0a2 100644 --- a/openconf/data/runtime_tuning.json +++ b/openconf/data/runtime_tuning.json @@ -47,26 +47,30 @@ "global_shake": 0.08, "ring_flip": 0.1, "crankshaft": 0.06, - "ring_kic": 0.06 + "ring_kic": 0.06, + "amide_flip": 0.06 }, "forced_periodic_move": "global_shake", "constrained_fallback_move": "single_rotor", "constrained_disabled_moves": [ "global_shake", "crankshaft", - "ring_kic" + "ring_kic", + "amide_flip" ], "availability_fallbacks": { "ring_flip": "single_rotor", "crankshaft": "single_rotor", - "ring_kic": "crankshaft" + "ring_kic": "crankshaft", + "amide_flip": "single_rotor" } }, "clash_check": { "exempt_moves": [ "ring_flip", "crankshaft", - "ring_kic" + "ring_kic", + "amide_flip" ] }, "low_flex_path": { diff --git a/openconf/dedupe.py b/openconf/dedupe.py index 44b7c0e..e57f2c4 100644 --- a/openconf/dedupe.py +++ b/openconf/dedupe.py @@ -47,6 +47,7 @@ def prism_dedupe( mol: Chem.Mol, conf_ids: list[int], use_heavy_atoms_only: bool = True, + max_deviation: float = 0.01, ) -> list[int]: """Deduplicate conformers using PRISM Pruner. @@ -54,6 +55,8 @@ def prism_dedupe( mol: molecule with conformers conf_ids: conformer IDs to process use_heavy_atoms_only: use only heavy atoms for comparison + max_deviation: MoI similarity threshold; conformers with all three principal + moments within this fraction of each other are considered duplicates. Returns: Identifiers to keep @@ -62,5 +65,5 @@ def prism_dedupe( return conf_ids coords, atoms = _mol_to_arrays(mol, conf_ids, use_heavy_atoms_only) - _, mask = prune_by_moment_of_inertia(coords, atoms) + _, mask = prune_by_moment_of_inertia(coords, atoms, max_deviation=max_deviation) return [conf_id for conf_id, keep in zip(conf_ids, mask, strict=True) if keep] diff --git a/openconf/pool.py b/openconf/pool.py index 07c5e01..3b7aa39 100644 --- a/openconf/pool.py +++ b/openconf/pool.py @@ -366,6 +366,7 @@ def dedupe(self) -> int: self.mol, conf_ids, use_heavy_atoms_only=self.config.use_heavy_atoms_only, + max_deviation=self.config.prism_max_deviation, ) # Remove duplicates diff --git a/openconf/propose/__init__.py b/openconf/propose/__init__.py index 50a1852..8ad12ee 100644 --- a/openconf/propose/__init__.py +++ b/openconf/propose/__init__.py @@ -1,5 +1,6 @@ """Conformer proposal strategies.""" from .hybrid import HybridProposer, run_hybrid_generation +from .low_mode import generate_low_mode_seeds -__all__ = ["HybridProposer", "run_hybrid_generation"] +__all__ = ["HybridProposer", "generate_low_mode_seeds", "run_hybrid_generation"] diff --git a/openconf/propose/amide_seeds.py b/openconf/propose/amide_seeds.py index 8b7ca5f..161b22d 100644 --- a/openconf/propose/amide_seeds.py +++ b/openconf/propose/amide_seeds.py @@ -12,6 +12,72 @@ _TORSION_WINDOW_DEG = 30.0 _TORSION_FORCE_K = 500.0 # kcal/mol/rad² _MIN_MACROCYCLE_SIZE = 10 +# Smallest ring in which an in-ring amide can plausibly invert cis/trans. Below +# this the ring geometry locks the amide and a 180° flip is unrecoverable. +_MIN_AMIDE_FLIP_RING_SIZE = 8 + + +def find_ring_amide_bonds( + mol: Chem.Mol, + min_ring_size: int = _MIN_AMIDE_FLIP_RING_SIZE, +) -> list[tuple[int, int, tuple[int, ...]]]: + """Return in-ring amide C-N bonds eligible for the amide-flip MC move. + + Locates non-aromatic amide bonds (a carbonyl carbon single-bonded to N, + both in the same ring) inside rings of at least ``min_ring_size`` atoms, + where a 180° rotation about the C-N bond can meaningfully toggle the amide + between cis and trans. Both secondary (N-H) and tertiary amides qualify — + unlike the cis/trans seed enumeration, which targets only tertiary amides. + + Each amide bond is reported once, associated with the largest ring that + contains it (the macrocycle, for fused systems), so the flip rotates the + most flexible arc available. + + Args: + mol: molecule to scan + min_ring_size: smallest ring size considered flippable + + Returns: + Tuples of (carbonyl_carbon_idx, nitrogen_idx, ring_atom_indices), where + ring_atom_indices is the ordered atom ring containing the bond + """ + ring_info = mol.GetRingInfo() + rings = sorted((r for r in ring_info.AtomRings() if len(r) >= min_ring_size), key=len, reverse=True) + + results: list[tuple[int, int, tuple[int, ...]]] = [] + seen: set[tuple[int, int]] = set() + + for ring in rings: + ring_set = set(ring) + for n_idx in ring: + n_atom = mol.GetAtomWithIdx(n_idx) + if n_atom.GetAtomicNum() != 7 or n_atom.GetIsAromatic(): + continue + for bond in n_atom.GetBonds(): + if not bond.IsInRing(): + continue + c_atom = bond.GetOtherAtom(n_atom) + c_idx = c_atom.GetIdx() + if c_atom.GetAtomicNum() != 6 or c_idx not in ring_set or c_atom.GetIsAromatic(): + continue + + bond_key = (min(n_idx, c_idx), max(n_idx, c_idx)) + if bond_key in seen: + continue + + has_carbonyl = any( + cb.GetBondTypeAsDouble() >= 1.9 + and cb.GetOtherAtom(c_atom).GetAtomicNum() == 8 + and not cb.GetOtherAtom(c_atom).GetIsAromatic() + for cb in c_atom.GetBonds() + ) + if not has_carbonyl: + continue + + seen.add(bond_key) + results.append((c_idx, n_idx, tuple(ring))) + + return results def find_ring_tertiary_amide_dihedrals(mol: Chem.Mol) -> list[tuple[int, int, int, int]]: diff --git a/openconf/propose/hybrid.py b/openconf/propose/hybrid.py index 7df4e3e..426f3f6 100644 --- a/openconf/propose/hybrid.py +++ b/openconf/propose/hybrid.py @@ -28,6 +28,7 @@ ) from .amide_seeds import find_ring_tertiary_amide_dihedrals, generate_amide_variant_seeds from .candidates import CandidateBatchWorkspace, ClashChecker, build_nonbonded_mask +from .low_mode import generate_low_mode_seeds from .moves import MoveExecutor from .seeding import ( SeedPlan, @@ -510,6 +511,7 @@ def _select_move_type(self, step: int) -> str: has_ring_flips=bool(self.rotor_model.ring_flips), has_crankshaft=bool(self._moves.crankable_rings), has_kic=bool(self._moves.macro_kic_data), + has_amide=bool(self._moves.amide_flips), has_rotors=bool(self.rotor_model.rotors), ) @@ -964,6 +966,41 @@ def run_hybrid_generation( for conf_id, energy in proposer.seed_from_input_conformer(input_positions): pool.insert(conf_id, energy, source="seed_input") + if effective_config.use_low_mode_following and seeds: + lm_ff_props = getattr(proposer.fast_minimizer, "_mmff_props", None) + if lm_ff_props is not None: + lm_start = time.perf_counter() + n_lm_accepted = 0 + source_seeds = sorted(seeds, key=lambda x: x[1])[: effective_config.low_mode_n_source_seeds] + # Materialize the low-mode children for every source seed before any + # pool insertion. The source seeds are themselves pool members, so + # inserting one seed's children could evict a later, not-yet-processed + # source seed and leave generate_low_mode_seeds referencing a freed + # conformer (Bad Conformer Id, or a hard abort during MMFF setup). + lm_children: list[tuple[int, float]] = [] + for seed_conf_id, _ in source_seeds: + lm_children.extend( + generate_low_mode_seeds( + mol, + lm_ff_props, + seed_conf_id, + proposer.fast_minimizer, + eigenvalue_threshold=effective_config.low_mode_eigenvalue_threshold, + max_modes=effective_config.low_mode_max_modes, + scan_step_size=effective_config.low_mode_scan_step_size, + scan_energy_threshold=effective_config.low_mode_scan_energy_threshold, + scan_max_steps=effective_config.low_mode_scan_max_steps, + ) + ) + for lm_conf_id, lm_energy in lm_children: + if pool.insert(lm_conf_id, lm_energy, source="low_mode"): + n_lm_accepted += 1 + else: + mol.RemoveConformer(lm_conf_id) + if stats: + stats["low_mode_time_s"] = time.perf_counter() - lm_start + stats["n_low_mode_seeds"] = n_lm_accepted + batch_size = effective_config.minimize_batch_size step = 0 stagnation = 0 diff --git a/openconf/propose/low_mode.py b/openconf/propose/low_mode.py new file mode 100644 index 0000000..ff172a2 --- /dev/null +++ b/openconf/propose/low_mode.py @@ -0,0 +1,296 @@ +"""Low mode following for conformer seeding. + +Numerically evaluates Hessians via gradient finite differences and scans +minimized geometries along low-eigenvalue eigenvectors to generate structurally +diverse starting points. Most effective for macrocycles and correlated flexible +systems where independent torsion moves fail to sample collective soft motions. +""" + +from typing import TYPE_CHECKING + +import numpy as np +from rdkit import Chem +from rdkit.Chem import AllChem + +if TYPE_CHECKING: + from ..relax import Minimizer + +_N_RIGID_MODES: int = 6 +_RIGID_BODY_GAP_WINDOW: int = _N_RIGID_MODES + 2 +_DEFAULT_FD_STEP: float = 0.005 +_DEFAULT_EIGENVALUE_THRESHOLD: float = 100.0 +_DEFAULT_MAX_MODES: int = 5 +_DEFAULT_SCAN_STEP_SIZE: float = 0.25 +_DEFAULT_SCAN_ENERGY_THRESHOLD: float = 2390.0 # kcal/mol ≈ 10 000 kJ/mol (paper value) +_DEFAULT_SCAN_MAX_STEPS: int = 10 + + +def _compute_hessian( + mol: Chem.Mol, + ff_props: object, + conf_id: int, + step: float = _DEFAULT_FD_STEP, +) -> np.ndarray: + """Compute numerical Hessian via central differences of MMFF gradients. + + Perturbs each Cartesian coordinate of each atom by ±step, evaluates + the MMFF gradient at each displaced geometry, and assembles the + 3N×3N Hessian matrix. Conformer positions are restored on exit. + + Args: + mol: molecule containing the conformer + ff_props: pre-prepared MMFFMoleculeProperties used to build force fields + at displaced geometries + conf_id: ID of the conformer at which to evaluate the Hessian; + should be at or near a local minimum for meaningful low modes + step: finite difference displacement in Å + + Returns: + Symmetric 3N×3N Hessian matrix in kcal/mol/Ų + """ + conf = mol.GetConformer(conf_id) + n_atoms = mol.GetNumAtoms() + n_dof = 3 * n_atoms + pos0 = conf.GetPositions().copy() + hessian = np.zeros((n_dof, n_dof)) + + for i in range(n_dof): + atom_i = i // 3 + coord_i = i % 3 + original = pos0[atom_i].copy() + + fwd = original.copy() + fwd[coord_i] += step + conf.SetAtomPosition(atom_i, fwd.tolist()) + ff_fwd = AllChem.MMFFGetMoleculeForceField(mol, ff_props, confId=conf_id) + if ff_fwd is None: + conf.SetAtomPosition(atom_i, original.tolist()) + continue + grad_plus = np.array(ff_fwd.CalcGrad()) + + bwd = original.copy() + bwd[coord_i] -= step + conf.SetAtomPosition(atom_i, bwd.tolist()) + ff_bwd = AllChem.MMFFGetMoleculeForceField(mol, ff_props, confId=conf_id) + if ff_bwd is None: + conf.SetAtomPosition(atom_i, original.tolist()) + continue + grad_minus = np.array(ff_bwd.CalcGrad()) + + conf.SetAtomPosition(atom_i, original.tolist()) + hessian[i] = (grad_plus - grad_minus) / (2.0 * step) + + return (hessian + hessian.T) * 0.5 + + +def _select_low_modes( + hessian: np.ndarray, + eigenvalue_threshold: float, + max_modes: int, +) -> np.ndarray: + """Select low-frequency eigenvectors of a Hessian. + + Diagonalizes the Hessian, locates the rigid-body / conformational boundary + via the largest eigenvalue gap within the first _RIGID_BODY_GAP_WINDOW + entries, skips all modes up to that boundary, and returns eigenvectors + whose eigenvalues are below eigenvalue_threshold, capped at max_modes. + + Using a gap rather than a fixed count of 6 handles linear molecules + (only 5 rigid-body modes) and numerical drift when minimization is not + fully converged. + + Args: + hessian: symmetric 3N×3N Hessian matrix + eigenvalue_threshold: upper eigenvalue bound (kcal/mol/Ų) for mode + selection; modes at or above this value are excluded + max_modes: maximum number of modes to return + + Returns: + Array of shape (3N, k) where k ≤ max_modes; columns are unit + eigenvectors in ascending eigenvalue order; shape (3N, 0) when + no conformational modes satisfy the threshold + """ + eigenvalues, eigenvectors = np.linalg.eigh(hessian) + window = min(_RIGID_BODY_GAP_WINDOW, len(eigenvalues) - 1) + n_skip = int(np.argmax(np.diff(eigenvalues[: window + 1]))) + 1 + conf_vals = eigenvalues[n_skip:] + conf_vecs = eigenvectors[:, n_skip:] + + mask = conf_vals < eigenvalue_threshold + if not np.any(mask): + return np.empty((hessian.shape[0], 0)) + return conf_vecs[:, mask][:, :max_modes] + + +def _scan_along_mode( + mol: Chem.Mol, + ff_props: object, + start_conf_id: int, + direction: np.ndarray, + step_size: float, + energy_threshold: float, + max_steps: int, +) -> np.ndarray: + """Scan from a minimized conformer along a unit direction vector. + + Takes discrete steps of step_size Å in direction, evaluating the MMFF + energy after each step and stopping as soon as the per-step energy increase + exceeds energy_threshold. Returns the positions from the last accepted step, + or the starting positions when the first step already exceeds the threshold. + + A temporary conformer is created and removed; the start conformer is + never modified. + + Args: + mol: molecule to scan; receives a temporary conformer during the call + ff_props: pre-prepared MMFFMoleculeProperties for energy evaluation + start_conf_id: ID of the minimized starting conformer + direction: unit displacement direction of shape (n_atoms, 3) in Å; + each step moves the geometry by step_size × direction in 3N space + step_size: distance (Å) to move in 3N Cartesian space per step + energy_threshold: maximum allowed per-step energy increase (kcal/mol); + scanning stops when ΔE in a single step exceeds this value + max_steps: maximum number of steps regardless of energy criterion + + Returns: + Positions of shape (n_atoms, 3) at the last accepted scan point + """ + n_atoms = mol.GetNumAtoms() + src_conf = mol.GetConformer(start_conf_id) + start_pos = src_conf.GetPositions().copy() + + ff0 = AllChem.MMFFGetMoleculeForceField(mol, ff_props, confId=start_conf_id) + if ff0 is None: + return start_pos + prev_energy = float(ff0.CalcEnergy()) + + working_conf = Chem.Conformer(src_conf) + working_id = mol.AddConformer(working_conf, assignId=True) + + current_pos = start_pos.copy() + accepted_pos = start_pos.copy() + + try: + for _ in range(max_steps): + current_pos = current_pos + step_size * direction + working = mol.GetConformer(working_id) + for i in range(n_atoms): + working.SetAtomPosition(i, current_pos[i].tolist()) + + ff = AllChem.MMFFGetMoleculeForceField(mol, ff_props, confId=working_id) + if ff is None: + break + + curr_energy = float(ff.CalcEnergy()) + if curr_energy - prev_energy > energy_threshold: + break + + accepted_pos = current_pos.copy() + prev_energy = curr_energy + finally: + mol.RemoveConformer(working_id) + + return accepted_pos + + +def generate_low_mode_seeds( + mol: Chem.Mol, + ff_props: object, + conf_id: int, + minimizer: "Minimizer", + *, + eigenvalue_threshold: float = _DEFAULT_EIGENVALUE_THRESHOLD, + max_modes: int = _DEFAULT_MAX_MODES, + scan_step_size: float = _DEFAULT_SCAN_STEP_SIZE, + scan_energy_threshold: float = _DEFAULT_SCAN_ENERGY_THRESHOLD, + scan_max_steps: int = _DEFAULT_SCAN_MAX_STEPS, + fd_step: float = _DEFAULT_FD_STEP, +) -> list[tuple[int, float]]: + """Generate conformers by scanning along low-frequency Hessian eigenvectors. + + Numerically evaluates the MMFF Hessian at the given minimized conformer, + identifies eigenvectors with eigenvalues below eigenvalue_threshold (soft + conformational modes), and for each such mode scans in both the positive and + negative directions using discrete steps of scan_step_size Å. Scanning along + a direction stops when the per-step energy increase exceeds scan_energy_threshold + or scan_max_steps is reached. Each scan endpoint is then minimized to a local + minimum. + + This mirrors the LMOD procedure of Kolossváry & Guida (JACS 1996): the scan + naturally traverses soft conformational barriers and terminates at the onset + of severe steric clashes, placing the starting geometry for minimization on + the far side of a barrier. + + New conformers are added to mol in place. Callers are responsible for + removing conformers that are not kept (e.g., pool rejects). + + Note: + Hessian evaluation requires 6N MMFF force-field constructions where N + is the atom count. This is the dominant cost; each scan step adds two + further force-field constructions (one per direction). + + Args: + mol: molecule containing the conformer; receives new conformers in place + ff_props: pre-prepared MMFFMoleculeProperties used to build force fields + conf_id: ID of a minimized conformer to compute modes from + minimizer: minimizer applied to each scan endpoint + eigenvalue_threshold: Hessian eigenvalue cutoff (kcal/mol/Ų); modes + below this value are treated as conformationally soft + max_modes: maximum number of low modes to scan per conformer + scan_step_size: distance to advance per scan step in Å (3N Euclidean + norm of the displacement vector); smaller values give finer + resolution of the stopping point + scan_energy_threshold: maximum per-step energy increase (kcal/mol) + before scanning stops; the default (~2390 kcal/mol ≈ 10 000 kJ/mol) + follows the paper and effectively allows the scan to pass through + conformational barriers, stopping only at severe steric clashes + scan_max_steps: upper bound on scan steps per direction regardless of + the energy criterion; acts as a safety cap on total displacement + fd_step: finite difference step size for numerical Hessian in Å + + Returns: + Pairs of (conformer_id, energy_kcal_mol) for each successfully + minimized scan endpoint; at most 2 results per mode (the two scan + senses); empty when no low modes satisfy the threshold or all + minimizations fail + """ + hessian = _compute_hessian(mol, ff_props, conf_id, fd_step) + low_vecs = _select_low_modes(hessian, eigenvalue_threshold, max_modes) + if low_vecs.shape[1] == 0: + return [] + + n_atoms = mol.GetNumAtoms() + start_pos = mol.GetConformer(conf_id).GetPositions().copy() + + results: list[tuple[int, float]] = [] + for col in range(low_vecs.shape[1]): + direction = low_vecs[:, col].reshape(n_atoms, 3) + + for sign in (1.0, -1.0): + final_pos = _scan_along_mode( + mol, + ff_props, + conf_id, + sign * direction, + scan_step_size, + scan_energy_threshold, + scan_max_steps, + ) + + if np.allclose(final_pos, start_pos, atol=1e-8): + continue + + new_conf = Chem.Conformer(mol.GetConformer(conf_id)) + new_conf_id = mol.AddConformer(new_conf, assignId=True) + displaced = mol.GetConformer(new_conf_id) + for i in range(n_atoms): + displaced.SetAtomPosition(i, final_pos[i].tolist()) + + energy = minimizer.minimize(mol, new_conf_id) + if not np.isfinite(energy): + mol.RemoveConformer(new_conf_id) + continue + + results.append((new_conf_id, energy)) + + return results diff --git a/openconf/propose/moves.py b/openconf/propose/moves.py index b061f64..c604a70 100644 --- a/openconf/propose/moves.py +++ b/openconf/propose/moves.py @@ -15,6 +15,7 @@ conformer_matches_specified_stereochemistry, specified_stereochemistry, ) +from .amide_seeds import find_ring_amide_bonds if TYPE_CHECKING: from collections.abc import Callable @@ -109,6 +110,7 @@ def __init__( self._ring_flip_allowed_arcs(ring_flip.ring_atoms, ring_flip.junction_atoms) for ring_flip in rotor_model.ring_flips ] + self.amide_flips = self._build_amide_flips(mol) self.operators: dict[str, MoveOperator] = { "single_rotor": self.apply_single_rotor_move, "multi_rotor": self.apply_multi_rotor_move, @@ -117,6 +119,7 @@ def __init__( "ring_flip": self.apply_ring_flip_move, "crankshaft": self.apply_crankshaft_move, "ring_kic": self.apply_ring_kic_move, + "amide_flip": self.apply_amide_flip_move, } @staticmethod @@ -179,6 +182,30 @@ def _build_crankable_rings( crankable.append((tuple(ring), subtrees)) return crankable + def _build_amide_flips(self, mol: Chem.Mol) -> list[tuple[int, int, np.ndarray]]: + """Precompute in-ring amide-flip rotation data. + + For each flippable in-ring amide bond, the moving set is the union of + the substituent subtrees of every ring atom except the carbonyl carbon + and nitrogen themselves. Those two atoms anchor the rotation axis, so + rotating the rest of the ring 180° about the C-N bond inverts the amide + cis/trans relationship while preserving every ring bond length exactly + (only the bond angles at C and N deform, which MMFF relaxes). + """ + amide_flips: list[tuple[int, int, np.ndarray]] = [] + for c_idx, n_idx, ring_atoms in find_ring_amide_bonds(mol): + ring_set = frozenset(ring_atoms) + moving: set[int] = set() + for ring_atom in ring_atoms: + if ring_atom in (c_idx, n_idx): + continue + moving |= self._compute_substituent_atoms(mol, ring_set, ring_atom) + if not moving: + continue + moving_arr = np.fromiter(moving, dtype=np.int64, count=len(moving)) + amide_flips.append((c_idx, n_idx, moving_arr)) + return amide_flips + def sample_angle(self, rotor_idx: int) -> float: """Sample a torsion angle for a rotor.""" angles_arr, weights_arr = self._rotor_angles[rotor_idx] @@ -276,6 +303,43 @@ def apply_crankshaft_move(self, conf_id: int) -> None: ) conf.SetPositions(all_pos) + def apply_amide_flip_move(self, conf_id: int) -> None: + """Invert an in-ring amide by rotating its ring arc 180° about the C-N bond. + + Picks a flippable in-ring amide and rotates the precomputed moving arc + about the carbonyl-carbon→nitrogen axis by ~180° (with the usual torsion + jitter). Because both anchors lie on the rotation axis, ring bond lengths + are preserved and only the angles at C and N strain; the subsequent MMFF + minimization relaxes them. Like the other ring moves this is exempt from + the static clash filter, since the flipped geometry is transiently + crowded before minimization. No-op when the molecule has no flippable + in-ring amides. + """ + if not self.amide_flips: + return + + c_idx, n_idx, moving_arr = random.choice(self.amide_flips) + conf = self.mol.GetConformer(conf_id) + all_pos = conf.GetPositions() + axis_origin = all_pos[c_idx] + axis_vec = all_pos[n_idx] - axis_origin + axis_norm = float(np.linalg.norm(axis_vec)) + if axis_norm < 1e-6: + return + axis_unit = axis_vec / axis_norm + + theta = np.deg2rad(180.0 + np.random.normal(0.0, self.config.torsion_jitter_deg)) + _rodrigues_rotate( + all_pos, + moving_arr, + axis_origin, + float(axis_unit[0]), + float(axis_unit[1]), + float(axis_unit[2]), + theta, + ) + conf.SetPositions(all_pos) + def apply_ring_flip_move(self, conf_id: int) -> None: """Apply ring flip or stereo-preserving ring pucker move.""" if not self.rotor_model.ring_flips: diff --git a/openconf/propose/stats.py b/openconf/propose/stats.py index 54dae54..c8196c8 100644 --- a/openconf/propose/stats.py +++ b/openconf/propose/stats.py @@ -38,9 +38,11 @@ def new_generation_stats() -> dict[str, GenerationStat]: "n_dedupe_calls": 0, "n_dedupe_removed": 0, "n_final_selected": 0, + "n_low_mode_seeds": 0, "seed_time_s": 0.0, "seed_embedding_time_s": 0.0, "seed_minimization_time_s": 0.0, + "low_mode_time_s": 0.0, "proposal_stage_time_s": 0.0, "parent_selection_time_s": 0.0, "move_selection_time_s": 0.0, diff --git a/openconf/tuning.py b/openconf/tuning.py index 8a64799..0b73a28 100644 --- a/openconf/tuning.py +++ b/openconf/tuning.py @@ -208,6 +208,7 @@ def resolve_move_probabilities( has_ring_flips: bool, has_crankshaft: bool, has_kic: bool = False, + has_amide: bool = False, has_rotors: bool = True, ) -> dict[str, float]: """Return move probabilities adjusted for current availability constraints.""" @@ -224,6 +225,7 @@ def resolve_move_probabilities( "ring_flip": has_ring_flips, "crankshaft": has_crankshaft, "ring_kic": has_kic, + "amide_flip": has_amide, } for move_type, is_available in availability.items(): if is_available or move_type not in probs: diff --git a/tests/test_low_mode.py b/tests/test_low_mode.py new file mode 100644 index 0000000..7aada98 --- /dev/null +++ b/tests/test_low_mode.py @@ -0,0 +1,457 @@ +"""Unit and integration tests for low mode following conformer seeding.""" + +import numpy as np +import pytest +from rdkit import Chem +from rdkit.Chem import AllChem + +from openconf import ConformerConfig, generate_conformers, preset_config +from openconf.perceive import prepare_molecule +from openconf.propose.low_mode import ( + _DEFAULT_EIGENVALUE_THRESHOLD, + _DEFAULT_MAX_MODES, + _DEFAULT_SCAN_ENERGY_THRESHOLD, + _DEFAULT_SCAN_MAX_STEPS, + _DEFAULT_SCAN_STEP_SIZE, + _N_RIGID_MODES, + _compute_hessian, + _scan_along_mode, + _select_low_modes, + generate_low_mode_seeds, +) +from openconf.relax import RDKitMMFFMinimizer + + +def _make_minimized_mol(smiles: str, seed: int = 1) -> tuple[Chem.Mol, object]: + """Return (mol, ff_props) for a fully H-added, MMFF-minimized molecule.""" + mol = prepare_molecule(Chem.MolFromSmiles(smiles)) + AllChem.EmbedMolecule(mol, randomSeed=seed) + props = AllChem.MMFFGetMoleculeProperties(mol, mmffVariant="MMFF94s") + assert props is not None, f"MMFF typing failed for {smiles}" + ff = AllChem.MMFFGetMoleculeForceField(mol, props, confId=0) + assert ff is not None + ff.Minimize(maxIts=500) + return mol, props + + +def _make_fast_minimizer(mol: Chem.Mol) -> RDKitMMFFMinimizer: + """Return an RDKitMMFFMinimizer prepared for mol.""" + m = RDKitMMFFMinimizer(max_iters=50) + m.prepare(mol) + return m + + +# --------------------------------------------------------------------------- +# _compute_hessian +# --------------------------------------------------------------------------- + + +def test_hessian_shape(): + """Hessian must have shape (3N, 3N) for N-atom molecule.""" + mol, props = _make_minimized_mol("CCC") + n = mol.GetNumAtoms() + H = _compute_hessian(mol, props, conf_id=0) + assert H.shape == (3 * n, 3 * n) + + +def test_hessian_is_symmetric(): + """Hessian must be symmetric to within numerical noise.""" + mol, props = _make_minimized_mol("CCC") + H = _compute_hessian(mol, props, conf_id=0) + assert np.allclose(H, H.T, atol=1e-6), "Hessian is not symmetric" + + +def test_hessian_restores_positions(): + """Conformer coordinates must be unchanged after Hessian computation.""" + mol, props = _make_minimized_mol("CCC") + pos_before = mol.GetConformer(0).GetPositions().copy() + _compute_hessian(mol, props, conf_id=0) + pos_after = mol.GetConformer(0).GetPositions() + assert np.allclose(pos_before, pos_after, atol=1e-10), "Positions were not restored" + + +def test_hessian_positive_eigenvalues_at_minimum(): + """Non-rigid eigenvalues should be positive at a converged minimum.""" + mol, props = _make_minimized_mol("CCC") + H = _compute_hessian(mol, props, conf_id=0) + eigenvalues = np.linalg.eigvalsh(H) + conformational = eigenvalues[_N_RIGID_MODES:] + neg = conformational[conformational < -0.5] + assert np.all(conformational > -0.5), f"Unexpected large negative eigenvalues: {neg}" + + +# --------------------------------------------------------------------------- +# _select_low_modes +# --------------------------------------------------------------------------- + + +def test_select_low_modes_returns_correct_shape(): + """Selected modes must form a (3N, k) matrix with k ≤ max_modes.""" + mol, props = _make_minimized_mol("CCC") + H = _compute_hessian(mol, props, conf_id=0) + n_dof = 3 * mol.GetNumAtoms() + vecs = _select_low_modes(H, _DEFAULT_EIGENVALUE_THRESHOLD, _DEFAULT_MAX_MODES) + assert vecs.shape[0] == n_dof + assert vecs.shape[1] <= _DEFAULT_MAX_MODES + + +def test_select_low_modes_empty_when_threshold_zero(): + """No modes returned when threshold is effectively zero.""" + mol, props = _make_minimized_mol("CCC") + H = _compute_hessian(mol, props, conf_id=0) + vecs = _select_low_modes(H, eigenvalue_threshold=0.0, max_modes=10) + assert vecs.shape[1] == 0 + + +def test_select_low_modes_respects_max_modes(): + """Result must contain at most max_modes columns.""" + mol, props = _make_minimized_mol("CCCCC") + H = _compute_hessian(mol, props, conf_id=0) + for cap in [1, 2, 3]: + vecs = _select_low_modes(H, eigenvalue_threshold=500.0, max_modes=cap) + assert vecs.shape[1] <= cap, f"Got {vecs.shape[1]} modes with cap={cap}" + + +def test_select_low_modes_columns_are_unit_vectors(): + """Returned eigenvectors must have unit norm (they come from eigh).""" + mol, props = _make_minimized_mol("CCC") + H = _compute_hessian(mol, props, conf_id=0) + vecs = _select_low_modes(H, _DEFAULT_EIGENVALUE_THRESHOLD, _DEFAULT_MAX_MODES) + if vecs.shape[1] == 0: + pytest.skip("No low modes found — threshold may need adjustment") + norms = np.linalg.norm(vecs, axis=0) + assert np.allclose(norms, 1.0, atol=1e-10), f"Mode norms not 1.0: {norms}" + + +# --------------------------------------------------------------------------- +# _scan_along_mode +# --------------------------------------------------------------------------- + + +def test_scan_along_mode_moves_from_start(): + """Scan must return positions that differ from the starting geometry.""" + mol, props = _make_minimized_mol("CCC") + H = _compute_hessian(mol, props, conf_id=0) + vecs = _select_low_modes(H, _DEFAULT_EIGENVALUE_THRESHOLD, 1) + if vecs.shape[1] == 0: + pytest.skip("No low modes found") + + n_atoms = mol.GetNumAtoms() + direction = vecs[:, 0].reshape(n_atoms, 3) + start_pos = mol.GetConformer(0).GetPositions().copy() + + final_pos = _scan_along_mode( + mol, props, 0, direction, _DEFAULT_SCAN_STEP_SIZE, + _DEFAULT_SCAN_ENERGY_THRESHOLD, _DEFAULT_SCAN_MAX_STEPS, + ) + assert not np.allclose(final_pos, start_pos), "Scan returned starting positions" + + +def test_scan_along_mode_restores_start_conformer(): + """The start conformer must be unchanged after scanning.""" + mol, props = _make_minimized_mol("CCC") + H = _compute_hessian(mol, props, conf_id=0) + vecs = _select_low_modes(H, _DEFAULT_EIGENVALUE_THRESHOLD, 1) + if vecs.shape[1] == 0: + pytest.skip("No low modes found") + + n_atoms = mol.GetNumAtoms() + direction = vecs[:, 0].reshape(n_atoms, 3) + pos_before = mol.GetConformer(0).GetPositions().copy() + n_confs_before = mol.GetNumConformers() + + _scan_along_mode( + mol, props, 0, direction, _DEFAULT_SCAN_STEP_SIZE, + _DEFAULT_SCAN_ENERGY_THRESHOLD, _DEFAULT_SCAN_MAX_STEPS, + ) + + assert mol.GetNumConformers() == n_confs_before, "Scan left a temporary conformer behind" + assert np.allclose(mol.GetConformer(0).GetPositions(), pos_before, atol=1e-10) + + +def test_scan_stops_at_energy_threshold(): + """With a near-zero energy threshold, scan must stop after the first step.""" + mol, props = _make_minimized_mol("CCC") + H = _compute_hessian(mol, props, conf_id=0) + vecs = _select_low_modes(H, _DEFAULT_EIGENVALUE_THRESHOLD, 1) + if vecs.shape[1] == 0: + pytest.skip("No low modes found") + + n_atoms = mol.GetNumAtoms() + direction = vecs[:, 0].reshape(n_atoms, 3) + start_pos = mol.GetConformer(0).GetPositions().copy() + + # Threshold of 0 means any energy increase (even numerical noise) stops the scan. + # The returned positions should still be the start positions (first step rejected). + final_pos = _scan_along_mode( + mol, props, 0, direction, _DEFAULT_SCAN_STEP_SIZE, + energy_threshold=0.0, max_steps=10, + ) + # Either no progress or exactly one step taken — at most a tiny displacement + displacement = float(np.linalg.norm(final_pos - start_pos)) + assert displacement < _DEFAULT_SCAN_STEP_SIZE + 1e-6 + + +# --------------------------------------------------------------------------- +# generate_low_mode_seeds +# --------------------------------------------------------------------------- + + +@pytest.fixture +def propane_mol_props() -> tuple[Chem.Mol, object]: + """Propane molecule with MMFF properties, minimized.""" + return _make_minimized_mol("CCC") + + +def test_generate_low_mode_seeds_returns_finite_energies( + propane_mol_props: tuple[Chem.Mol, object], +) -> None: + """All returned seeds must have finite MMFF energies.""" + mol, props = propane_mol_props + minimizer = _make_fast_minimizer(mol) + seeds = generate_low_mode_seeds(mol, props, conf_id=0, minimizer=minimizer) + for _, energy in seeds: + assert np.isfinite(energy), f"Non-finite energy in low mode seed: {energy}" + + +def test_generate_low_mode_seeds_adds_conformers_to_mol( + propane_mol_props: tuple[Chem.Mol, object], +) -> None: + """Each returned seed conformer must exist in the molecule.""" + mol, props = propane_mol_props + minimizer = _make_fast_minimizer(mol) + n_before = mol.GetNumConformers() + seeds = generate_low_mode_seeds(mol, props, conf_id=0, minimizer=minimizer) + existing_ids = {c.GetId() for c in mol.GetConformers()} + for conf_id, _ in seeds: + assert conf_id in existing_ids, f"Seed conf_id {conf_id} not found in mol" + assert mol.GetNumConformers() >= n_before + + +def test_generate_low_mode_seeds_empty_when_threshold_zero( + propane_mol_props: tuple[Chem.Mol, object], +) -> None: + """No seeds generated when eigenvalue threshold excludes all conformational modes.""" + mol, props = propane_mol_props + minimizer = _make_fast_minimizer(mol) + seeds = generate_low_mode_seeds( + mol, props, conf_id=0, minimizer=minimizer, eigenvalue_threshold=0.0 + ) + assert seeds == [] + + +def test_generate_low_mode_seeds_at_most_two_per_mode( + propane_mol_props: tuple[Chem.Mol, object], +) -> None: + """With max_modes=k, at most 2k seeds are returned (two directions per mode).""" + mol, props = propane_mol_props + minimizer = _make_fast_minimizer(mol) + for cap in [1, 2]: + n_before = mol.GetNumConformers() + seeds = generate_low_mode_seeds( + mol, props, conf_id=0, minimizer=minimizer, max_modes=cap + ) + assert len(seeds) <= 2 * cap, f"Got {len(seeds)} seeds with cap={cap}" + for conf_id, _ in seeds: + mol.RemoveConformer(conf_id) + assert mol.GetNumConformers() == n_before + + +def test_generate_low_mode_seeds_scans_both_directions() -> None: + """Both the + and − scan directions must be explored for each mode.""" + mol, props = _make_minimized_mol("CCCCCC") + minimizer = _make_fast_minimizer(mol) + + seeds = generate_low_mode_seeds( + mol, props, conf_id=0, minimizer=minimizer, max_modes=1 + ) + # For one mode with two directions we expect 0 or 2 seeds (not 1), + # because both directions are always attempted. + assert len(seeds) in (0, 2), f"Expected 0 or 2 seeds for 1 mode, got {len(seeds)}" + for conf_id, _ in seeds: + mol.RemoveConformer(conf_id) + + +def test_generate_low_mode_seeds_perturbs_geometry() -> None: + """Seeds must be displaced from the starting geometry after minimization.""" + mol, props = _make_minimized_mol("CCCCCC") + minimizer = _make_fast_minimizer(mol) + pos0 = mol.GetConformer(0).GetPositions().copy() + + seeds = generate_low_mode_seeds(mol, props, conf_id=0, minimizer=minimizer) + if not seeds: + pytest.skip("No low modes below threshold for this molecule") + + any_moved = False + for conf_id, _ in seeds: + pos_new = mol.GetConformer(conf_id).GetPositions() + rmsd = float(np.sqrt(np.mean(np.sum((pos_new - pos0) ** 2, axis=1)))) + if rmsd > 0.01: + any_moved = True + mol.RemoveConformer(conf_id) + + assert any_moved, "All low-mode seeds minimized back to the starting geometry" + + +# --------------------------------------------------------------------------- +# Config validation +# --------------------------------------------------------------------------- + + +def test_config_default_low_mode_off(): + """use_low_mode_following must be False by default.""" + assert ConformerConfig().use_low_mode_following is False + + +def test_config_low_mode_threshold_must_be_positive(): + """Negative or zero eigenvalue threshold must raise ValueError.""" + with pytest.raises(ValueError, match="low_mode_eigenvalue_threshold"): + ConformerConfig(low_mode_eigenvalue_threshold=0.0) + + +def test_config_low_mode_max_modes_must_be_at_least_one(): + """max_modes < 1 must raise ValueError.""" + with pytest.raises(ValueError, match="low_mode_max_modes"): + ConformerConfig(low_mode_max_modes=0) + + +def test_config_scan_step_size_must_be_positive(): + """Non-positive scan step size must raise ValueError.""" + with pytest.raises(ValueError, match="low_mode_scan_step_size"): + ConformerConfig(low_mode_scan_step_size=0.0) + + +def test_config_scan_energy_threshold_must_be_positive(): + """Non-positive energy threshold must raise ValueError.""" + with pytest.raises(ValueError, match="low_mode_scan_energy_threshold"): + ConformerConfig(low_mode_scan_energy_threshold=0.0) + + +def test_config_scan_max_steps_must_be_at_least_one(): + """scan_max_steps < 1 must raise ValueError.""" + with pytest.raises(ValueError, match="low_mode_scan_max_steps"): + ConformerConfig(low_mode_scan_max_steps=0) + + +def test_config_n_source_seeds_must_be_at_least_one(): + """low_mode_n_source_seeds < 1 must raise ValueError.""" + with pytest.raises(ValueError, match="low_mode_n_source_seeds"): + ConformerConfig(low_mode_n_source_seeds=0) + + +def test_n_source_seeds_limits_hessian_evaluations(): + """With n_source_seeds=1, only one seed triggers low mode scanning.""" + config = ConformerConfig( + use_low_mode_following=True, + low_mode_n_source_seeds=1, + n_seeds=5, + n_steps=10, + max_out=5, + random_seed=0, + collect_stats=True, + do_final_refine=False, + ) + ens = generate_conformers("CCCC", config=config) + stats = ens.generation_stats + # With n_source_seeds=1 and max_modes=5, at most 2*5=10 low mode seeds + # can be generated (two directions per mode); this is an upper bound check. + assert int(stats.get("n_low_mode_seeds", 0)) <= 10 + + +# --------------------------------------------------------------------------- +# Integration: stats populated when enabled +# --------------------------------------------------------------------------- + + +def test_low_mode_following_stat_populated(): + """With use_low_mode_following=True and collect_stats=True, stat keys appear.""" + config = ConformerConfig( + use_low_mode_following=True, + n_steps=20, + n_seeds=3, + max_out=5, + random_seed=42, + collect_stats=True, + do_final_refine=False, + ) + ens = generate_conformers("CCCC", config=config) + stats = ens.generation_stats + assert stats + assert "low_mode_time_s" in stats + assert "n_low_mode_seeds" in stats + assert float(stats["low_mode_time_s"]) >= 0.0 + assert int(stats["n_low_mode_seeds"]) >= 0 + + +def test_low_mode_following_survives_pool_overflow(): + """Source seeds must not be evicted before their low modes are scanned. + + With a small pool and several source seeds, the low-mode children can + overflow the pool partway through. Eviction of a not-yet-processed source + seed used to leave generate_low_mode_seeds referencing a freed conformer + (Bad Conformer Id or a hard MMFF abort). Generation must complete cleanly. + """ + config = ConformerConfig( + use_low_mode_following=True, + low_mode_n_source_seeds=5, + n_seeds=8, + n_steps=10, + max_out=5, + pool_max=5, + random_seed=0, + collect_stats=True, + do_final_refine=False, + ) + ens = generate_conformers("CCCCCCCC", config=config) + assert ens.n_conformers > 0 + assert int(ens.generation_stats.get("n_low_mode_seeds", 0)) >= 0 + + +def test_low_mode_following_stat_absent_when_disabled(): + """With use_low_mode_following=False, low_mode_time_s must be 0.""" + config = ConformerConfig( + use_low_mode_following=False, + n_steps=10, + n_seeds=2, + max_out=3, + random_seed=0, + collect_stats=True, + do_final_refine=False, + ) + ens = generate_conformers("CCC", config=config) + stats = ens.generation_stats + assert float(stats.get("low_mode_time_s", 0.0)) == 0.0 + assert int(stats.get("n_low_mode_seeds", 0)) == 0 + + +# --------------------------------------------------------------------------- +# Macrocycle preset +# --------------------------------------------------------------------------- + + +def test_macrocycle_preset_config_values(): + """Macrocycle preset must have wide energy window and low-mode following on.""" + config = preset_config("macrocycle") + assert config.energy_window_kcal == 100.0 + assert config.use_low_mode_following is True + assert config.parent_strategy == "softmax" + assert config.max_out == 200 + assert config.final_select == "diverse" + + +def test_macrocycle_preset_generates_conformers(): + """Macrocycle preset must produce at least one conformer on a simple macrocycle.""" + import dataclasses + + config = dataclasses.replace( + preset_config("macrocycle"), + max_out=5, + n_steps=20, + n_seeds=4, + random_seed=7, + do_final_refine=False, + ) + ens = generate_conformers("C1CCCCCCCCCCC1", config=config) + assert ens.n_conformers > 0 + + diff --git a/tests/test_macrocycles.py b/tests/test_macrocycles.py index 3ebf7c6..240a0c9 100644 --- a/tests/test_macrocycles.py +++ b/tests/test_macrocycles.py @@ -330,3 +330,116 @@ def test_ring_torsion_diversity(name: str, smiles: str, min_sigs: int): ring_atoms = list(ens.mol.GetRingInfo().AtomRings()[0]) sigs = {_torsion_signature(_ring_torsions(ens.mol, r.conf_id, ring_atoms)) for r in ens.records} assert len(sigs) >= min_sigs, f"{name}: found {len(sigs)} torsion families, need ≥ {min_sigs}" + + +# --------------------------------------------------------------------------- +# Unit tests for the amide-flip move +# --------------------------------------------------------------------------- + +# 13-membered macrolactams: N-methyl (tertiary) and N-H (secondary). +_TERTIARY_MACROLACTAM = "O=C1CCCCCCCCCCN1C" +_SECONDARY_MACROLACTAM = "O=C1CCCCCCCCCCN1" + + +def _amide_proposer(smiles: str, seed: int = 1): + from openconf.propose.hybrid import HybridProposer + from openconf.torsionlib import TorsionLibrary + + mol = prepare_molecule(Chem.MolFromSmiles(smiles)) + AllChem.EmbedMolecule(mol, randomSeed=seed) + rm = build_rotor_model(mol) + proposer = HybridProposer(mol, rm, TorsionLibrary(), ConformerConfig(random_seed=0)) + return mol, proposer + + +def test_amide_flip_detected_on_macrolactam(): + """A tertiary in-ring amide in a macrocycle registers exactly one flip site.""" + mol, proposer = _amide_proposer(_TERTIARY_MACROLACTAM) + flips = proposer._moves.amide_flips + assert len(flips) == 1 + c_idx, n_idx, moving = flips[0] + assert mol.GetAtomWithIdx(c_idx).GetAtomicNum() == 6 + assert mol.GetAtomWithIdx(n_idx).GetAtomicNum() == 7 + # Anchors (C, N) must never be in the moving set. + assert c_idx not in moving + assert n_idx not in moving + assert len(moving) > 0 + + +def test_amide_flip_detected_on_secondary_macrolactam(): + """Secondary (N-H) in-ring amides also qualify for the flip move.""" + _, proposer = _amide_proposer(_SECONDARY_MACROLACTAM) + assert len(proposer._moves.amide_flips) == 1 + + +def test_amide_flip_absent_on_cycloalkane(): + """A ring with no amide bond registers no flip sites.""" + _, proposer = _amide_proposer("C1CCCCCCCCCCC1") + assert proposer._moves.amide_flips == [] + + +def test_amide_flip_absent_on_small_lactam(): + """Caprolactam (7-membered) is below the flippable-ring threshold.""" + _, proposer = _amide_proposer("O=C1CCCCCN1C") + assert proposer._moves.amide_flips == [] + + +def test_amide_flip_preserves_ring_bonds(): + """Rotating about the C-N axis keeps every ring bond length fixed (closure-preserving).""" + from openconf.propose.hybrid import _copy_conformer + + mol, proposer = _amide_proposer(_TERTIARY_MACROLACTAM) + assert proposer._moves.amide_flips, "expected a flippable in-ring amide" + ring = mol.GetRingInfo().AtomRings()[0] + orig_id = mol.GetConformers()[0].GetId() + + def bond_lengths(cid: int) -> list[float]: + pos = mol.GetConformer(cid).GetPositions() + return [float(np.linalg.norm(pos[ring[i]] - pos[ring[(i + 1) % len(ring)]])) for i in range(len(ring))] + + before = bond_lengths(orig_id) + for _ in range(20): + new_id = _copy_conformer(mol, orig_id) + proposer._moves.apply_amide_flip_move(new_id) + after = bond_lengths(new_id) + for b, a in zip(before, after, strict=True): + assert abs(a - b) < 1e-6, f"Ring bond length drifted under amide flip: {b} -> {a}" + mol.RemoveConformer(new_id) + + +def test_amide_flip_inverts_amide_dihedral(): + """The O=C-N-C dihedral changes by roughly 180 degrees after a flip.""" + from openconf.propose.amide_seeds import find_ring_tertiary_amide_dihedrals + from openconf.propose.hybrid import _copy_conformer + + mol, proposer = _amide_proposer(_TERTIARY_MACROLACTAM) + quad = find_ring_tertiary_amide_dihedrals(mol)[0] + orig_id = mol.GetConformers()[0].GetId() + + # Average |delta| over several flips should sit near 180 (the move targets 180 + jitter). + deltas = [] + for _ in range(15): + new_id = _copy_conformer(mol, orig_id) + before = rdMolTransforms.GetDihedralDeg(mol.GetConformer(new_id), *quad) + proposer._moves.apply_amide_flip_move(new_id) + after = rdMolTransforms.GetDihedralDeg(mol.GetConformer(new_id), *quad) + diff = abs((after - before + 180.0) % 360.0 - 180.0) + deltas.append(diff) + mol.RemoveConformer(new_id) + + mean_delta = float(np.mean(deltas)) + assert mean_delta > 120.0, f"Amide flips did not substantially invert the dihedral: mean |delta|={mean_delta:.1f}" + + +def test_amide_flip_generation_runs(): + """End-to-end generation with an amide-flip-heavy schedule produces conformers.""" + config = ConformerConfig( + max_out=8, + n_steps=40, + n_seeds=6, + random_seed=7, + do_final_refine=False, + move_probs={"amide_flip": 0.6, "single_rotor": 0.2, "crankshaft": 0.2}, + ) + ens = generate_conformers(_TERTIARY_MACROLACTAM, config=config) + assert ens.n_conformers > 0 diff --git a/uv.lock b/uv.lock index 563b701..2147439 100644 --- a/uv.lock +++ b/uv.lock @@ -260,7 +260,7 @@ wheels = [ [[package]] name = "openconf" -version = "0.0.14" +version = "0.0.15" source = { editable = "." } dependencies = [ { name = "numpy" }, From 3d279303bb2b5a354a4aae2b05e578058e6ba364 Mon Sep 17 00:00:00 2001 From: ncasetti7 Date: Mon, 15 Jun 2026 14:22:52 -0400 Subject: [PATCH 2/4] Fix ty type errors in find_ring_amide_bonds --- openconf/propose/amide_seeds.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/openconf/propose/amide_seeds.py b/openconf/propose/amide_seeds.py index 161b22d..0cfd056 100644 --- a/openconf/propose/amide_seeds.py +++ b/openconf/propose/amide_seeds.py @@ -42,14 +42,16 @@ def find_ring_amide_bonds( ring_atom_indices is the ordered atom ring containing the bond """ ring_info = mol.GetRingInfo() - rings = sorted((r for r in ring_info.AtomRings() if len(r) >= min_ring_size), key=len, reverse=True) + # Pre-build (ring_set, ring_atoms) pairs in a list comprehension so ty can + # resolve the element type from AtomRings() without going through sorted(). + eligible = [(set(r), tuple(r)) for r in ring_info.AtomRings() if len(r) >= min_ring_size] + rings = sorted(eligible, key=lambda x: len(x[1]), reverse=True) results: list[tuple[int, int, tuple[int, ...]]] = [] seen: set[tuple[int, int]] = set() - for ring in rings: - ring_set = set(ring) - for n_idx in ring: + for ring_set, ring_atoms in rings: + for n_idx in ring_set: n_atom = mol.GetAtomWithIdx(n_idx) if n_atom.GetAtomicNum() != 7 or n_atom.GetIsAromatic(): continue @@ -75,7 +77,7 @@ def find_ring_amide_bonds( continue seen.add(bond_key) - results.append((c_idx, n_idx, tuple(ring))) + results.append((c_idx, n_idx, ring_atoms)) return results From 49b938c4e0dabc07a8c4b96b9a0fe15be39385e4 Mon Sep 17 00:00:00 2001 From: ncasetti7 Date: Tue, 16 Jun 2026 10:15:53 -0400 Subject: [PATCH 3/4] Replace gap heuristic with T/R projection; normalize move probs --- openconf/data/runtime_tuning.json | 10 ++-- openconf/propose/low_mode.py | 97 ++++++++++++++++++++++++------- tests/test_low_mode.py | 17 +++--- 3 files changed, 88 insertions(+), 36 deletions(-) diff --git a/openconf/data/runtime_tuning.json b/openconf/data/runtime_tuning.json index b56c0a2..82d7b3c 100644 --- a/openconf/data/runtime_tuning.json +++ b/openconf/data/runtime_tuning.json @@ -44,11 +44,11 @@ "single_rotor": 0.3, "multi_rotor": 0.24, "correlated": 0.16, - "global_shake": 0.08, - "ring_flip": 0.1, - "crankshaft": 0.06, - "ring_kic": 0.06, - "amide_flip": 0.06 + "global_shake": 0.06, + "ring_flip": 0.09, + "crankshaft": 0.05, + "ring_kic": 0.05, + "amide_flip": 0.05 }, "forced_periodic_move": "global_shake", "constrained_fallback_move": "single_rotor", diff --git a/openconf/propose/low_mode.py b/openconf/propose/low_mode.py index ff172a2..a6bd078 100644 --- a/openconf/propose/low_mode.py +++ b/openconf/propose/low_mode.py @@ -15,8 +15,6 @@ if TYPE_CHECKING: from ..relax import Minimizer -_N_RIGID_MODES: int = 6 -_RIGID_BODY_GAP_WINDOW: int = _N_RIGID_MODES + 2 _DEFAULT_FD_STEP: float = 0.005 _DEFAULT_EIGENVALUE_THRESHOLD: float = 100.0 _DEFAULT_MAX_MODES: int = 5 @@ -83,26 +81,81 @@ def _compute_hessian( return (hessian + hessian.T) * 0.5 +def _null_space(a: np.ndarray) -> np.ndarray: + """Orthonormal basis for the null space of a, via SVD.""" + _, s, vt = np.linalg.svd(a, full_matrices=True) + rcond = max(a.shape) * np.finfo(float).eps * (s[0] if len(s) > 0 else 1.0) + rank = int(np.sum(s > rcond)) + return vt[rank:].T + + +def _build_vibrational_basis(mol: Chem.Mol, conf_id: int) -> np.ndarray: + """Orthonormal basis for the vibrational subspace, shape (3N, 3N−6). + + Builds the six rigid-body vectors (3 translations + 3 rotations) in + Cartesian coordinates, then returns their null space — the vibrational + subspace. Near-zero rotation vectors (linear molecules have one) are + dropped before the SVD so the basis always spans exactly 3N−6 or 3N−5 + dimensions. + + Args: + mol: molecule supplying atom positions and masses + conf_id: conformer from which positions are taken + + Returns: + Array of shape (3N, 3N−6) or (3N, 3N−5) for linear molecules + """ + n_atoms = mol.GetNumAtoms() + n_dof = 3 * n_atoms + pt = Chem.GetPeriodicTable() + masses = np.array([pt.GetAtomicWeight(a.GetAtomicNum()) for a in mol.GetAtoms()]) + pos = mol.GetConformer(conf_id).GetPositions() # (N, 3) + + com = np.sum(masses[:, None] * pos, axis=0) / np.sum(masses) + r = pos - com # COM-centred positions (N, 3) + + # Translation: uniform displacement along each Cartesian axis + trans = np.zeros((3, n_dof)) + for axis in range(3): + trans[axis, axis::3] = 1.0 + trans /= np.linalg.norm(trans, axis=1, keepdims=True) + + # Rotation: (r × e_axis) per atom, flattened into a 3N vector + rot = np.zeros((3, n_dof)) + for k in range(n_atoms): + rx, ry, rz = r[k] + rot[0, 3 * k : 3 * k + 3] = [0.0, rz, -ry] # r × e_x + rot[1, 3 * k : 3 * k + 3] = [-rz, 0.0, rx] # r × e_y + rot[2, 3 * k : 3 * k + 3] = [ ry, -rx, 0.0] # r × e_z + + tr_vecs = list(trans) + for d in rot: + norm = float(np.linalg.norm(d)) + if norm > 1e-8 * np.sqrt(float(n_atoms)): + tr_vecs.append(d / norm) + + return _null_space(np.stack(tr_vecs)) + + def _select_low_modes( hessian: np.ndarray, + mol: Chem.Mol, + conf_id: int, eigenvalue_threshold: float, max_modes: int, ) -> np.ndarray: - """Select low-frequency eigenvectors of a Hessian. - - Diagonalizes the Hessian, locates the rigid-body / conformational boundary - via the largest eigenvalue gap within the first _RIGID_BODY_GAP_WINDOW - entries, skips all modes up to that boundary, and returns eigenvectors - whose eigenvalues are below eigenvalue_threshold, capped at max_modes. + """Select low-frequency eigenvectors by projecting out rigid-body modes. - Using a gap rather than a fixed count of 6 handles linear molecules - (only 5 rigid-body modes) and numerical drift when minimization is not - fully converged. + Builds the Cartesian translation/rotation basis for the conformer, computes + its null space (the vibrational subspace), projects the Hessian into that + subspace, and diagonalises. Eigenvectors whose eigenvalues are below + eigenvalue_threshold are returned as Cartesian unit-norm displacement vectors. Args: - hessian: symmetric 3N×3N Hessian matrix - eigenvalue_threshold: upper eigenvalue bound (kcal/mol/Ų) for mode - selection; modes at or above this value are excluded + hessian: symmetric 3N×3N Hessian matrix in kcal/mol/Ų + mol: molecule providing atom positions and masses + conf_id: conformer ID used to build the translation/rotation basis + eigenvalue_threshold: upper eigenvalue bound (kcal/mol/Ų) max_modes: maximum number of modes to return Returns: @@ -110,16 +163,16 @@ def _select_low_modes( eigenvectors in ascending eigenvalue order; shape (3N, 0) when no conformational modes satisfy the threshold """ - eigenvalues, eigenvectors = np.linalg.eigh(hessian) - window = min(_RIGID_BODY_GAP_WINDOW, len(eigenvalues) - 1) - n_skip = int(np.argmax(np.diff(eigenvalues[: window + 1]))) + 1 - conf_vals = eigenvalues[n_skip:] - conf_vecs = eigenvectors[:, n_skip:] + d_vib = _build_vibrational_basis(mol, conf_id) + h_vib = d_vib.T @ hessian @ d_vib + eigenvalues, eigenvectors = np.linalg.eigh(h_vib) + # d_vib has orthonormal columns so d_vib @ eigenvectors is orthonormal in 3N space + vecs = d_vib @ eigenvectors - mask = conf_vals < eigenvalue_threshold + mask = eigenvalues < eigenvalue_threshold if not np.any(mask): return np.empty((hessian.shape[0], 0)) - return conf_vecs[:, mask][:, :max_modes] + return vecs[:, mask][:, :max_modes] def _scan_along_mode( @@ -255,7 +308,7 @@ def generate_low_mode_seeds( minimizations fail """ hessian = _compute_hessian(mol, ff_props, conf_id, fd_step) - low_vecs = _select_low_modes(hessian, eigenvalue_threshold, max_modes) + low_vecs = _select_low_modes(hessian, mol, conf_id, eigenvalue_threshold, max_modes) if low_vecs.shape[1] == 0: return [] diff --git a/tests/test_low_mode.py b/tests/test_low_mode.py index 7aada98..8669f42 100644 --- a/tests/test_low_mode.py +++ b/tests/test_low_mode.py @@ -13,7 +13,6 @@ _DEFAULT_SCAN_ENERGY_THRESHOLD, _DEFAULT_SCAN_MAX_STEPS, _DEFAULT_SCAN_STEP_SIZE, - _N_RIGID_MODES, _compute_hessian, _scan_along_mode, _select_low_modes, @@ -75,7 +74,7 @@ def test_hessian_positive_eigenvalues_at_minimum(): mol, props = _make_minimized_mol("CCC") H = _compute_hessian(mol, props, conf_id=0) eigenvalues = np.linalg.eigvalsh(H) - conformational = eigenvalues[_N_RIGID_MODES:] + conformational = eigenvalues[6:] neg = conformational[conformational < -0.5] assert np.all(conformational > -0.5), f"Unexpected large negative eigenvalues: {neg}" @@ -90,7 +89,7 @@ def test_select_low_modes_returns_correct_shape(): mol, props = _make_minimized_mol("CCC") H = _compute_hessian(mol, props, conf_id=0) n_dof = 3 * mol.GetNumAtoms() - vecs = _select_low_modes(H, _DEFAULT_EIGENVALUE_THRESHOLD, _DEFAULT_MAX_MODES) + vecs = _select_low_modes(H, mol, 0, _DEFAULT_EIGENVALUE_THRESHOLD, _DEFAULT_MAX_MODES) assert vecs.shape[0] == n_dof assert vecs.shape[1] <= _DEFAULT_MAX_MODES @@ -99,7 +98,7 @@ def test_select_low_modes_empty_when_threshold_zero(): """No modes returned when threshold is effectively zero.""" mol, props = _make_minimized_mol("CCC") H = _compute_hessian(mol, props, conf_id=0) - vecs = _select_low_modes(H, eigenvalue_threshold=0.0, max_modes=10) + vecs = _select_low_modes(H, mol, 0, eigenvalue_threshold=0.0, max_modes=10) assert vecs.shape[1] == 0 @@ -108,7 +107,7 @@ def test_select_low_modes_respects_max_modes(): mol, props = _make_minimized_mol("CCCCC") H = _compute_hessian(mol, props, conf_id=0) for cap in [1, 2, 3]: - vecs = _select_low_modes(H, eigenvalue_threshold=500.0, max_modes=cap) + vecs = _select_low_modes(H, mol, 0, eigenvalue_threshold=500.0, max_modes=cap) assert vecs.shape[1] <= cap, f"Got {vecs.shape[1]} modes with cap={cap}" @@ -116,7 +115,7 @@ def test_select_low_modes_columns_are_unit_vectors(): """Returned eigenvectors must have unit norm (they come from eigh).""" mol, props = _make_minimized_mol("CCC") H = _compute_hessian(mol, props, conf_id=0) - vecs = _select_low_modes(H, _DEFAULT_EIGENVALUE_THRESHOLD, _DEFAULT_MAX_MODES) + vecs = _select_low_modes(H, mol, 0, _DEFAULT_EIGENVALUE_THRESHOLD, _DEFAULT_MAX_MODES) if vecs.shape[1] == 0: pytest.skip("No low modes found — threshold may need adjustment") norms = np.linalg.norm(vecs, axis=0) @@ -132,7 +131,7 @@ def test_scan_along_mode_moves_from_start(): """Scan must return positions that differ from the starting geometry.""" mol, props = _make_minimized_mol("CCC") H = _compute_hessian(mol, props, conf_id=0) - vecs = _select_low_modes(H, _DEFAULT_EIGENVALUE_THRESHOLD, 1) + vecs = _select_low_modes(H, mol, 0, _DEFAULT_EIGENVALUE_THRESHOLD, 1) if vecs.shape[1] == 0: pytest.skip("No low modes found") @@ -151,7 +150,7 @@ def test_scan_along_mode_restores_start_conformer(): """The start conformer must be unchanged after scanning.""" mol, props = _make_minimized_mol("CCC") H = _compute_hessian(mol, props, conf_id=0) - vecs = _select_low_modes(H, _DEFAULT_EIGENVALUE_THRESHOLD, 1) + vecs = _select_low_modes(H, mol, 0, _DEFAULT_EIGENVALUE_THRESHOLD, 1) if vecs.shape[1] == 0: pytest.skip("No low modes found") @@ -173,7 +172,7 @@ def test_scan_stops_at_energy_threshold(): """With a near-zero energy threshold, scan must stop after the first step.""" mol, props = _make_minimized_mol("CCC") H = _compute_hessian(mol, props, conf_id=0) - vecs = _select_low_modes(H, _DEFAULT_EIGENVALUE_THRESHOLD, 1) + vecs = _select_low_modes(H, mol, 0, _DEFAULT_EIGENVALUE_THRESHOLD, 1) if vecs.shape[1] == 0: pytest.skip("No low modes found") From 4c9218ac47cec6bc495400dd2256cec782e93894 Mon Sep 17 00:00:00 2001 From: ncasetti7 Date: Tue, 16 Jun 2026 10:22:42 -0400 Subject: [PATCH 4/4] formatting fix --- openconf/propose/low_mode.py | 8 +++---- tests/test_low_mode.py | 41 +++++++++++++++++++++--------------- 2 files changed, 28 insertions(+), 21 deletions(-) diff --git a/openconf/propose/low_mode.py b/openconf/propose/low_mode.py index a6bd078..be81817 100644 --- a/openconf/propose/low_mode.py +++ b/openconf/propose/low_mode.py @@ -120,13 +120,13 @@ def _build_vibrational_basis(mol: Chem.Mol, conf_id: int) -> np.ndarray: trans[axis, axis::3] = 1.0 trans /= np.linalg.norm(trans, axis=1, keepdims=True) - # Rotation: (r × e_axis) per atom, flattened into a 3N vector + # Rotation: (r cross e_axis) per atom, flattened into a 3N vector rot = np.zeros((3, n_dof)) for k in range(n_atoms): rx, ry, rz = r[k] - rot[0, 3 * k : 3 * k + 3] = [0.0, rz, -ry] # r × e_x - rot[1, 3 * k : 3 * k + 3] = [-rz, 0.0, rx] # r × e_y - rot[2, 3 * k : 3 * k + 3] = [ ry, -rx, 0.0] # r × e_z + rot[0, 3 * k : 3 * k + 3] = [0.0, rz, -ry] # r cross e_x + rot[1, 3 * k : 3 * k + 3] = [-rz, 0.0, rx] # r cross e_y + rot[2, 3 * k : 3 * k + 3] = [ry, -rx, 0.0] # r cross e_z tr_vecs = list(trans) for d in rot: diff --git a/tests/test_low_mode.py b/tests/test_low_mode.py index 8669f42..38060b3 100644 --- a/tests/test_low_mode.py +++ b/tests/test_low_mode.py @@ -140,8 +140,13 @@ def test_scan_along_mode_moves_from_start(): start_pos = mol.GetConformer(0).GetPositions().copy() final_pos = _scan_along_mode( - mol, props, 0, direction, _DEFAULT_SCAN_STEP_SIZE, - _DEFAULT_SCAN_ENERGY_THRESHOLD, _DEFAULT_SCAN_MAX_STEPS, + mol, + props, + 0, + direction, + _DEFAULT_SCAN_STEP_SIZE, + _DEFAULT_SCAN_ENERGY_THRESHOLD, + _DEFAULT_SCAN_MAX_STEPS, ) assert not np.allclose(final_pos, start_pos), "Scan returned starting positions" @@ -160,8 +165,13 @@ def test_scan_along_mode_restores_start_conformer(): n_confs_before = mol.GetNumConformers() _scan_along_mode( - mol, props, 0, direction, _DEFAULT_SCAN_STEP_SIZE, - _DEFAULT_SCAN_ENERGY_THRESHOLD, _DEFAULT_SCAN_MAX_STEPS, + mol, + props, + 0, + direction, + _DEFAULT_SCAN_STEP_SIZE, + _DEFAULT_SCAN_ENERGY_THRESHOLD, + _DEFAULT_SCAN_MAX_STEPS, ) assert mol.GetNumConformers() == n_confs_before, "Scan left a temporary conformer behind" @@ -183,8 +193,13 @@ def test_scan_stops_at_energy_threshold(): # Threshold of 0 means any energy increase (even numerical noise) stops the scan. # The returned positions should still be the start positions (first step rejected). final_pos = _scan_along_mode( - mol, props, 0, direction, _DEFAULT_SCAN_STEP_SIZE, - energy_threshold=0.0, max_steps=10, + mol, + props, + 0, + direction, + _DEFAULT_SCAN_STEP_SIZE, + energy_threshold=0.0, + max_steps=10, ) # Either no progress or exactly one step taken — at most a tiny displacement displacement = float(np.linalg.norm(final_pos - start_pos)) @@ -233,9 +248,7 @@ def test_generate_low_mode_seeds_empty_when_threshold_zero( """No seeds generated when eigenvalue threshold excludes all conformational modes.""" mol, props = propane_mol_props minimizer = _make_fast_minimizer(mol) - seeds = generate_low_mode_seeds( - mol, props, conf_id=0, minimizer=minimizer, eigenvalue_threshold=0.0 - ) + seeds = generate_low_mode_seeds(mol, props, conf_id=0, minimizer=minimizer, eigenvalue_threshold=0.0) assert seeds == [] @@ -247,9 +260,7 @@ def test_generate_low_mode_seeds_at_most_two_per_mode( minimizer = _make_fast_minimizer(mol) for cap in [1, 2]: n_before = mol.GetNumConformers() - seeds = generate_low_mode_seeds( - mol, props, conf_id=0, minimizer=minimizer, max_modes=cap - ) + seeds = generate_low_mode_seeds(mol, props, conf_id=0, minimizer=minimizer, max_modes=cap) assert len(seeds) <= 2 * cap, f"Got {len(seeds)} seeds with cap={cap}" for conf_id, _ in seeds: mol.RemoveConformer(conf_id) @@ -261,9 +272,7 @@ def test_generate_low_mode_seeds_scans_both_directions() -> None: mol, props = _make_minimized_mol("CCCCCC") minimizer = _make_fast_minimizer(mol) - seeds = generate_low_mode_seeds( - mol, props, conf_id=0, minimizer=minimizer, max_modes=1 - ) + seeds = generate_low_mode_seeds(mol, props, conf_id=0, minimizer=minimizer, max_modes=1) # For one mode with two directions we expect 0 or 2 seeds (not 1), # because both directions are always attempted. assert len(seeds) in (0, 2), f"Expected 0 or 2 seeds for 1 mode, got {len(seeds)}" @@ -452,5 +461,3 @@ def test_macrocycle_preset_generates_conformers(): ) ens = generate_conformers("C1CCCCCCCCCCC1", config=config) assert ens.n_conformers > 0 - -