{
  "skill": "alterlab-datamol",
  "evals": [
    {
      "id": "standardize-and-descriptors",
      "prompt": "I have a CSV of ~3000 vendor compounds with a SMILES column. Standardize them (disconnect metals, normalize, reionize), drop anything that fails to parse, and compute MW, logP, HBD, HBA and TPSA for the survivors so I can filter by Lipinski.",
      "expected_output": "Invokes alterlab-datamol: reads the CSV with dm.read_csv, standardizes each molecule via dm.standardize_mol (disconnect_metals/normalize/reionize) while dropping None parses from dm.to_mol, then runs dm.descriptors.batch_compute_many_descriptors with n_jobs=-1 (and an explicit batch_size, which joblib >= 1.6 requires) to get datamol's keys mw/clogp/n_lipinski_hbd/n_lipinski_hba/tpsa, returning a clean DataFrame ready for a Rule-of-Five mask.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "behavior", "value": "Uses datamol (dm.*) for parsing, standardization, and dm.descriptors.batch_compute_many_descriptors rather than raw RDKit, handles failed SMILES parses, and references the descriptors by datamol's own key names (clogp, n_lipinski_hbd, n_lipinski_hba) — not RDKit's logp/hbd/hba." }
      ]
    },
    {
      "id": "butina-clustering-diversity",
      "prompt": "Cluster these 800 hits by structural similarity at a Tanimoto distance cutoff of 0.2, then pick the 50 most diverse representatives for follow-up.",
      "expected_output": "Invokes alterlab-datamol: runs dm.cluster_mols with cutoff=0.2 and n_jobs=-1 (Butina clustering, noting it builds a full distance matrix suitable at this ~800 scale), unpacks its (cluster_indices, cluster_mols) return value to report cluster membership, and uses dm.pick_diverse / dm.pick_centroids (each returning (indices, mols)) to select 50 representatives.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "behavior", "value": "Uses dm.cluster_mols (Butina) at the 0.2 cutoff and a datamol diversity picker, and notes the full-distance-matrix scale limit." }
      ]
    },
    {
      "id": "scaffold-grouping-sar",
      "prompt": "Extract the Bemis-Murcko scaffold for each molecule in my series, group the compounds by scaffold, and show me the 10 most common scaffolds.",
      "expected_output": "Invokes alterlab-datamol: computes Bemis-Murcko scaffolds with dm.to_scaffold_murcko, canonicalizes each via dm.to_smiles, groups molecules per scaffold and uses a Counter to surface the 10 most frequent scaffold SMILES.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "dm.to_scaffold_murcko" }
      ]
    },
    {
      "id": "conformers-and-viz",
      "prompt": "Generate 3D conformers for this molecule with ETKDGv3, minimize them, filter near-duplicates by RMSD, and save an aligned grid image of the molecules to a PNG.",
      "expected_output": "Invokes alterlab-datamol: calls dm.conformers.generate (method='ETKDGv3', minimize_energy=True, rms_cutoff to drop duplicates), then renders with dm.viz.to_image (align=True, outfile PNG with use_svg=False, since SVG is the default). 3D conformer embedding plus visualization are datamol's wheelhouse.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "behavior", "value": "Uses dm.conformers.generate with ETKDGv3 and dm.viz.to_image to write the PNG grid." }
      ]
    },
    {
      "id": "near-miss-molfeat",
      "prompt": "I want to turn these SMILES into ML-ready feature matrices and benchmark ECFP vs MACCS vs ChemBERTa embeddings for a QSAR random forest. Build the featurizers and a scikit-learn pipeline.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-molfeat. The user wants ML featurizer objects, pretrained ChemBERTa embeddings, and a scikit-learn-compatible transformer pipeline to compare representations, which is molfeat's territory, not datamol's everyday descriptor/cleanup wrapper.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-molfeat" }
      ]
    },
    {
      "id": "near-miss-rdkit",
      "prompt": "I need fine-grained control over a custom RDKit reaction with explicit atom-map sanitization flags and a hand-tuned conformer embedding using my own ETKDG parameters and bounds matrix. Show me the raw RDKit calls.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-rdkit. The user explicitly wants low-level, custom-parameter RDKit control rather than datamol's sensible-default wrapper, which the datamol skill itself recommends deferring to for advanced control.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-rdkit" }
      ]
    },
    {
      "id": "near-miss-medchem",
      "prompt": "Before I order these 2,000 screening hits, flag PAINS and other structural alerts, apply the NIBR filters, and drop anything failing Rule of Five or Veber.",
      "expected_output": "Does NOT invoke alterlab-datamol; defers to alterlab-medchem. Rule-based drug-likeness triage with PAINS / structural-alert catalogs and NIBR filters is medchem's job (RuleFilters, CommonAlertsFilters, NIBRFilters); datamol only supplies the underlying molecule handling and descriptors.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-medchem" }
      ]
    }
  ]
}

# Datamol Conformers Module Reference

The `datamol.conformers` module provides tools for generating and analyzing 3D molecular conformations.

## Conformer Generation

### `dm.conformers.generate(mol, n_confs=None, rms_cutoff=None, minimize_energy=False, method=None, add_hs=True, ...)`
Generate 3D molecular conformers.
- **Parameters**:
  - `mol`: Input molecule
  - `n_confs`: Number of conformers to generate (auto-determined based on rotatable bonds if None)
  - `rms_cutoff`: RMS threshold in Ångströms for filtering similar conformers (removes duplicates)
  - `minimize_energy`: Apply UFF energy minimization (default: **False** — pass `True` to minimize)
  - `method`: Embedding method - options:
    - `'ETDG'` - Experimental Torsion Distance Geometry
    - `'ETKDG'` - ETDG with additional basic knowledge
    - `'ETKDGv2'` - Enhanced version 2
    - `'ETKDGv3'` - Enhanced version 3 (used when `method=None`; recommended)
  - `add_hs`: Add hydrogens before embedding (default: True, critical for quality)
  - `random_seed`: Set for reproducibility
- **Returns**: Molecule with embedded conformers
- **Example**:
  ```python
  mol = dm.to_mol("CCO")
  mol_3d = dm.conformers.generate(mol, n_confs=10, rms_cutoff=0.5)
  conformers = mol_3d.GetConformers()  # Access all conformers
  ```

## Conformer Clustering

