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

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 26 additions & 20 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -279,7 +282,7 @@ ensemble.to_sdf("output.sdf")
<details><summary>Full config equivalent</summary>

```python
from openconf import ConformerConfig, PrismConfig, generate_conformers
from openconf import ConformerConfig, generate_conformers

config = ConformerConfig(
max_out=250,
Expand All @@ -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")
Expand Down Expand Up @@ -354,7 +356,7 @@ ensemble = generate_conformers_from_pose(mol, constrained_atoms=ring_atoms, conf
<details><summary>Full analogue preset equivalent</summary>

```python
from openconf import ConformerConfig, PrismConfig, generate_conformers_from_pose
from openconf import ConformerConfig, generate_conformers_from_pose

config = ConformerConfig(
max_out=50,
Expand All @@ -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)
```
Expand All @@ -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
Expand All @@ -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:
Expand All @@ -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
Expand All @@ -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
Expand Down
43 changes: 26 additions & 17 deletions SCIENCE.md
Original file line number Diff line number Diff line change
Expand Up @@ -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`.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down
Loading
Loading