Datamol

{
  "description": "\"Pythonic wrapper around RDKit with simplified interface and sensible defaults. Preferred for standard drug discovery: SMILES parsing, standardization, descriptors, fingerprints, clustering, 3D conformers, parallel processing. Returns native rdkit.Chem.Mol objects. For advanced control or custom parameters, use rdkit directly.\"",
  "references": {
    "files": [
      "references/conformers_module.md",
      "references/core_api.md",
      "references/descriptors_viz.md",
      "references/fragments_scaffolds.md",
      "references/io_module.md",
      "references/reactions_data.md"
    ]
  },
  "content": "Guide users to install datamol:\r\n\r\n```bash\r\nconda install -c conda-forge datamol\r\n\r\n\r\n### 1. Basic Molecule Handling\r\n\r\n**Creating molecules from SMILES**:\r\n```python\r\nimport datamol as dm\r\n\r\nmol = dm.to_mol(\"CCO\")  # Ethanol\r\n\r\nsmiles_list = [\"CCO\", \"c1ccccc1\", \"CC(=O)O\"]\r\nmols = [dm.to_mol(smi) for smi in smiles_list]\r\n\r\nmol = dm.to_mol(\"invalid_smiles\")  # Returns None\r\nif mol is None:\r\n    print(\"Failed to parse SMILES\")\r\n```\r\n\r\n**Converting molecules to SMILES**:\r\n```python\r\nsmiles = dm.to_smiles(mol)\r\n\r\nsmiles = dm.to_smiles(mol, isomeric=True)\r\n\r\ninchi = dm.to_inchi(mol)\r\ninchikey = dm.to_inchikey(mol)\r\nselfies = dm.to_selfies(mol)\r\n```\r\n\r\n**Standardization and sanitization** (always recommend for user-provided molecules):\r\n```python\r\nmol = dm.sanitize_mol(mol)\r\n\r\nmol = dm.standardize_mol(\r\n    mol,\r\n    disconnect_metals=True,\r\n    normalize=True,\r\n    reionize=True\r\n)\r\n\r\nclean_smiles = dm.standardize_smiles(smiles)\r\n```\r\n\r\n### 2. Reading and Writing Molecular Files\r\n\r\nRefer to `references/io_module.md` for comprehensive I/O documentation.\r\n\r\n**Reading files**:\r\n```python\r\ndf = dm.read_sdf(\"compounds.sdf\", mol_column='mol')\r\n\r\ndf = dm.read_smi(\"molecules.smi\", smiles_column='smiles', mol_column='mol')\r\n\r\ndf = dm.read_csv(\"data.csv\", smiles_column=\"SMILES\", mol_column=\"mol\")\r\n\r\ndf = dm.read_excel(\"compounds.xlsx\", sheet_name=0, mol_column=\"mol\")\r\n\r\ndf = dm.open_df(\"file.sdf\")  # Works with .sdf, .csv, .xlsx, .parquet, .json\r\n```\r\n\r\n**Writing files**:\r\n```python\r\ndm.to_sdf(mols, \"output.sdf\")\r\ndm.to_sdf(df, \"output.sdf\", mol_column=\"mol\")\r\n\r\ndm.to_smi(mols, \"output.smi\")\r\n\r\ndm.to_xlsx(df, \"output.xlsx\", mol_columns=[\"mol\"])\r\n```\r\n\r\n**Remote file support** (S3, GCS, HTTP):\r\n```python\r\ndf = dm.read_sdf(\"s3://bucket/compounds.sdf\")\r\ndf = dm.read_csv(\"https://example.com/data.csv\")\r\n\r\ndm.to_sdf(mols, \"s3://bucket/output.sdf\")\r\n```\r\n\r\n### 3. Molecular Descriptors and Properties\r\n\r\nRefer to `references/descriptors_viz.md` for detailed descriptor documentation.\r\n\r\n**Computing descriptors for a single molecule**:\r\n```python\r\ndescriptors = dm.descriptors.compute_many_descriptors(mol)\r\n```\r\n\r\n**Batch descriptor computation** (recommended for datasets):\r\n```python\r\ndesc_df = dm.descriptors.batch_compute_many_descriptors(\r\n    mols,\r\n    n_jobs=-1,      # Use all CPU cores\r\n    progress=True   # Show progress bar\r\n)\r\n```\r\n\r\n**Specific descriptors**:\r\n```python\r\nn_aromatic = dm.descriptors.n_aromatic_atoms(mol)\r\naromatic_ratio = dm.descriptors.n_aromatic_atoms_proportion(mol)\r\n\r\nn_stereo = dm.descriptors.n_stereo_centers(mol)\r\nn_unspec = dm.descriptors.n_stereo_centers_unspecified(mol)\r\n\r\nn_rigid = dm.descriptors.n_rigid_bonds(mol)\r\n```\r\n\r\n**Drug-likeness filtering (Lipinski's Rule of Five)**:\r\n```python\r\ndef is_druglike(mol):\r\n    desc = dm.descriptors.compute_many_descriptors(mol)\r\n    return (\r\n        desc['mw'] <= 500 and\r\n        desc['logp'] <= 5 and\r\n        desc['hbd'] <= 5 and\r\n        desc['hba'] <= 10\r\n    )\r\n\r\ndruglike_mols = [mol for mol in mols if is_druglike(mol)]\r\n```\r\n\r\n### 4. Molecular Fingerprints and Similarity\r\n\r\n**Generating fingerprints**:\r\n```python\r\nfp = dm.to_fp(mol, fp_type='ecfp', radius=2, n_bits=2048)\r\n\r\nfp_maccs = dm.to_fp(mol, fp_type='maccs')\r\nfp_topological = dm.to_fp(mol, fp_type='topological')\r\nfp_atompair = dm.to_fp(mol, fp_type='atompair')\r\n```\r\n\r\n**Similarity calculations**:\r\n```python\r\ndistance_matrix = dm.pdist(mols, n_jobs=-1)\r\n\r\ndistances = dm.cdist(query_mols, library_mols, n_jobs=-1)\r\n\r\nfrom scipy.spatial.distance import squareform\r\ndist_matrix = squareform(dm.pdist(mols))\r\n```\r\n\r\n### 5. Clustering and Diversity Selection\r\n\r\nRefer to `references/core_api.md` for clustering details.\r\n\r\n**Butina clustering**:\r\n```python\r\nclusters = dm.cluster_mols(\r\n    mols,\r\n    cutoff=0.2,    # Tanimoto distance threshold (0=identical, 1=completely different)\r\n    n_jobs=-1      # Parallel processing\r\n)\r\n\r\nfor i, cluster in enumerate(clusters):\r\n    print(f\"Cluster {i}: {len(cluster)} molecules\")\r\n    cluster_mols = [mols[idx] for idx in cluster]\r\n```\r\n\r\n**Important**: Butina clustering builds a full distance matrix - suitable for ~1000 molecules, not for 10,000+.\r\n\r\n**Diversity selection**:\r\n```python\r\ndiverse_mols = dm.pick_diverse(\r\n    mols,\r\n    npick=100  # Select 100 diverse molecules\r\n)\r\n\r\ncentroids = dm.pick_centroids(\r\n    mols,\r\n    npick=50   # Select 50 representative molecules\r\n)\r\n```\r\n\r\n### 6. Scaffold Analysis\r\n\r\nRefer to `references/fragments_scaffolds.md` for complete scaffold documentation.\r\n\r\n**Extracting Murcko scaffolds**:\r\n```python\r\nscaffold = dm.to_scaffold_murcko(mol)\r\nscaffold_smiles = dm.to_smiles(scaffold)\r\n```\r\n\r\n**Scaffold-based analysis**:\r\n```python\r\nfrom collections import Counter\r\n\r\nscaffolds = [dm.to_scaffold_murcko(mol) for mol in mols]\r\nscaffold_smiles = [dm.to_smiles(s) for s in scaffolds]\r\n\r\nscaffold_counts = Counter(scaffold_smiles)\r\nmost_common = scaffold_counts.most_common(10)\r\n\r\nscaffold_groups = {}\r\nfor mol, scaf_smi in zip(mols, scaffold_smiles):\r\n    if scaf_smi not in scaffold_groups:\r\n        scaffold_groups[scaf_smi] = []\r\n    scaffold_groups[scaf_smi].append(mol)\r\n```\r\n\r\n**Scaffold-based train/test splitting** (for ML):\r\n```python\r\nscaffold_to_mols = {}\r\nfor mol, scaf in zip(mols, scaffold_smiles):\r\n    if scaf not in scaffold_to_mols:\r\n        scaffold_to_mols[scaf] = []\r\n    scaffold_to_mols[scaf].append(mol)\r\n\r\nimport random\r\nscaffolds = list(scaffold_to_mols.keys())\r\nrandom.shuffle(scaffolds)\r\nsplit_idx = int(0.8 * len(scaffolds))\r\ntrain_scaffolds = scaffolds[:split_idx]\r\ntest_scaffolds = scaffolds[split_idx:]\r\n\r\ntrain_mols = [mol for scaf in train_scaffolds for mol in scaffold_to_mols[scaf]]\r\ntest_mols = [mol for scaf in test_scaffolds for mol in scaffold_to_mols[scaf]]\r\n```\r\n\r\n### 7. Molecular Fragmentation\r\n\r\nRefer to `references/fragments_scaffolds.md` for fragmentation details.\r\n\r\n**BRICS fragmentation** (16 bond types):\r\n```python\r\nfragments = dm.fragment.brics(mol)\r\n```\r\n\r\n**RECAP fragmentation** (11 bond types):\r\n```python\r\nfragments = dm.fragment.recap(mol)\r\n```\r\n\r\n**Fragment analysis**:\r\n```python\r\nfrom collections import Counter\r\n\r\nall_fragments = []\r\nfor mol in mols:\r\n    frags = dm.fragment.brics(mol)\r\n    all_fragments.extend(frags)\r\n\r\nfragment_counts = Counter(all_fragments)\r\ncommon_frags = fragment_counts.most_common(20)\r\n\r\ndef fragment_score(mol, reference_fragments):\r\n    mol_frags = dm.fragment.brics(mol)\r\n    overlap = mol_frags.intersection(reference_fragments)\r\n    return len(overlap) / len(mol_frags) if mol_frags else 0\r\n```\r\n\r\n### 8. 3D Conformer Generation\r\n\r\nRefer to `references/conformers_module.md` for detailed conformer documentation.\r\n\r\n**Generating conformers**:\r\n```python\r\nmol_3d = dm.conformers.generate(\r\n    mol,\r\n    n_confs=50,           # Number to generate (auto if None)\r\n    rms_cutoff=0.5,       # Filter similar conformers (Ångströms)\r\n    minimize_energy=True,  # Minimize with UFF force field\r\n    method='ETKDGv3'      # Embedding method (recommended)\r\n)\r\n\r\nn_conformers = mol_3d.GetNumConformers()\r\nconf = mol_3d.GetConformer(0)  # Get first conformer\r\npositions = conf.GetPositions()  # Nx3 array of atom coordinates\r\n```\r\n\r\n**Conformer clustering**:\r\n```python\r\nclusters = dm.conformers.cluster(\r\n    mol_3d,\r\n    rms_cutoff=1.0,\r\n    centroids=False\r\n)\r\n\r\ncentroids = dm.conformers.return_centroids(mol_3d, clusters)\r\n```\r\n\r\n**SASA calculation**:\r\n```python\r\nsasa_values = dm.conformers.sasa(mol_3d, n_jobs=-1)\r\n\r\nconf = mol_3d.GetConformer(0)\r\nsasa = conf.GetDoubleProp('rdkit_free_sasa')\r\n```\r\n\r\n### 9. Visualization\r\n\r\nRefer to `references/descriptors_viz.md` for visualization documentation.\r\n\r\n**Basic molecule grid**:\r\n```python\r\ndm.viz.to_image(\r\n    mols[:20],\r\n    legends=[dm.to_smiles(m) for m in mols[:20]],\r\n    n_cols=5,\r\n    mol_size=(300, 300)\r\n)\r\n\r\ndm.viz.to_image(mols, outfile=\"molecules.png\")\r\n\r\ndm.viz.to_image(mols, outfile=\"molecules.svg\", use_svg=True)\r\n```\r\n\r\n**Aligned visualization** (for SAR analysis):\r\n```python\r\ndm.viz.to_image(\r\n    similar_mols,\r\n    align=True,  # Enable MCS alignment\r\n    legends=activity_labels,\r\n    n_cols=4\r\n)\r\n```\r\n\r\n**Highlighting substructures**:\r\n```python\r\ndm.viz.to_image(\r\n    mol,\r\n    highlight_atom=[0, 1, 2, 3],  # Atom indices\r\n    highlight_bond=[0, 1, 2]      # Bond indices\r\n)\r\n```\r\n\r\n**Conformer visualization**:\r\n```python\r\ndm.viz.conformers(\r\n    mol_3d,\r\n    n_confs=10,\r\n    align_conf=True,\r\n    n_cols=3\r\n)\r\n```\r\n\r\n### 10. Chemical Reactions\r\n\r\nRefer to `references/reactions_data.md` for reactions documentation.\r\n\r\n**Applying reactions**:\r\n```python\r\nfrom rdkit.Chem import rdChemReactions\r\n\r\nrxn_smarts = '[C:1](=[O:2])[OH:3]>>[C:1](=[O:2])[Cl:3]'\r\nrxn = rdChemReactions.ReactionFromSmarts(rxn_smarts)\r\n\r\nreactant = dm.to_mol(\"CC(=O)O\")  # Acetic acid\r\nproduct = dm.reactions.apply_reaction(\r\n    rxn,\r\n    (reactant,),\r\n    sanitize=True\r\n)\r\n\r\nproduct_smiles = dm.to_smiles(product)\r\n```\r\n\r\n**Batch reaction application**:\r\n```python\r\n\r\n### Complete Pipeline: Data Loading → Filtering → Analysis\r\n\r\n```python\r\nimport datamol as dm\r\nimport pandas as pd\r\n\r\ndf = dm.read_sdf(\"compounds.sdf\")\r\n\r\ndf['mol'] = df['mol'].apply(lambda m: dm.standardize_mol(m) if m else None)\r\ndf = df[df['mol'].notna()]  # Remove failed molecules\r\n\r\ndesc_df = dm.descriptors.batch_compute_many_descriptors(\r\n    df['mol'].tolist(),\r\n    n_jobs=-1,\r\n    progress=True\r\n)\r\n\r\ndruglike = (\r\n    (desc_df['mw'] <= 500) &\r\n    (desc_df['logp'] <= 5) &\r\n    (desc_df['hbd'] <= 5) &\r\n    (desc_df['hba'] <= 10)\r\n)\r\nfiltered_df = df[druglike]\r\n\r\ndiverse_mols = dm.pick_diverse(\r\n    filtered_df['mol'].tolist(),\r\n    npick=100\r\n)\r\n\r\ndm.viz.to_image(\r\n    diverse_mols,\r\n    legends=[dm.to_smiles(m) for m in diverse_mols],\r\n    outfile=\"diverse_compounds.png\",\r\n    n_cols=10\r\n)\r\n```\r\n\r\n### Structure-Activity Relationship (SAR) Analysis\r\n\r\n```python\r\nscaffolds = [dm.to_scaffold_murcko(mol) for mol in mols]\r\nscaffold_smiles = [dm.to_smiles(s) for s in scaffolds]\r\n\r\nsar_df = pd.DataFrame({\r\n    'mol': mols,\r\n    'scaffold': scaffold_smiles,\r\n    'activity': activities  # User-provided activity data\r\n})\r\n\r\nfor scaffold, group in sar_df.groupby('scaffold'):\r\n    if len(group) >= 3:  # Need multiple examples\r\n        print(f\"\\nScaffold: {scaffold}\")\r\n        print(f\"Count: {len(group)}\")\r\n        print(f\"Activity range: {group['activity'].min():.2f} - {group['activity'].max():.2f}\")\r\n\r\n        # Visualize with activities as legends\r\n        dm.viz.to_image(\r\n            group['mol'].tolist(),\r\n            legends=[f\"Activity: {act:.2f}\" for act in group['activity']],\r\n            align=True  # Align by common substructure\r\n        )\r\n```\r\n\r\n### Virtual Screening Pipeline\r\n\r\n```python\r\nquery_fps = [dm.to_fp(mol) for mol in query_actives]\r\nlibrary_fps = [dm.to_fp(mol) for mol in library_mols]\r\n\r\nfrom scipy.spatial.distance import cdist\r\nimport numpy as np\r\n\r\ndistances = dm.cdist(query_actives, library_mols, n_jobs=-1)\r\n\r\nmin_distances = distances.min(axis=0)\r\nsimilarities = 1 - min_distances  # Convert distance to similarity\r\n\r\ntop_indices = np.argsort(similarities)[::-1][:100]  # Top 100\r\ntop_hits = [library_mols[i] for i in top_indices]\r\ntop_scores = [similarities[i] for i in top_indices]\r\n\r\n\r\n```python\r\ndef safe_to_mol(smiles):\r\n    try:\r\n        mol = dm.to_mol(smiles)\r\n        if mol is not None:\r\n            mol = dm.standardize_mol(mol)\r\n        return mol\r\n    except Exception as e:\r\n        print(f\"Failed to process {smiles}: {e}\")\r\n        return None\r\n\r\n\r\n```python\r\nX = np.array([dm.to_fp(mol) for mol in mols])\r\n\r\ndesc_df = dm.descriptors.batch_compute_many_descriptors(mols, n_jobs=-1)\r\nX = desc_df.values\r\n\r\nfrom sklearn.ensemble import RandomForestRegressor\r\nmodel = RandomForestRegressor()\r\nmodel.fit(X, y_target)",
  "name": "datamol",
  "id": "scientific-pkg-datamol",
  "sections": {
    "Reference Documentation": "For detailed API documentation, consult these reference files:\r\n\r\n- **`references/core_api.md`**: Core namespace functions (conversions, standardization, fingerprints, clustering)\r\n- **`references/io_module.md`**: File I/O operations (read/write SDF, CSV, Excel, remote files)\r\n- **`references/conformers_module.md`**: 3D conformer generation, clustering, SASA calculations\r\n- **`references/descriptors_viz.md`**: Molecular descriptors and visualization functions\r\n- **`references/fragments_scaffolds.md`**: Scaffold extraction, BRICS/RECAP fragmentation\r\n- **`references/reactions_data.md`**: Chemical reactions and toy datasets",
    "Additional Resources": "- **Datamol Documentation**: https://docs.datamol.io/\r\n- **RDKit Documentation**: https://www.rdkit.org/docs/\r\n- **GitHub Repository**: https://github.com/datamol-io/datamol",
    "Overview": "Datamol is a Python library that provides a lightweight, Pythonic abstraction layer over RDKit for molecular cheminformatics. Simplify complex molecular operations with sensible defaults, efficient parallelization, and modern I/O capabilities. All molecular objects are native `rdkit.Chem.Mol` instances, ensuring full compatibility with the RDKit ecosystem.\r\n\r\n**Key capabilities**:\r\n- Molecular format conversion (SMILES, SELFIES, InChI)\r\n- Structure standardization and sanitization\r\n- Molecular descriptors and fingerprints\r\n- 3D conformer generation and analysis\r\n- Clustering and diversity selection\r\n- Scaffold and fragment analysis\r\n- Chemical reaction application\r\n- Visualization and alignment\r\n- Batch processing with parallelization\r\n- Cloud storage support via fsspec",
    "Troubleshooting": "**Issue**: Molecule parsing fails\r\n- **Solution**: Use `dm.standardize_smiles()` first or try `dm.fix_mol()`\r\n\r\n**Issue**: Memory errors with clustering\r\n- **Solution**: Use `dm.pick_diverse()` instead of full clustering for large sets\r\n\r\n**Issue**: Slow conformer generation\r\n- **Solution**: Reduce `n_confs` or increase `rms_cutoff` to generate fewer conformers\r\n\r\n**Issue**: Remote file access fails\r\n- **Solution**: Ensure fsspec and appropriate cloud provider libraries are installed (s3fs, gcsfs, etc.)",
    "Common Workflows and Patterns": "dm.viz.to_image(\r\n    top_hits[:20],\r\n    legends=[f\"Sim: {score:.3f}\" for score in top_scores[:20]],\r\n    outfile=\"screening_hits.png\"\r\n)\r\n```",
    "Installation and Setup": "pip install datamol\r\n```\r\n\r\n**Import convention**:\r\n```python\r\nimport datamol as dm\r\n```",
    "Parallelization": "Datamol includes built-in parallelization for many operations. Use `n_jobs` parameter:\r\n- `n_jobs=1`: Sequential (no parallelization)\r\n- `n_jobs=-1`: Use all available CPU cores\r\n- `n_jobs=4`: Use 4 cores\r\n\r\n**Functions supporting parallelization**:\r\n- `dm.read_sdf(..., n_jobs=-1)`\r\n- `dm.descriptors.batch_compute_many_descriptors(..., n_jobs=-1)`\r\n- `dm.cluster_mols(..., n_jobs=-1)`\r\n- `dm.pdist(..., n_jobs=-1)`\r\n- `dm.conformers.sasa(..., n_jobs=-1)`\r\n\r\n**Progress bars**: Many batch operations support `progress=True` parameter.",
    "Core Workflows": "products = []\r\nfor mol in reactant_mols:\r\n    try:\r\n        prod = dm.reactions.apply_reaction(rxn, (mol,))\r\n        if prod is not None:\r\n            products.append(prod)\r\n    except Exception as e:\r\n        print(f\"Reaction failed: {e}\")\r\n```",
    "Integration with Machine Learning": "predictions = model.predict(X_test)\r\n```",
    "Error Handling": "valid_mols = []\r\nfor smiles in smiles_list:\r\n    mol = safe_to_mol(smiles)\r\n    if mol is not None:\r\n        valid_mols.append(mol)\r\n```",
    "Best Practices": "1. **Always standardize molecules** from external sources:\r\n   ```python\r\n   mol = dm.standardize_mol(mol, disconnect_metals=True, normalize=True, reionize=True)\r\n   ```\r\n\r\n2. **Check for None values** after molecule parsing:\r\n   ```python\r\n   mol = dm.to_mol(smiles)\r\n   if mol is None:\r\n       # Handle invalid SMILES\r\n   ```\r\n\r\n3. **Use parallel processing** for large datasets:\r\n   ```python\r\n   result = dm.operation(..., n_jobs=-1, progress=True)\r\n   ```\r\n\r\n4. **Leverage fsspec** for cloud storage:\r\n   ```python\r\n   df = dm.read_sdf(\"s3://bucket/compounds.sdf\")\r\n   ```\r\n\r\n5. **Use appropriate fingerprints** for similarity:\r\n   - ECFP (Morgan): General purpose, structural similarity\r\n   - MACCS: Fast, smaller feature space\r\n   - Atom pairs: Considers atom pairs and distances\r\n\r\n6. **Consider scale limitations**:\r\n   - Butina clustering: ~1,000 molecules (full distance matrix)\r\n   - For larger datasets: Use diversity selection or hierarchical methods\r\n\r\n7. **Scaffold splitting for ML**: Ensure proper train/test separation by scaffold\r\n\r\n8. **Align molecules** when visualizing SAR series"
  }
}

