bio-conformer-generation

What this skill does

## Version Compatibility

Reference examples tested with: RDKit 2024.09+, xtb 6.7+, CREST 3.0+, OpenMM 8.1+ for follow-up MD.

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show <package>` then `help(module.function)` to check signatures
- CLI: `xtb --version`; `crest --version`

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# Conformer Generation

Generate 3D conformer ensembles for molecules from 2D structures. The choice of method depends on molecule size, flexibility, and downstream use: ETKDG / ETKDGv3 (Riniker & Landrum 2015 J Chem Inf Model 55:2562-2574; knowledge-enhanced distance geometry) is the modern default for drug-like molecules, MMFF94/UFF for fast energy minimization, CREST + GFN2-xTB for high-accuracy semi-empirical sampling of macrocycles and peptides. A single conformer is rarely sufficient: descriptor variance across the ensemble can exceed the descriptor signal, and docking pose accuracy degrades if the starting conformer is non-bioactive.

For docking pose validation, see `chemoinformatics/pose-validation`. For free-energy methods (which require ensemble sampling), see `chemoinformatics/free-energy-calculations`.

## Conformer Method Taxonomy

| Method | Cost / mol | Quality | Use case | Fails when |
|--------|-----------|---------|----------|------------|
| ETKDGv3 + MMFF94 | <1s | Good for drug-like | Default; docking input; descriptors | Macrocycles, peptides, transition metals |
| ETKDGv3 + UFF | <1s | Lower-quality MMFF94 alternative | Fallback when MMFF94 fails to parameterize | Same as MMFF94 |
| Omega (OpenEye) | 1s | Industry-standard commercial | Commercial pipelines | License cost |
| Confab (Open Babel) | 5s | Systematic torsion search | Patent expiration | Quality limited |
| RDKit ETKDGv3 + macrocycle preferences | 10-60s | Drug-like macrocycles | Macrocyclic peptides | Still limited; CREST better |
| CREST + GFN2-xTB | minutes | High-accuracy semi-empirical | Macrocycles, peptides, conformer ensembles for QSAR | Computationally expensive; metal centers |
| CREST + GFN-FF | seconds | GFN2 quality at FF speed | Quick screening | Limited element coverage |
| GeoMol (Ganea 2021) | <0.1s GPU | ML-fast, ETKDGv3-quality | Large library 3D conformers | ML training distribution |
| TorsionNet (Gogineni 2020) | <0.1s GPU | ML-fast | Drug-like | ML training distribution |
| MD sampling (OpenMM) | hours | High-quality dynamic | Free energy, induced fit | Computational cost |

**Decision:** For drug-like molecules (<500 Da, <8 rotatable bonds), **ETKDGv3 + MMFF94** with 20-100 conformers is the modern default. For macrocycles, peptides, or molecules with >12 rotatable bonds, **CREST + GFN2-xTB** captures the conformational diversity. For ML-scale (>1M molecules), **GeoMol** trades accuracy for speed.

## Decision Tree by Scenario

| Scenario | Method | Conformer count | Energy window |
|----------|--------|-----------------|----------------|
| Single docking pose (initial 3D) | ETKDGv3 + MMFF94 | 1 | n/a |
| Multi-conformer docking | ETKDGv3 + MMFF94 | 10-50 | 10 kcal/mol |
| 3D QSAR descriptor input | ETKDGv3 + MMFF94 | 50-200 | 5 kcal/mol |
| Pharmacophore search | ETKDGv3 + MMFF94 | 100-500 | 5 kcal/mol |
| Macrocycle / peptide | CREST + GFN2-xTB | 50-200 (auto from CREST) | 5-8 kcal/mol |
| FEP input | CREST + GFN2-xTB then MD relax | 1-3 representative | 3 kcal/mol |
| Bioactive conformer search | ETKDGv3 + MMFF94 then dock with rescore | 100-500 | 10 kcal/mol |
| Shape similarity / ROCS | ETKDGv3 + MMFF94 | 50-200 | 10 kcal/mol |
| Conformer-dependent descriptors | ETKDGv3 ensemble + Boltzmann avg | 20-100 | 5 kcal/mol |

## ETKDGv3 (Modern Default)

ETKDGv3 (Riniker & Landrum 2015) incorporates experimental torsion preferences into distance geometry: starts from random embeddings, refines by satisfying experimentally-derived bond, angle, and torsion preferences.

**Goal:** Generate an ensemble of 3D conformers from a SMILES with the modern default embedding algorithm.

**Approach:** Add explicit hydrogens, configure ETKDGv3 params (random seed, max attempts, random coords), and embed multiple conformers via `EmbedMultipleConfs`.

```python
from rdkit import Chem
from rdkit.Chem import AllChem