### `dm.conformers.cluster(mol, rms_cutoff=1.0, already_aligned=False, centroids=True)`
Group conformers by RMS distance.
- **Parameters**:
  - `rms_cutoff`: Clustering threshold in Ångströms (default: 1.0)
  - `already_aligned`: Whether conformers are pre-aligned
  - `centroids`: Return centroid conformers (True, default) or cluster groups (False)
- **Returns**: Cluster information or centroid conformers
- **Use case**: Identify distinct conformational families

### `dm.conformers.return_centroids(mol, conf_clusters, centroids=True)`
Extract representative conformers from clusters.
- **Parameters**:
  - `conf_clusters`: Sequence of index clusters (e.g. from RDKit `Butina.ClusterData`); `cluster()` calls this internally and already returns its result
  - `centroids`: Return single molecule (True) or list of molecules (False)
- **Returns**: Centroid conformer(s)

## Conformer Analysis

### `dm.conformers.rmsd(mol)`
Calculate pairwise RMSD matrix across all conformers.
- **Requirements**: Minimum 2 conformers
- **Returns**: NxN matrix of RMSD values
- **Use case**: Quantify conformer diversity

### `dm.conformers.sasa(mol, n_jobs=1, ...)`
Calculate Solvent Accessible Surface Area (SASA) using FreeSASA.
- **Parameters**:
  - `n_jobs`: Parallelization for multiple conformers
- **Returns**: Array of SASA values (one per conformer)
- **Storage**: Values stored in each conformer as property `'rdkit_free_sasa'`
- **Example**:
  ```python
  sasa_values = dm.conformers.sasa(mol_3d)
  # Or access from conformer properties
  conf = mol_3d.GetConformer(0)
  sasa = conf.GetDoubleProp('rdkit_free_sasa')
  ```

## Low-Level Conformer Manipulation

### `dm.conformers.center_of_mass(mol, conf_id=-1, use_atoms=True, round_coord=None)`
Calculate molecular center.
- **Parameters**:
  - `conf_id`: Conformer index (-1 for first conformer)
  - `use_atoms`: Use atomic masses (True) or geometric center (False)
  - `round_coord`: Decimal precision for rounding
- **Returns**: 3D coordinates of center
- **Use case**: Centering molecules for visualization or alignment

### `dm.conformers.get_coords(mol, conf_id=-1)`
Retrieve atomic coordinates from a conformer.
- **Returns**: Nx3 numpy array of atomic positions
- **Example**:
  ```python
  positions = dm.conformers.get_coords(mol_3d, conf_id=0)
  # positions.shape: (num_atoms, 3)
  ```

### `dm.conformers.translate(mol, conf_id=-1, transform_matrix=None)`
Reposition conformer using transformation matrix.
- **Modification**: Operates in-place
- **Use case**: Aligning or repositioning molecules

## Workflow Example

```python
import datamol as dm

# 1. Create molecule and generate conformers
mol = dm.to_mol("CC(C)CCO")  # Isopentanol
mol_3d = dm.conformers.generate(
    mol,
    n_confs=50,           # Generate 50 initial conformers
    rms_cutoff=0.5,       # Filter similar conformers
    minimize_energy=True   # Minimize energy
)

# 2. Analyze conformers
n_conformers = mol_3d.GetNumConformers()
print(f"Generated {n_conformers} unique conformers")

# 3. Calculate SASA
sasa_values = dm.conformers.sasa(mol_3d)

# 4. Cluster conformers
clusters = dm.conformers.cluster(mol_3d, rms_cutoff=1.0, centroids=False)

# 5. Get representative conformers
centroids = dm.conformers.return_centroids(mol_3d, clusters)

# 6. Access 3D coordinates
coords = dm.conformers.get_coords(mol_3d, conf_id=0)
```

## Key Concepts

- **Distance Geometry**: Method for generating 3D structures from connectivity information
- **ETKDG**: Uses experimental torsion angle preferences and additional chemical knowledge
- **RMS Cutoff**: Lower values = more unique conformers; higher values = fewer, more distinct conformers
- **Energy Minimization**: Relaxes structures to nearest local energy minimum
- **Hydrogens**: Critical for accurate 3D geometry - always include during embedding

# Datamol Core API Reference

This document covers the main functions available in the datamol namespace.

## Molecule Creation and Conversion

### `to_mol(mol, ...)`
Convert SMILES string or other molecular representations to RDKit molecule objects.
- **Parameters**: Accepts SMILES strings, InChI, or other molecular formats
- **Returns**: `rdkit.Chem.Mol` object
- **Common usage**: `mol = dm.to_mol("CCO")`

### `from_inchi(inchi)`
Convert InChI string to molecule object.

### `from_smarts(smarts)`
Convert SMARTS pattern to molecule object.

### `from_selfies(selfies)`
Convert SELFIES string to molecule object.

### `copy_mol(mol)`
Create a copy of a molecule object to avoid modifying the original.

## Molecule Export

### `to_smiles(mol, ...)`
Convert molecule object to SMILES string.
- **Common parameters**: `canonical=True`, `isomeric=True`

### `to_inchi(mol, ...)`
Convert molecule to InChI string representation.

### `to_inchikey(mol)`
Convert molecule to InChI key (fixed-length hash).

### `to_smarts(mol)`
Convert molecule to SMARTS pattern.

### `to_selfies(mol)`
Convert molecule to SELFIES (Self-Referencing Embedded Strings) format.

## Sanitization and Standardization

### `sanitize_mol(mol, ...)`
Enhanced version of RDKit's sanitize operation using mol→SMILES→mol conversion and aromatic nitrogen fixing.
- **Purpose**: Fix common molecular structure issues
- **Returns**: Sanitized molecule or None if sanitization fails

### `standardize_mol(mol, disconnect_metals=False, normalize=True, reionize=True, ...)`
Apply comprehensive standardization procedures including:
- Metal disconnection
- Normalization (charge corrections)
- Reionization
- Fragment handling (largest fragment selection)

### `standardize_smiles(smiles, ...)`
Apply SMILES standardization procedures directly to a SMILES string.

### `fix_mol(mol)`
Attempt to fix molecular structure issues automatically.

### `fix_valence(mol)`
Correct valence errors in molecular structures.

## Molecular Properties

### `reorder_atoms(mol, ...)`
Ensure consistent atom ordering for the same molecule regardless of original SMILES representation.
- **Purpose**: Maintain reproducible feature generation

### `remove_hs(mol, ...)`
Remove hydrogen atoms from molecular structure.

### `add_hs(mol, ...)`
Add explicit hydrogen atoms to molecular structure.

## Fingerprints and Similarity