---
name: datamol
description: "Pythonic wrapper around RDKit with simplified interface and sensible defaults. Preferred for standard drug discovery: SMILES parsing, standardization, descriptors, fingerprints, clustering, 3D conformers, parallel processing. Returns native rdkit.Chem.Mol objects. For advanced control or custom parameters, use rdkit directly."
---

# Datamol Cheminformatics Skill

## Overview

Datamol is a Python library that provides a lightweight, Pythonic abstraction layer over RDKit for molecular cheminformatics. Simplify complex molecular operations with sensible defaults, efficient parallelization, and modern I/O capabilities. All molecular objects are native `rdkit.Chem.Mol` instances, ensuring full compatibility with the RDKit ecosystem.

**Key capabilities**:
- Molecular format conversion (SMILES, SELFIES, InChI)
- Structure standardization and sanitization
- Molecular descriptors and fingerprints
- 3D conformer generation and analysis
- Clustering and diversity selection
- Scaffold and fragment analysis
- Chemical reaction application
- Visualization and alignment
- Batch processing with parallelization
- Cloud storage support via fsspec

## Installation and Setup

Guide users to install datamol:

```bash
# Via conda/mamba (recommended)
conda install -c conda-forge datamol

# Via pip
pip install datamol
```

**Import convention**:
```python
import datamol as dm
```