def gen_conformers(smiles, n_conf=20, seed=42):
    mol = Chem.MolFromSmiles(smiles)
    mol = Chem.AddHs(mol)
    params = AllChem.ETKDGv3()
    params.randomSeed = seed
    params.useRandomCoords = True
    params.maxAttempts = 1000
    ids = AllChem.EmbedMultipleConfs(mol, numConfs=n_conf, params=params)
    return mol, list(ids)
```

`useRandomCoords=True` improves convergence for macrocycles and heavily-rotated molecules. `maxAttempts=1000` handles difficult embeddings.

## Force-Field Optimization

After embedding, minimize each conformer to a local minimum.

**Goal:** Reduce strain in each embedded conformer to a stable local minimum and record the resulting energies.

**Approach:** Build MMFF94s force-field parameters, minimize each conformer in place, and collect energies; fall back to UFF when MMFF94 cannot parameterize the molecule.

```python
def optimize_conformers(mol, conf_ids, force_field='mmff94'):
    energies = []
    if force_field == 'mmff94':
        mmff_props = AllChem.MMFFGetMoleculeProperties(mol, mmffVariant='MMFF94s')
        for cid in conf_ids:
            ff = AllChem.MMFFGetMoleculeForceField(mol, mmff_props, confId=cid)
            ff.Minimize()
            energies.append(ff.CalcEnergy())
    else:  # UFF fallback
        for cid in conf_ids:
            ff = AllChem.UFFGetMoleculeForceField(mol, confId=cid)
            ff.Minimize()
            energies.append(ff.CalcEnergy())
    return energies
```

**MMFF94 vs MMFF94s:** MMFF94s is the "standard" set with simpler aromatic nitrogen handling; preferred for most drug-like.

**UFF (Universal Force Field):** Lower quality but handles any element including transition metals. Use as fallback when MMFF94 cannot parameterize (uncommon elements, charged species).

## RMSD Pruning

Remove near-duplicate conformers within a chosen RMSD cutoff to keep the ensemble diverse:

```python
import numpy as np

def prune_conformers_rmsd(mol, conf_ids, rmsd_cutoff=0.5):
    n = len(conf_ids)
    keep = []
    for i, cid in enumerate(conf_ids):
        is_unique = True
        for kept_cid in keep:
            rmsd = AllChem.GetBestRMS(mol, mol, cid, kept_cid)
            if rmsd < rmsd_cutoff:
                is_unique = False
                break
        if is_unique:
            keep.append(cid)
    return keep
```

**Typical RMSD cutoff (Source / Rationale):**

| Cutoff | Use case | Source |
|--------|----------|--------|
| 0.5 Å | Drug-like ensemble for descriptors / docking | Empirical: below this conformers represent same minimum (Hawkins 2007) |
| 1.0 Å | Drug-like ensemble for pharmacophore | Standard ROCS / pharmacophore practice |
| 1.5-2.0 Å | Macrocycles / peptides | Higher conformational freedom; Tan 2018 macrocycle benchmarks |
| 2.0+ Å | Cluster-centroid representative ensembles | Coarse representative sampling |

## Energy Window Filtering

Remove conformers above an energy cutoff (high-energy conformers are unlikely to be bioactive):

```python
def filter_by_energy(mol, conf_ids, energies, window_kcal=10.0):
    min_e = min(energies)
    keep = []
    for cid, e in zip(conf_ids, energies):
        if e - min_e <= window_kcal:
            keep.append(cid)
    return keep
```

**Window choice:**
- 3 kcal/mol: very strict, only near-global-min conformers (FEP, MD setup)
- 5 kcal/mol: typical for 3D QSAR, pharmacophore
- 10 kcal/mol: typical for docking input (bioactive conformer may be higher)
- 25 kcal/mol: macrocycles, no filter (bioactive conformer can be high-energy when bound)

## Macrocycle Handling

Macrocycles (>=12 atom rings) have distinct conformational issues: ETKDGv3 defau