### `to_fp(mol, fp_type='ecfp', ...)`
Generate molecular fingerprints for similarity calculations.
- **Fingerprint types**:
  - `'ecfp'` - Extended Connectivity Fingerprints (Morgan)
  - `'fcfp'` - Functional Connectivity Fingerprints
  - `'maccs'` - MACCS keys
  - `'topological'` - Topological fingerprints
  - `'atompair'` - Atom pair fingerprints
- **Common parameters**: `radius`, `fpSize` — extra kwargs are passed straight to RDKit's `rdFingerprintGenerator` (so `fpSize`, not `n_bits`/`nBits`); datamol's `ecfp` default is `radius=3` (ECFP6), `fpSize=2048`
- **Returns**: Numpy array (default `as_array=True`) or RDKit fingerprint object

### `pdist(mols, n_jobs=1, squareform=True, **fp_args)`
Calculate pairwise Tanimoto distances between all molecules in a list.
- **Supports**: Parallel processing via `n_jobs` parameter
- **Returns**: Square N×N distance matrix by default; `squareform=False` gives the condensed vector

### `cdist(mols1, mols2, ...)`
Calculate Tanimoto distances between two sets of molecules.

## Clustering and Diversity

### `cluster_mols(mols, cutoff=0.2, feature_fn=None, n_jobs=1)`
Cluster molecules using Butina clustering algorithm.
- **Parameters**:
  - `cutoff`: Distance threshold (default 0.2)
  - `feature_fn`: Custom function for molecular features
  - `n_jobs`: Parallelization (-1 for all cores)
- **Important**: Builds full distance matrix - suitable for ~1000 structures, not for 10,000+
- **Returns**: A tuple `(cluster_indices, cluster_mols)` — one tuple of indices and one list of Mols per cluster

### `pick_diverse(mols, npick, ...)`
Select diverse subset of molecules based on fingerprint diversity (MaxMin).
- **Returns**: `(picked_indices, picked_mols)`

### `pick_centroids(mols, npick=0, threshold=0.5, method='sphere', ...)`
Select centroid molecules representing clusters (sphere exclusion by default; `maxmin` or RDKit hierarchical methods also accepted).
- **Returns**: `(centroid_indices, centroid_mols)`

## Graph Operations

### `to_graph(mol)`
Convert molecule to graph representation for graph-based analysis.

### `get_all_path_between(mol, start, end)`
Find all paths between two atoms in molecular structure.

## DataFrame Integration

### `to_df(mols, smiles_column='smiles', mol_column='mol')`
Convert list of molecules to pandas DataFrame.

### `from_df(df, smiles_column='smiles', mol_column='mol')`
Convert pandas DataFrame to list of molecules.

# Datamol Descriptors and Visualization Reference

## Descriptors Module (`datamol.descriptors`)

The descriptors module provides tools for computing molecular properties and descriptors.

### Specialized Descriptor Functions

#### `dm.descriptors.n_aromatic_atoms(mol)`
Calculate the number of aromatic atoms.
- **Returns**: Integer count
- **Use case**: Aromaticity analysis

#### `dm.descriptors.n_aromatic_atoms_proportion(mol)`
Calculate ratio of aromatic atoms to total heavy atoms.
- **Returns**: Float between 0 and 1
- **Use case**: Quantifying aromatic character

#### `dm.descriptors.n_charged_atoms(mol)`
Count atoms with nonzero formal charge.
- **Returns**: Integer count
- **Use case**: Charge distribution analysis

#### `dm.descriptors.n_rigid_bonds(mol)`
Count non-rotatable bonds (neither single bonds nor ring bonds).
- **Returns**: Integer count
- **Use case**: Molecular flexibility assessment

#### `dm.descriptors.n_stereo_centers(mol)`
Count stereogenic centers (chiral centers).
- **Returns**: Integer count
- **Use case**: Stereochemistry analysis

#### `dm.descriptors.n_stereo_centers_unspecified(mol)`
Count stereocenters lacking stereochemical specification.
- **Returns**: Integer count
- **Use case**: Identifying incomplete stereochemistry

### Batch Descriptor Computation

#### `dm.descriptors.compute_many_descriptors(mol, properties_fn=None, add_properties=True)`
Compute multiple molecular properties for a single molecule.
- **Parameters**:
  - `properties_fn`: Custom list of descriptor functions
  - `add_properties`: Include additional computed properties