## Core Workflows

### 1. Basic Molecule Handling

**Creating molecules from SMILES**:
```python
import datamol as dm

# Single molecule
mol = dm.to_mol("CCO")  # Ethanol

# From list of SMILES
smiles_list = ["CCO", "c1ccccc1", "CC(=O)O"]
mols = [dm.to_mol(smi) for smi in smiles_list]

# Error handling
mol = dm.to_mol("invalid_smiles")  # Returns None
if mol is None:
    print("Failed to parse SMILES")
```

**Converting molecules to SMILES**:
```python
# Canonical SMILES
smiles = dm.to_smiles(mol)

# Isomeric SMILES (includes stereochemistry)
smiles = dm.to_smiles(mol, isomeric=True)

# Other formats
inchi = dm.to_inchi(mol)
inchikey = dm.to_inchikey(mol)
selfies = dm.to_selfies(mol)
```

**Standardization and sanitization** (always recommend for user-provided molecules):
```python
# Sanitize molecule
mol = dm.sanitize_mol(mol)

# Full standardization (recommended for datasets)
mol = dm.standardize_mol(
    mol,
    disconnect_metals=True,
    normalize=True,
    reionize=True
)

# For SMILES strings directly
clean_smiles = dm.standardize_smiles(smiles)
```

### 2. Reading and Writing Molecular Files

Refer to `references/io_module.md` for comprehensive I/O documentation.

**Reading files**:
```python
# SDF files (most common in chemistry)
df = dm.read_sdf("compounds.sdf", mol_column='mol')

# SMILES files
df = dm.read_smi("molecules.smi", smiles_column='smiles', mol_column='mol')

# CSV with SMILES column
df = dm.read_csv("data.csv", smiles_column="SMILES", mol_column="mol")

# Excel files
df = dm.read_excel("compounds.xlsx", sheet_name=0, mol_column="mol")

# Universal reader (auto-detects format)
df = dm.open_df("file.sdf")  # Works with .sdf, .csv, .xlsx, .parquet, .json
```

**Writing files**:
```python
# Save as SDF
dm.to_sdf(mols, "output.sdf")
# Or from DataFrame
dm.to_sdf(df, "output.sdf", mol_column="mol")

# Save as SMILES file
dm.to_smi(mols, "output.smi")

# Excel with rendered molecule images
dm.to_xlsx(df, "output.xlsx", mol_columns=["mol"])
```

**Remote file support** (S3, GCS, HTTP):
```python
# Read from cloud storage
df = dm.read_sdf("s3://bucket/compounds.sdf")
df = dm.read_csv("https://example.com/data.csv")

# Write to cloud storage
dm.to_sdf(mols, "s3://bucket/output.sdf")
```