- **Returns**: Dictionary of descriptor name → value pairs (~22 keys)
- **Default keys** (datamol's own names, verified on 0.13.0): `mw`, `fsp3`,
  `n_lipinski_hba`, `n_lipinski_hbd`, `n_rings`, `n_hetero_atoms`,
  `n_heavy_atoms`, `n_rotatable_bonds`, `n_radical_electrons`, `tpsa`, `qed`,
  `clogp`, `sas`, plus aliphatic/aromatic/saturated carbocycle/heterocycle/ring
  counts (`n_aromatic_heterocycles` etc. — 0.13 fixed the old `..._heterocyles`
  misspelling, so the dict keys changed).
- **Naming gotcha**: logP is `clogp`; H-bond donors/acceptors are
  `n_lipinski_hbd` / `n_lipinski_hba`. There is **no** `logp`, `hbd`, `hba`, or
  `n_aromatic_atoms` key (use the standalone `n_aromatic_atoms(mol)` function for
  that count).
- **Example**:
  ```python
  mol = dm.to_mol("CCO")
  descriptors = dm.descriptors.compute_many_descriptors(mol)
  # {'mw': 46.04, 'clogp': -0.0, 'n_lipinski_hbd': 1, 'n_lipinski_hba': 1,
  #  'tpsa': 20.23, ...}
  ```

#### `dm.descriptors.batch_compute_many_descriptors(mols, properties_fn=None, add_properties=True, n_jobs=1, batch_size=None, progress=False)`
Compute descriptors for multiple molecules in parallel.
- **Parameters**:
  - `mols`: List of molecules
  - `n_jobs`: Number of parallel jobs (-1 for all cores)
  - `batch_size`: Chunk size for parallel processing — set it explicitly (e.g. `256`) whenever `n_jobs != 1`; the `None` default is rejected by joblib ≥ 1.6
  - `progress`: Show progress bar
- **Returns**: Pandas DataFrame with one row per molecule
- **Example**:
  ```python
  mols = [dm.to_mol(smi) for smi in smiles_list]
  df = dm.descriptors.batch_compute_many_descriptors(
      mols,
      n_jobs=-1,
      batch_size=256,
      progress=True
  )
  ```

### RDKit Descriptor Access

#### `dm.descriptors.any_rdkit_descriptor(name)`
Retrieve any descriptor function from RDKit by name.
- **Parameters**: `name` - Descriptor function name (e.g., 'MolWt', 'TPSA')
- **Returns**: RDKit descriptor function
- **Available descriptors**: From `rdkit.Chem.Descriptors` and `rdkit.Chem.rdMolDescriptors`
- **Example**:
  ```python
  tpsa_fn = dm.descriptors.any_rdkit_descriptor('TPSA')
  tpsa_value = tpsa_fn(mol)
  ```

### Common Use Cases

**Drug-likeness Filtering (Lipinski's Rule of Five)**:
```python
descriptors = dm.descriptors.compute_many_descriptors(mol)
is_druglike = (
    descriptors['mw'] <= 500 and
    descriptors['clogp'] <= 5 and
    descriptors['n_lipinski_hbd'] <= 5 and
    descriptors['n_lipinski_hba'] <= 10
)
```

**ADME Property Analysis**:
```python
df = dm.descriptors.batch_compute_many_descriptors(compound_library)
# Filter by TPSA for blood-brain barrier penetration
bbb_candidates = df[df['tpsa'] < 90]
```

---

## Visualization Module (`datamol.viz`)

The viz module provides tools for rendering molecules and conformers as images.

### Main Visualization Function

#### `dm.viz.to_image(mols, legends=None, n_cols=4, use_svg=True, mol_size=(300, 300), highlight_atom=None, highlight_bond=None, outfile=None, max_mols=32, copy=True, indices=False, align=False, ...)`
Generate image grid from molecules.
- **Parameters**:
  - `mols`: Single molecule or list of molecules
  - `legends`: String or list of strings as labels (one per molecule)
  - `n_cols`: Number of molecules per row (default: 4)
  - `use_svg`: Output SVG (True, **default**) or PNG (False) — set `use_svg=False` when writing a `.png` file
  - `mol_size`: Tuple (width, height) or single int for square images
  - `highlight_atom`: Atom indices to highlight (list or dict)
  - `highlight_bond`: Bond indices to highlight (list or dict)
  - `outfile`: Save path (local or remote, supports fsspec)
  - `max_mols`: Maximum number of molecules to display
  - `indices`: Draw atom indices on structures (default: False)
  - `align`: Align molecules using MCS (Maximum Common Substructure)
- **Returns**: Image object (can be displayed in Jupyter) or saves to file
- **Example**:
  ```python
  # Basic grid
  dm.viz.to_image(mols[:10], legends=[dm.to_smiles(m) for m in mols[:10]])

  # Save to file (PNG needs use_svg=False; SVG is the default)
  dm.viz.to_image(mols, outfile="molecules.png", n_cols=5, use_svg=False)

  # Highlight substructure
  dm.viz.to_image(mol, highlight_atom=[0, 1, 2], highlight_bond=[0, 1])

  # Aligned visualization
  dm.viz.to_image(mols, align=True, legends=activity_labels)
  ```

### Conformer Visualization

#### `dm.viz.conformers(mol, n_confs=None, align_conf=True, n_cols=3, sync_views=True, remove_hs=True, ...)`
Display multiple conformers in grid layout.
- **Parameters**:
  - `mol`: Molecule with embedded conformers
  - `n_confs`: Number or list of conformer indices to display (None = all)
  - `align_conf`: Align conformers for comparison (default: True)
  - `n_cols`: Grid columns (default: 3)
  - `sync_views`: Synchronize 3D views when interactive (default: True)
  - `remove_hs`: Remove hydrogens for clarity (default: True)
- **Returns**: Grid of conformer visualizations
- **Use case**: Comparing conformational diversity
- **Example**:
  ```python
  mol_3d = dm.conformers.generate(mol, n_confs=20)
  dm.viz.conformers(mol_3d, n_confs=10, align_conf=True)
  ```

### Circle Grid Visualization

#### `dm.viz.circle_grid(center_mol, circle_mols, mol_size=200, circle_margin=50, act_mapper=None, ...)`
Create concentric ring visualization with central molecule.
- **Parameters**:
  - `center_mol`: Molecule at center
  - `circle_mols`: List of molecule lists (one list per ring)
  - `mol_size`: Image size per molecule
  - `circle_margin`: Spacing between rings (default: 50)
  - `act_mapper`: Activity mapping dictionary for color-coding
- **Returns**: Circular grid image
- **Use case**: Visualizing molecular neighborhoods, SAR analysis, similarity networks
- **Example**:
  ```python
  # Show a reference molecule surrounded by similar compounds
  dm.viz.circle_grid(
      center_mol=reference,
      circle_mols=[nearest_neighbors, second_tier]
  )
  ```

### Visualization Best Practices

1. **Use legends for clarity**: Always label molecules with SMILES, IDs, or activity values
2. **Align related molecules**: Use `align=True` in `to_image()` for SAR analysis
3. **Adjust grid size**: Set `n_cols` based on molecule count and display width
4. **Use SVG for publications**: Set `use_svg=True` for scalable vector graphics
5. **Highlight substructures**: Use `highlight_atom` and `highlight_bond` to emphasize features
6. **Save large grids**: Use `outfile` parameter to save rather than display in memory

# Datamol Fragments and Scaffolds Reference

## Scaffolds Module (`datamol.scaffold`)

Scaffolds represent the core structure of molecules, useful for identifying structural families and analyzing structure-activity relationships (SAR).

### Murcko Scaffolds

#### `dm.to_scaffold_murcko(mol)`
Extract Bemis-Murcko scaffold (molecular framework).
- **Method**: Removes side chains, retaining ring systems and linkers
- **Returns**: Molecule object representing the scaffold
- **Use case**: Identify core structures across compound series
- **Example**:
  ```python
  mol = dm.to_mol("c1ccc(cc1)CCN")  # Phenethylamine
  scaffold = dm.to_scaffold_murcko(mol)
  scaffold_smiles = dm.to_smiles(scaffold)
  # Returns: 'c1ccccc1' — the ethylamine chain is a side chain (not a linker
  # between two ring systems), so it is stripped, leaving only the ring.
  ```

**Workflow for scaffold analysis**:
```python
# Extract scaffolds from compound library
scaffolds = [dm.to_scaffold_murcko(mol) for mol in mols]
scaffold_smiles = [dm.to_smiles(s) for s in scaffolds]

# Count scaffold frequency
from collections import Counter
scaffold_counts = Counter(scaffold_smiles)
most_common = scaffold_counts.most_common(10)
```

### Fuzzy Scaffolds

#### `dm.scaffold.fuzzy_scaffolding(mol, ...)`
Generate fuzzy scaffolds with enforceable groups that must appear in the core.
- **Purpose**: More flexible scaffold definition allowing specified functional groups
- **Use case**: Custom scaffold definitions beyond Murcko rules

### Applications

**Scaffold-based splitting** (for ML model validation):
```python
# Group compounds by scaffold
scaffold_to_mols = {}
for mol, scaffold in zip(mols, scaffolds):
    smi = dm.to_smiles(scaffold)
    if smi not in scaffold_to_mols:
        scaffold_to_mols[smi] = []
    scaffold_to_mols[smi].append(mol)

# Ensure train/test sets have different scaffolds
```

**SAR analysis**:
```python
# Group by scaffold and analyze activity
for scaffold_smi, molecules in scaffold_to_mols.items():
    activities = [get_activity(mol) for mol in molecules]
    print(f"Scaffold: {scaffold_smi}, Mean activity: {np.mean(activities)}")
```

---

## Fragments Module (`datamol.fragment`)

Molecular fragmentation breaks molecules into smaller pieces based on chemical rules, useful for fragment-based drug design and substructure analysis.

### BRICS Fragmentation

#### `dm.fragment.brics(mol, ...)`
Fragment molecule using BRICS (Breaking Retrosynthetically Interesting Chemical Substructures).
- **Method**: Dissects based on 16 chemically meaningful bond types
- **Consideration**: Considers chemical environment and surrounding substructures
- **Returns**: List of RDKit `Mol` fragments — the parent molecule comes first unless `remove_parent=True`; the default `fix=True` caps dummy atoms, so pass `fix=False` to keep attachment points
- **Use case**: Retrosynthetic analysis, fragment-based design
- **Example**:
  ```python
  mol = dm.to_mol("c1ccccc1CCN")
  fragments = dm.fragment.brics(mol, remove_parent=True, fix=False)
  frag_smiles = {dm.to_smiles(f) for f in fragments}
  # e.g. '[4*]CCN', '[16*]c1ccccc1' — [n*] marks attachment points
  ```

### RECAP Fragmentation

#### `dm.fragment.recap(mol, ...)`
Fragment molecule using RECAP (Retrosynthetic Combinatorial Analysis Procedure).
- **Method**: Dissects based on 11 predefined bond types
- **Rules**:
  - Leaves alkyl groups smaller than 5 carbons intact
  - Preserves cyclic bonds
- **Returns**: List of RDKit `Mol` fragments (parent first unless `remove_parent=True`; `fix=False` keeps `*` attachment points)
- **Use case**: Combinatorial library design
- **Example**:
  ```python
  mol = dm.to_mol("CCCCCc1ccccc1")
  fragments = dm.fragment.recap(mol)
  ```

### MMPA Fragmentation

#### `dm.fragment.mmpa_frag(mol, ...)`
Fragment for Matched Molecular Pair Analysis.
- **Purpose**: Generate fragments suitable for identifying molecular pairs
- **Use case**: Analyzing how small structural changes affect properties
- **Example**:
  ```python
  fragments = dm.fragment.mmpa_frag(mol)
  # Used to find pairs of molecules differing by single transformation
  ```

### Comparison of Methods

| Method | Bond Types | Preserves Cycles | Best For |
|--------|-----------|------------------|----------|
| BRICS  | 16        | Yes              | Retrosynthetic analysis, fragment recombination |
| RECAP  | 11        | Yes              | Combinatorial library design |
| MMPA   | Variable  | Depends          | Structure-activity relationship analysis |

### Fragmentation Workflow

```python
import datamol as dm

# 1. Fragment a molecule
mol = dm.to_mol("CC(=O)Oc1ccccc1C(=O)O")  # Aspirin
brics_frags = dm.fragment.brics(mol)
recap_frags = dm.fragment.recap(mol)

# 2. Analyze fragment frequency across library (count SMILES, not Mol objects)
all_fragments = []
for mol in molecule_library:
    frags = dm.fragment.brics(mol, remove_parent=True, fix=False)
    all_fragments.extend(dm.to_smiles(f) for f in frags)

# 3. Identify common fragments
from collections import Counter
fragment_counts = Counter(all_fragments)
common_fragments = fragment_counts.most_common(20)

# 4. Convert fragments back to molecules (remove attachment points)
def clean_fragment(frag_smiles):
    # Cap every [n*] / * attachment point with hydrogen
    import re
    clean = re.sub(r"\[\d*\*\]|\*", "[H]", frag_smiles)
    return dm.to_mol(clean)
```

### Advanced: Fragment-Based Virtual Screening

```python
# Build fragment library from known actives (as canonical SMILES)
def brics_smiles(mol):
    return {dm.to_smiles(f) for f in dm.fragment.brics(mol, remove_parent=True, fix=False)}

active_fragments = set()
for active_mol in active_compounds:
    active_fragments.update(brics_smiles(active_mol))

# Screen compounds for presence of active fragments
def score_by_fragments(mol, fragment_set):
    mol_frags = brics_smiles(mol)
    overlap = mol_frags.intersection(fragment_set)
    return len(overlap) / len(mol_frags) if mol_frags else 0.0

# Score screening library
scores = [score_by_fragments(mol, active_fragments) for mol in screening_lib]
```

### Key Concepts

- **Attachment Points**: Marked with [1*], [2*], etc. in fragment SMILES
- **Retrosynthetic**: Fragmentation mimics synthetic disconnections
- **Chemically Meaningful**: Breaks occur at typical synthetic bonds
- **Recombination**: Fragments can theoretically be recombined into valid molecules

# Datamol I/O Module Reference

The `datamol.io` module provides comprehensive file handling for molecular data across multiple formats.

## Reading Molecular Files

### `dm.read_sdf(urlpath, sanitize=True, as_df=False, smiles_column='smiles', mol_column=None, remove_hs=True, n_jobs=1, ...)`
Read Structure-Data File (SDF) format.
- **Parameters**:
  - `filename`: Path to SDF file (supports local and remote paths via fsspec)
  - `sanitize`: Apply sanitization to molecules
  - `remove_hs`: Remove explicit hydrogens
  - `as_df`: Return a DataFrame (True) or a list of molecules (False, **default**)
  - `mol_column`: Name of molecule column in the DataFrame (only with `as_df=True`)
  - `n_jobs`: Enable parallel processing
- **Returns**: List of molecules, or a DataFrame with `as_df=True`
- **Example**: `mols = dm.read_sdf("compounds.sdf")`; `df = dm.read_sdf("compounds.sdf", as_df=True, mol_column="mol")`

### `dm.read_smi(urlpath)`
Read a SMILES file (SMILES optionally followed by an ID/name).
- **Returns**: List of molecules (no DataFrame option; wrap with `dm.to_df(mols)` if needed)
- **Example**: `mols = dm.read_smi("molecules.smi")`

### `dm.read_csv(filename, smiles_column='smiles', mol_column=None, ...)`
Read CSV file with optional automatic SMILES-to-molecule conversion.
- **Parameters**:
  - `smiles_column`: Column containing SMILES strings
  - `mol_column`: If specified, creates molecule objects from SMILES column
- **Example**: `df = dm.read_csv("data.csv", smiles_column="SMILES", mol_column="mol")`

### `dm.read_excel(filename, sheet_name=0, smiles_column='smiles', mol_column=None, ...)`
Read Excel files with molecule handling.
- **Parameters**:
  - `sheet_name`: Sheet to read (index or name)
  - Other parameters similar to `read_csv`
- **Example**: `df = dm.read_excel("compounds.xlsx", sheet_name="Sheet1")`

### `dm.read_molblock(molblock, sanitize=True, remove_hs=True)`
Parse MOL block string (molecular structure text representation).

### `dm.read_mol2file(filename, sanitize=True, remove_hs=True, cleanupSubstructures=True)`
Read Mol2 format files.

### `dm.read_pdbfile(filename, sanitize=True, remove_hs=True, proximityBonding=True)`
Read Protein Data Bank (PDB) format files.

### `dm.read_pdbblock(pdbblock, sanitize=True, remove_hs=True, proximityBonding=True)`
Parse PDB block string.

### `dm.open_df(filename, ...)`
Universal DataFrame reader - automatically detects format.
- **Supported formats**: CSV, Excel, Parquet, JSON, SDF
- **Example**: `df = dm.open_df("data.csv")` or `df = dm.open_df("molecules.sdf")`

## Writing Molecular Files

### `dm.to_sdf(mols, filename, mol_column=None, ...)`
Write molecules to SDF file.
- **Input types**:
  - List of molecules
  - DataFrame with molecule column
  - Sequence of molecules
- **Parameters**:
  - `mol_column`: Column name if input is DataFrame
- **Example**:
  ```python
  dm.to_sdf(mols, "output.sdf")
  # or from DataFrame
  dm.to_sdf(df, "output.sdf", mol_column="mol")
  ```

### `dm.to_smi(mols, filename, mol_column=None, ...)`
Write molecules to SMILES file with optional validation.
- **Format**: SMILES strings with optional molecule names/IDs

### `dm.to_xlsx(mols, urlpath, smiles_column='smiles', mol_column='mol', mol_size=[300, 300])`
Export molecules or a DataFrame to Excel with rendered molecular images.
- **Parameters**:
  - `mol_column`: The (single) column containing molecules to render as images
- **Special feature**: Automatically renders molecules as images in Excel cells
- **Example**: `dm.to_xlsx(df, "molecules.xlsx", mol_column="mol")`

### `dm.to_molblock(mol, ...)`
Convert molecule to MOL block string.

### `dm.to_pdbblock(mol, ...)`
Convert molecule to PDB block string.

### `dm.save_df(df, filename, ...)`
Save DataFrame in multiple formats (CSV, Excel, Parquet, JSON).

## Remote File Support

All I/O functions support remote file paths through fsspec integration:
- **Supported protocols**: S3 (AWS), GCS (Google Cloud), Azure, HTTP/HTTPS
- **Example**:
  ```python
  dm.read_sdf("s3://bucket/compounds.sdf")
  dm.read_csv("https://example.com/data.csv")
  ```

## Key Parameters Across Functions

- **`sanitize`**: Apply molecule sanitization (default: True)
- **`remove_hs`**: Remove explicit hydrogens (default: True)
- **`as_df`**: Return DataFrame vs list (`read_sdf` defaults to a list, `as_df=False`)
- **`n_jobs`**: Parallel processing (`-1` = all cores; `0`, `1`, or `None` = sequential)
- **`mol_column`**: Name of molecule column in DataFrames
- **`smiles_column`**: Name of SMILES column in DataFrames

# Datamol Reactions and Data Modules Reference

## Reactions Module (`datamol.reactions`)

The reactions module enables programmatic application of chemical transformations using SMARTS reaction patterns.

### Applying Chemical Reactions

#### `dm.reactions.apply_reaction(rxn, reactants, product_index=None, single_product_group=False, as_smiles=False, rm_attach=False, sanitize=True)`
Apply a chemical reaction to reactant molecules.
- **Parameters**:
  - `rxn`: Reaction object (from SMARTS pattern)
  - `reactants`: Tuple of reactant molecules
  - `as_smiles`: Return SMILES strings (True) or molecule objects (False)
  - `sanitize`: Sanitize product molecules
  - `single_product_group`: Return one product set (True) or all product sets (False, default)
  - `rm_attach`: Remove attachment point markers (default False)
  - `product_index`: Which product to return from each set (default None = all)
- **Returns**: With the defaults, a list of product lists; with `single_product_group=True, product_index=0`, a single Mol (or SMILES with `as_smiles=True`)
- **Example**:
  ```python
  from rdkit import Chem

  # Define reaction: alcohol + carboxylic acid → ester
  rxn = Chem.rdChemReactions.ReactionFromSmarts(
      '[C:1][OH:2].[C:3](=[O:4])[OH:5]>>[C:1][O:2][C:3](=[O:4])'
  )

  # Apply to reactants
  alcohol = dm.to_mol("CCO")
  acid = dm.to_mol("CC(=O)O")
  product = dm.reactions.apply_reaction(rxn, (alcohol, acid),
                                        single_product_group=True, product_index=0)
  ```

### Creating Reactions

Reactions are typically created from SMARTS patterns using RDKit:
```python
from rdkit.Chem import rdChemReactions

# Reaction pattern: [reactant1].[reactant2]>>[product]
rxn = rdChemReactions.ReactionFromSmarts(
    '[1*][*:1].[1*][*:2]>>[*:1][*:2]'
)
```

### Validation Functions

The module includes functions to:
- **Check if molecule is reactant**: Verify if molecule matches reactant pattern
- **Validate reaction**: Check if reaction is synthetically reasonable
- **Process reaction files**: Load reactions from files or databases

### Common Reaction Patterns

**Amide formation**:
```python
# Amine + carboxylic acid → amide
amide_rxn = rdChemReactions.ReactionFromSmarts(
    '[N:1].[C:2](=[O:3])[OH]>>[N:1][C:2](=[O:3])'
)
```

**Suzuki coupling**:
```python
# Aryl halide + boronic acid → biaryl
suzuki_rxn = rdChemReactions.ReactionFromSmarts(
    '[c:1][Br].[c:2][B]([OH])[OH]>>[c:1][c:2]'
)
```

**Functional group transformations**:
```python
# Alcohol → ester
esterification = rdChemReactions.ReactionFromSmarts(
    '[C:1][OH:2].[C:3](=[O:4])[Cl]>>[C:1][O:2][C:3](=[O:4])'
)
```

### Workflow Example

```python
import datamol as dm
from rdkit.Chem import rdChemReactions

# 1. Define reaction
rxn_smarts = '[C:1](=[O:2])[OH:3]>>[C:1](=[O:2])[Cl:3]'  # Acid → acid chloride
rxn = rdChemReactions.ReactionFromSmarts(rxn_smarts)

# 2. Apply to molecule library
acids = [dm.to_mol(smi) for smi in acid_smiles_list]
acid_chlorides = []

for acid in acids:
    try:
        product = dm.reactions.apply_reaction(
            rxn,
            (acid,),  # Single reactant as tuple
            single_product_group=True,
            product_index=0,
            sanitize=True
        )
        acid_chlorides.append(product)
    except Exception as e:
        print(f"Reaction failed: {e}")

# 3. Validate products
valid_products = [p for p in acid_chlorides if p is not None]
```

### Key Concepts

- **SMARTS**: SMiles ARbitrary Target Specification - pattern language for reactions
- **Atom Mapping**: Numbers like [C:1] preserve atom identity through reaction
- **Attachment Points**: [1*] represents generic connection points
- **Reaction Validation**: Not all SMARTS reactions are chemically reasonable

---

## Data Module (`datamol.data`)

The data module provides convenient access to curated molecular datasets for testing and learning.

### Available Datasets

#### `dm.data.cdk2(as_df=True, mol_column='mol')`
RDKit CDK2 dataset - kinase inhibitor data.
- **Parameters**:
  - `as_df`: Return as DataFrame (True) or list of molecules (False)
  - `mol_column`: Name for molecule column
- **Returns**: Dataset with molecular structures and activity data
- **Use case**: Small dataset for algorithm testing
- **Example**:
  ```python
  cdk2_df = dm.data.cdk2(as_df=True)
  print(cdk2_df.shape)
  print(cdk2_df.columns)
  ```

#### `dm.data.freesolv()`
FreeSolv dataset - experimental and calculated hydration free energies.
- **Contents**: 642 molecules with:
  - IUPAC names
  - SMILES strings
  - Experimental hydration free energy values
  - Calculated values
- **Warning**: "Only meant to be used as a toy dataset for pedagogic and testing purposes"
- **Not suitable for**: Benchmarking or production model training
- **Example**:
  ```python
  freesolv_df = dm.data.freesolv()
  # Columns: iupac, smiles, expt (kcal/mol), calc (kcal/mol)
  ```

#### `dm.data.solubility(as_df=True, mol_column='mol')`
RDKit solubility dataset with train/test splits.
- **Contents**: Aqueous solubility data with pre-defined splits
- **Columns**: Includes 'split' column with 'train' or 'test' values
- **Use case**: Testing ML workflows with proper train/test separation
- **Example**:
  ```python
  sol_df = dm.data.solubility(as_df=True)
  # Columns: mol, ID, NAME, SOL, SOL_classification, smiles, split
  # The numeric solubility target is the 'SOL' column.

  # Split into train/test
  train_df = sol_df[sol_df['split'] == 'train']
  test_df = sol_df[sol_df['split'] == 'test']

  # Featurize: to_fp takes ONE mol, so stack per-row into a matrix.
  import numpy as np
  X_train = np.array([dm.to_fp(m) for m in train_df['mol']])
  y_train = train_df['SOL'].values
  ```

### Usage Guidelines

**For testing and tutorials**:
```python
# Quick dataset for testing code
df = dm.data.cdk2()
mols = df['mol'].tolist()

# Test descriptor calculation
descriptors_df = dm.descriptors.batch_compute_many_descriptors(mols)

# Test clustering (returns (cluster_indices, cluster_mols))
cluster_idx, cluster_mols = dm.cluster_mols(mols, cutoff=0.3)
```

**For learning workflows**:
```python
# Complete ML pipeline example
sol_df = dm.data.solubility()

# Preprocessing
train = sol_df[sol_df['split'] == 'train']
test = sol_df[sol_df['split'] == 'test']

# Featurization: to_fp is per-molecule, so build the matrix row by row
import numpy as np
X_train = np.array([dm.to_fp(m) for m in train['mol']])
X_test = np.array([dm.to_fp(m) for m in test['mol']])

# Model training (example) — 'SOL' is the numeric solubility target column
from sklearn.ensemble import RandomForestRegressor
model = RandomForestRegressor()
model.fit(X_train, train['SOL'])
predictions = model.predict(X_test)
```

### Important Notes

- **Toy Datasets**: Designed for pedagogical purposes, not production use
- **Small Size**: Limited number of compounds suitable for quick tests
- **Pre-processed**: Data already cleaned and formatted
- **Citations**: Check dataset documentation for proper attribution if publishing

### Best Practices

1. **Use for development only**: Don't draw scientific conclusions from toy datasets
2. **Validate on real data**: Always test production code on actual project data
3. **Proper attribution**: Cite original data sources if using in publications
4. **Understand limitations**: Know the scope and quality of each dataset

# Datamol Worked Workflow Recipes

End-to-end, copy-ready pipelines and multi-step recipes extracted from the skill body so `SKILL.md` stays lean; each block is self-contained and assumes `import datamol as dm`.

## Complete Pipeline: Data Loading → Filtering → Analysis

```python
import datamol as dm
import pandas as pd

# 1. Load molecules (read_sdf returns a list unless as_df=True)
df = dm.read_sdf("compounds.sdf", as_df=True, mol_column="mol")

# 2. Standardize
df['mol'] = df['mol'].apply(lambda m: dm.standardize_mol(m) if m else None)
df = df[df['mol'].notna()]  # Remove failed molecules

# 3. Compute descriptors
desc_df = dm.descriptors.batch_compute_many_descriptors(
    df['mol'].tolist(),
    n_jobs=-1,
    batch_size=256,   # required with n_jobs != 1 under joblib >= 1.6
    progress=True
)

# 4. Filter by drug-likeness (batch_compute_many_descriptors uses the same keys
#    as compute_many_descriptors: clogp, n_lipinski_hbd, n_lipinski_hba)
druglike = (
    (desc_df['mw'] <= 500) &
    (desc_df['clogp'] <= 5) &
    (desc_df['n_lipinski_hbd'] <= 5) &
    (desc_df['n_lipinski_hba'] <= 10)
)
filtered_df = df[druglike.values]

# 5. Cluster and select diverse subset (returns (indices, mols))
diverse_idx, diverse_mols = dm.pick_diverse(
    filtered_df['mol'].tolist(),
    npick=100
)

# 6. Visualize results (PNG output needs use_svg=False)
dm.viz.to_image(
    diverse_mols,
    legends=[dm.to_smiles(m) for m in diverse_mols],
    outfile="diverse_compounds.png",
    use_svg=False,
    n_cols=10
)
```

## Structure-Activity Relationship (SAR) Analysis

```python
# Group by scaffold
scaffolds = [dm.to_scaffold_murcko(mol) for mol in mols]
scaffold_smiles = [dm.to_smiles(s) for s in scaffolds]

# Create DataFrame with activities
sar_df = pd.DataFrame({
    'mol': mols,
    'scaffold': scaffold_smiles,
    'activity': activities  # User-provided activity data
})

# Analyze each scaffold series
for scaffold, group in sar_df.groupby('scaffold'):
    if len(group) >= 3:  # Need multiple examples
        print(f"\nScaffold: {scaffold}")
        print(f"Count: {len(group)}")
        print(f"Activity range: {group['activity'].min():.2f} - {group['activity'].max():.2f}")

        # Visualize with activities as legends
        dm.viz.to_image(
            group['mol'].tolist(),
            legends=[f"Activity: {act:.2f}" for act in group['activity']],
            align=True  # Align by common substructure
        )
```

## Virtual Screening Pipeline

```python
# 1. Calculate Tanimoto distances between query actives and the library.
#    dm.cdist takes the molecules directly (it fingerprints internally),
#    returning an (n_query, n_library) distance matrix.
import numpy as np

distances = dm.cdist(query_actives, library_mols, n_jobs=-1)

# 3. Find closest matches (min distance to any query)
min_distances = distances.min(axis=0)
similarities = 1 - min_distances  # Convert distance to similarity

# 4. Rank and select top hits
top_indices = np.argsort(similarities)[::-1][:100]  # Top 100
top_hits = [library_mols[i] for i in top_indices]
top_scores = [similarities[i] for i in top_indices]

# 5. Visualize hits
dm.viz.to_image(
    top_hits[:20],
    legends=[f"Sim: {score:.3f}" for score in top_scores[:20]],
    outfile="screening_hits.png"
)
```

## Scaffold-Based Analysis

```python
# Group compounds by scaffold
from collections import Counter

scaffolds = [dm.to_scaffold_murcko(mol) for mol in mols]
scaffold_smiles = [dm.to_smiles(s) for s in scaffolds]

# Count scaffold frequency
scaffold_counts = Counter(scaffold_smiles)
most_common = scaffold_counts.most_common(10)

# Create scaffold-to-molecules mapping
scaffold_groups = {}
for mol, scaf_smi in zip(mols, scaffold_smiles):
    if scaf_smi not in scaffold_groups:
        scaffold_groups[scaf_smi] = []
    scaffold_groups[scaf_smi].append(mol)
```

## Scaffold-Based Train/Test Splitting (for ML)

```python
# Ensure train and test sets have different scaffolds
scaffold_to_mols = {}
for mol, scaf in zip(mols, scaffold_smiles):
    if scaf not in scaffold_to_mols:
        scaffold_to_mols[scaf] = []
    scaffold_to_mols[scaf].append(mol)

# Split scaffolds into train/test
import random
scaffolds = list(scaffold_to_mols.keys())
random.shuffle(scaffolds)
split_idx = int(0.8 * len(scaffolds))
train_scaffolds = scaffolds[:split_idx]
test_scaffolds = scaffolds[split_idx:]

# Get molecules for each split
train_mols = [mol for scaf in train_scaffolds for mol in scaffold_to_mols[scaf]]
test_mols = [mol for scaf in test_scaffolds for mol in scaffold_to_mols[scaf]]
```

## Fragment Frequency and Overlap Scoring

```python
# Find common fragments across compound library
from collections import Counter

def brics_smiles(mol):
    # brics() returns Mol objects (parent first); compare fragments as SMILES
    return {dm.to_smiles(f) for f in dm.fragment.brics(mol, remove_parent=True, fix=False)}

all_fragments = []
for mol in mols:
    all_fragments.extend(brics_smiles(mol))

fragment_counts = Counter(all_fragments)
common_frags = fragment_counts.most_common(20)

# Fragment-based scoring
def fragment_score(mol, reference_fragments):
    mol_frags = brics_smiles(mol)
    overlap = mol_frags.intersection(reference_fragments)
    return len(overlap) / len(mol_frags) if mol_frags else 0
```

## Integration with Machine Learning

```python
# Feature generation
X = np.array([dm.to_fp(mol) for mol in mols])

# Or descriptors
desc_df = dm.descriptors.batch_compute_many_descriptors(mols, n_jobs=-1, batch_size=256)
X = desc_df.values

# Train model
from sklearn.ensemble import RandomForestRegressor
model = RandomForestRegressor()
model.fit(X, y_target)

# Predict
predictions = model.predict(X_test)
```

## Robust Error Handling Wrappers

```python
# Safe molecule creation
def safe_to_mol(smiles):
    try:
        mol = dm.to_mol(smiles)
        if mol is not None:
            mol = dm.standardize_mol(mol)
        return mol
    except Exception as e:
        print(f"Failed to process {smiles}: {e}")
        return None

# Safe batch processing
valid_mols = []
for smiles in smiles_list:
    mol = safe_to_mol(smiles)
    if mol is not None:
        valid_mols.append(mol)
```