### 3. Molecular Descriptors and Properties

Refer to `references/descriptors_viz.md` for detailed descriptor documentation.

**Computing descriptors for a single molecule**:
```python
# Get standard descriptor set
descriptors = dm.descriptors.compute_many_descriptors(mol)
# Returns: {'mw': 46.07, 'logp': -0.03, 'hbd': 1, 'hba': 1,
#           'tpsa': 20.23, 'n_aromatic_atoms': 0, ...}
```

**Batch descriptor computation** (recommended for datasets):
```python
# Compute for all molecules in parallel
desc_df = dm.descriptors.batch_compute_many_descriptors(
    mols,
    n_jobs=-1,      # Use all CPU cores
    progress=True   # Show progress bar
)
```

**Specific descriptors**:
```python
# Aromaticity
n_aromatic = dm.descriptors.n_aromatic_atoms(mol)
aromatic_ratio = dm.descriptors.n_aromatic_atoms_proportion(mol)

# Stereochemistry
n_stereo = dm.descriptors.n_stereo_centers(mol)
n_unspec = dm.descriptors.n_stereo_centers_unspecified(mol)

# Flexibility
n_rigid = dm.descriptors.n_rigid_bonds(mol)
```

**Drug-likeness filtering (Lipinski's Rule of Five)**:
```python
# Filter compounds
def is_druglike(mol):
    desc = dm.descriptors.compute_many_descriptors(mol)
    return (
        desc['mw'] <= 500 and
        desc['logp'] <= 5 and
        desc['hbd'] <= 5 and
        desc['hba'] <= 10
    )

druglike_mols = [mol for mol in mols if is_druglike(mol)]
```

### 4. Molecular Fingerprints and Similarity

**Generating fingerprints**:
```python
# ECFP (Extended Connectivity Fingerprint, default)
fp = dm.to_fp(mol, fp_type='ecfp', radius=2, n_bits=2048)

# Other fingerprint types
fp_maccs = dm.to_fp(mol, fp_type='maccs')
fp_topological = dm.to_fp(mol, fp_type='topological')
fp_atompair = dm.to_fp(mol, fp_type='atompair')
```

**Similarity calculations**:
```python
# Pairwise distances within a set
distance_matrix = dm.pdist(mols, n_jobs=-1)

# Distances between two sets
distances = dm.cdist(query_mols, library_mols, n_jobs=-1)

# Find most similar molecules
from scipy.spatial.distance import squareform
dist_matrix = squareform(dm.pdist(mols))
# Lower distance = higher similarity (Tanimoto distance = 1 - Tanimoto similarity)
```

### 5. Clustering and Diversity Selection

Refer to `references/core_api.md` for clustering details.

**Butina clustering**:
```python
# Cluster molecules by structural similarity
clusters = dm.cluster_mols(
    mols,
    cutoff=0.2,    # Tanimoto distance threshold (0=identical, 1=completely different)
    n_jobs=-1      # Parallel processing
)

# Each cluster is a list of molecule indices
for i, cluster in enumerate(clusters):
    print(f"Cluster {i}: {len(cluster)} molecules")
    cluster_mols = [mols[idx] for idx in cluster]
```

**Important**: Butina clustering builds a full distance matrix - suitable for ~1000 molecules, not for 10,000+.

**Diversity selection**:
```python
# Pick diverse subset
diverse_mols = dm.pick_diverse(
    mols,
    npick=100  # Select 100 diverse molecules
)

# Pick cluster centroids
centroids = dm.pick_centroids(
    mols,
    npick=50   # Select 50 representative molecules
)
```

### 6. Scaffold Analysis

Refer to `references/fragments_scaffolds.md` for complete scaffold documentation.

**Extracting Murcko scaffolds**:
```python
# Get Bemis-Murcko scaffold (core structure)
scaffold = dm.to_scaffold_murcko(mol)
scaffold_smiles = dm.to_smiles(scaffold)
```

**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]]
```

### 7. Molecular Fragmentation

Refer to `references/fragments_scaffolds.md` for fragmentation details.

**BRICS fragmentation** (16 bond types):
```python
# Fragment molecule
fragments = dm.fragment.brics(mol)
# Returns: set of fragment SMILES with attachment points like '[1*]CCN'
```

**RECAP fragmentation** (11 bond types):
```python
fragments = dm.fragment.recap(mol)
```

**Fragment analysis**:
```python
# Find common fragments across compound library
from collections import Counter

all_fragments = []
for mol in mols:
    frags = dm.fragment.brics(mol)
    all_fragments.extend(frags)

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

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

### 8. 3D Conformer Generation

Refer to `references/conformers_module.md` for detailed conformer documentation.

**Generating conformers**:
```python
# Generate 3D conformers
mol_3d = dm.conformers.generate(
    mol,
    n_confs=50,           # Number to generate (auto if None)
    rms_cutoff=0.5,       # Filter similar conformers (Ångströms)
    minimize_energy=True,  # Minimize with UFF force field
    method='ETKDGv3'      # Embedding method (recommended)
)

# Access conformers
n_conformers = mol_3d.GetNumConformers()
conf = mol_3d.GetConformer(0)  # Get first conformer
positions = conf.GetPositions()  # Nx3 array of atom coordinates
```

**Conformer clustering**:
```python
# Cluster conformers by RMSD
clusters = dm.conformers.cluster(
    mol_3d,
    rms_cutoff=1.0,
    centroids=False
)

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

**SASA calculation**:
```python
# Calculate solvent accessible surface area
sasa_values = dm.conformers.sasa(mol_3d, n_jobs=-1)

# Access SASA from conformer properties
conf = mol_3d.GetConformer(0)
sasa = conf.GetDoubleProp('rdkit_free_sasa')
```

### 9. Visualization

Refer to `references/descriptors_viz.md` for visualization documentation.

**Basic molecule grid**:
```python
# Visualize molecules
dm.viz.to_image(
    mols[:20],
    legends=[dm.to_smiles(m) for m in mols[:20]],
    n_cols=5,
    mol_size=(300, 300)
)

# Save to file
dm.viz.to_image(mols, outfile="molecules.png")

# SVG for publications
dm.viz.to_image(mols, outfile="molecules.svg", use_svg=True)
```

**Aligned visualization** (for SAR analysis):
```python
# Align molecules by common substructure
dm.viz.to_image(
    similar_mols,
    align=True,  # Enable MCS alignment
    legends=activity_labels,
    n_cols=4
)
```

**Highlighting substructures**:
```python
# Highlight specific atoms and bonds
dm.viz.to_image(
    mol,
    highlight_atom=[0, 1, 2, 3],  # Atom indices
    highlight_bond=[0, 1, 2]      # Bond indices
)
```

**Conformer visualization**:
```python
# Display multiple conformers
dm.viz.conformers(
    mol_3d,
    n_confs=10,
    align_conf=True,
    n_cols=3
)
```

### 10. Chemical Reactions

Refer to `references/reactions_data.md` for reactions documentation.

**Applying reactions**:
```python
from rdkit.Chem import rdChemReactions

# Define reaction from SMARTS
rxn_smarts = '[C:1](=[O:2])[OH:3]>>[C:1](=[O:2])[Cl:3]'
rxn = rdChemReactions.ReactionFromSmarts(rxn_smarts)

# Apply to molecule
reactant = dm.to_mol("CC(=O)O")  # Acetic acid
product = dm.reactions.apply_reaction(
    rxn,
    (reactant,),
    sanitize=True
)

# Convert to SMILES
product_smiles = dm.to_smiles(product)
```

**Batch reaction application**:
```python
# Apply reaction to library
products = []
for mol in reactant_mols:
    try:
        prod = dm.reactions.apply_reaction(rxn, (mol,))
        if prod is not None:
            products.append(prod)
    except Exception as e:
        print(f"Reaction failed: {e}")
```

## Parallelization

Datamol includes built-in parallelization for many operations. Use `n_jobs` parameter:
- `n_jobs=1`: Sequential (no parallelization)
- `n_jobs=-1`: Use all available CPU cores
- `n_jobs=4`: Use 4 cores

**Functions supporting parallelization**:
- `dm.read_sdf(..., n_jobs=-1)`
- `dm.descriptors.batch_compute_many_descriptors(..., n_jobs=-1)`
- `dm.cluster_mols(..., n_jobs=-1)`
- `dm.pdist(..., n_jobs=-1)`
- `dm.conformers.sasa(..., n_jobs=-1)`

**Progress bars**: Many batch operations support `progress=True` parameter.

## Common Workflows and Patterns

### Complete Pipeline: Data Loading → Filtering → Analysis

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

# 1. Load molecules
df = dm.read_sdf("compounds.sdf")

# 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,
    progress=True
)

# 4. Filter by drug-likeness
druglike = (
    (desc_df['mw'] <= 500) &
    (desc_df['logp'] <= 5) &
    (desc_df['hbd'] <= 5) &
    (desc_df['hba'] <= 10)
)
filtered_df = df[druglike]

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

# 6. Visualize results
dm.viz.to_image(
    diverse_mols,
    legends=[dm.to_smiles(m) for m in diverse_mols],
    outfile="diverse_compounds.png",
    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. Generate fingerprints for query and library
query_fps = [dm.to_fp(mol) for mol in query_actives]
library_fps = [dm.to_fp(mol) for mol in library_mols]

# 2. Calculate similarities
from scipy.spatial.distance import cdist
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"
)
```

## Reference Documentation

For detailed API documentation, consult these reference files:

- **`references/core_api.md`**: Core namespace functions (conversions, standardization, fingerprints, clustering)
- **`references/io_module.md`**: File I/O operations (read/write SDF, CSV, Excel, remote files)
- **`references/conformers_module.md`**: 3D conformer generation, clustering, SASA calculations
- **`references/descriptors_viz.md`**: Molecular descriptors and visualization functions
- **`references/fragments_scaffolds.md`**: Scaffold extraction, BRICS/RECAP fragmentation
- **`references/reactions_data.md`**: Chemical reactions and toy datasets

## Best Practices

1. **Always standardize molecules** from external sources:
   ```python
   mol = dm.standardize_mol(mol, disconnect_metals=True, normalize=True, reionize=True)
   ```

2. **Check for None values** after molecule parsing:
   ```python
   mol = dm.to_mol(smiles)
   if mol is None:
       # Handle invalid SMILES
   ```

3. **Use parallel processing** for large datasets:
   ```python
   result = dm.operation(..., n_jobs=-1, progress=True)
   ```

4. **Leverage fsspec** for cloud storage:
   ```python
   df = dm.read_sdf("s3://bucket/compounds.sdf")
   ```

5. **Use appropriate fingerprints** for similarity:
   - ECFP (Morgan): General purpose, structural similarity
   - MACCS: Fast, smaller feature space
   - Atom pairs: Considers atom pairs and distances

6. **Consider scale limitations**:
   - Butina clustering: ~1,000 molecules (full distance matrix)
   - For larger datasets: Use diversity selection or hierarchical methods

7. **Scaffold splitting for ML**: Ensure proper train/test separation by scaffold

8. **Align molecules** when visualizing SAR series

## Error Handling

```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)
```

## 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)
X = desc_df.values

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

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

## Troubleshooting

**Issue**: Molecule parsing fails
- **Solution**: Use `dm.standardize_smiles()` first or try `dm.fix_mol()`

**Issue**: Memory errors with clustering
- **Solution**: Use `dm.pick_diverse()` instead of full clustering for large sets

**Issue**: Slow conformer generation
- **Solution**: Reduce `n_confs` or increase `rms_cutoff` to generate fewer conformers

**Issue**: Remote file access fails
- **Solution**: Ensure fsspec and appropriate cloud provider libraries are installed (s3fs, gcsfs, etc.)

## Additional Resources

- **Datamol Documentation**: https://docs.datamol.io/
- **RDKit Documentation**: https://www.rdkit.org/docs/
- **GitHub Repository**: https://github.com/datamol-io/datamol