Rdkit

RDKit - Cheminformatics and Drug Discovery

RDKit is the industry-standard open-source toolkit for cheminformatics. It provides comprehensive tools for molecular manipulation, descriptor calculation, fingerprinting, substructure searching, and 3D molecular modeling. RDKit is used extensively in pharmaceutical companies for drug discovery and virtual screening.

When to Use

• Reading and writing chemical file formats (SMILES, SDF, MOL2, PDB).
• Calculating molecular descriptors and drug-like properties (Lipinski's Rule of Five).
• Generating molecular fingerprints for similarity searching.
• Substructure searching and chemical pattern matching (SMARTS).
• 3D conformer generation and molecular alignment.
• Virtual screening of compound libraries.
• Pharmacophore modeling and shape similarity.
• QSAR (Quantitative Structure-Activity Relationship) modeling.
• Reaction enumeration and retrosynthesis.
• Visualizing chemical structures in 2D and 3D.
• Building machine learning models for molecular property prediction.

Reference Documentation

Official docs: https://www.rdkit.org/docs/ RDKit Book: https://www.rdkit.org/docs/RDKit_Book.html GitHub: https://github.com/rdkit/rdkit Search patterns: rdkit.Chem, rdkit.Chem.Descriptors, rdkit.Chem.AllChem, rdkit.DataStructs

Core Principles

Molecular Representation

RDKit represents molecules as graphs where atoms are nodes and bonds are edges. The core object is Mol, which can be created from SMILES, SDF files, or built programmatically.

SMILES (Simplified Molecular Input Line Entry System)

A text-based notation for chemical structures. Example: CCO is ethanol, c1ccccc1 is benzene. RDKit can parse and generate SMILES strings.

Fingerprints for Similarity

Molecular fingerprints are binary vectors encoding structural features. They enable fast similarity searching and clustering of large compound libraries.

Lazy Evaluation

Many RDKit operations are lazy - properties are computed only when needed. This makes operations on large libraries very efficient.

Quick Reference

Installation

# Via conda (recommended)
conda install -c conda-forge rdkit

# Via pip
pip install rdkit

# For visualization
pip install rdkit pillow

Standard Imports

from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors, Draw, Lipinski
from rdkit.Chem import rdFingerprintGenerator
from rdkit import DataStructs
import numpy as np
import pandas as pd

Basic Pattern - SMILES to Molecule

from rdkit import Chem

# 1. Create molecule from SMILES
smiles = "CC(=O)OC1=CC=CC=C1C(=O)O"  # Aspirin
mol = Chem.MolFromSmiles(smiles)

# 2. Check if molecule is valid
if mol is None:
    print("Invalid SMILES")
else:
    print(f"Molecular formula: {Chem.rdMolDescriptors.CalcMolFormula(mol)}")
    print(f"Molecular weight: {Descriptors.MolWt(mol):.2f}")

Basic Pattern - Calculate Properties

from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski

mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O")

# Calculate drug-like properties
mw = Descriptors.MolWt(mol)
logp = Descriptors.MolLogP(mol)
hbd = Lipinski.NumHDonors(mol)
hba = Lipinski.NumHAcceptors(mol)

print(f"MW: {mw:.2f}, LogP: {logp:.2f}, HBD: {hbd}, HBA: {hba}")

# Check Lipinski's Rule of Five
lipinski_pass = (mw >[C:1](=O)OC')

# Apply reaction
products = rxn.RunReactants((mol,))

if products:
    product = products[0][0]
    print(f"Product: {Chem.MolToSmiles(product)}")  # Methyl acetate

3D Conformer Generation

Generate 3D Coordinates

from rdkit import Chem
from rdkit.Chem import AllChem

# Create molecule
mol = Chem.MolFromSmiles("CCO")

# Add hydrogens (required for 3D)
mol = Chem.AddHs(mol)

# Generate 3D coordinates
result = AllChem.EmbedMolecule(mol, randomSeed=42)

if result == 0:  # Success
    print("3D coordinates generated")
else:
    print("Failed to generate 3D coordinates")

# Optimize geometry with MMFF force field
AllChem.MMFFOptimizeMolecule(mol)

# Get atomic positions
conf = mol.GetConformer()
for i in range(mol.GetNumAtoms()):
    pos = conf.GetAtomPosition(i)
    print(f"Atom {i}: ({pos.x:.3f}, {pos.y:.3f}, {pos.z:.3f})")

Multiple Conformers

from rdkit import Chem
from rdkit.Chem import AllChem

mol = Chem.MolFromSmiles("CCCC")  # Butane
mol = Chem.AddHs(mol)

# Generate multiple conformers
conf_ids = AllChem.EmbedMultipleConfs(
    mol,
    numConfs=10,
    randomSeed=42,
    pruneRmsThresh=0.5  # Remove similar conformers
)

print(f"Generated {len(conf_ids)} conformers")

# Optimize each conformer
for conf_id in conf_ids:
    AllChem.MMFFOptimizeMolecule(mol, confId=conf_id)

# Get energies
props = AllChem.MMFFGetMoleculeProperties(mol)
for conf_id in conf_ids:
    ff = AllChem.MMFFGetMoleculeForceField(mol, props, confId=conf_id)
    energy = ff.CalcEnergy()
    print(f"Conformer {conf_id}: {energy:.2f} kcal/mol")

Molecular Alignment

from rdkit import Chem
from rdkit.Chem import AllChem

# Reference molecule
ref_mol = Chem.MolFromSmiles("c1ccccc1C")  # Toluene
ref_mol = Chem.AddHs(ref_mol)
AllChem.EmbedMolecule(ref_mol)

# Probe molecule
probe_mol = Chem.MolFromSmiles("c1ccccc1CC")  # Ethylbenzene
probe_mol = Chem.AddHs(probe_mol)
AllChem.EmbedMolecule(probe_mol)

# Align probe to reference
rmsd = AllChem.AlignMol(probe_mol, ref_mol)
print(f"RMSD: {rmsd:.3f} Å")

# Get aligned coordinates
# Now probe_mol has coordinates aligned to ref_mol

Molecular Visualization

2D Drawings

from rdkit import Chem
from rdkit.Chem import Draw
import matplotlib.pyplot as plt

# Single molecule
mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O")  # Aspirin
img = Draw.MolToImage(mol, size=(300, 300))
plt.imshow(img)
plt.axis('off')
plt.show()

# Multiple molecules
mols = [Chem.MolFromSmiles(s) for s in ["CCO", "c1ccccc1", "CC(=O)O"]]
legends = ["Ethanol", "Benzene", "Acetic acid"]

img = Draw.MolsToGridImage(
    mols,
    molsPerRow=3,
    subImgSize=(200, 200),
    legends=legends
)
plt.imshow(img)
plt.axis('off')
plt.show()

Highlight Substructures

from rdkit import Chem
from rdkit.Chem import Draw

mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O")  # Aspirin

# Highlight carboxylic acid group
pattern = Chem.MolFromSmarts("C(=O)O")
match = mol.GetSubstructMatch(pattern)

# Draw with highlighted atoms
img = Draw.MolToImage(mol, highlightAtoms=match, size=(300, 300))

Save to File

from rdkit import Chem
from rdkit.Chem import Draw

mol = Chem.MolFromSmiles("CCO")

# Save as PNG
Draw.MolToFile(mol, "molecule.png", size=(300, 300))

# Save as SVG (vector graphics)
from rdkit.Chem.Draw import rdMolDraw2D

drawer = rdMolDraw2D.MolDraw2DSVG(300, 300)
drawer.DrawMolecule(mol)
drawer.FinishDrawing()
svg = drawer.GetDrawingText()

with open("molecule.svg", "w") as f:
    f.write(svg)

Practical Workflows

1. Virtual Screening Pipeline

from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors, Lipinski
from rdkit import DataStructs
import pandas as pd

def screen_library(library_file, reference_smiles, similarity_threshold=0.7):
    """Screen compound library for similar, drug-like molecules."""
    
    # Reference molecule and fingerprint
    ref_mol = Chem.MolFromSmiles(reference_smiles)
    ref_fp = AllChem.GetMorganFingerprintAsBitVect(ref_mol, radius=2)
    
    # Results
    hits = []
    
    # Read library
    suppl = Chem.SDMolSupplier(library_file)
    
    for i, mol in enumerate(suppl):
        if mol is None:
            continue
        
        # Step 1: Drug-likeness filter (Lipinski)
        mw = Descriptors.MolWt(mol)
        logp = Descriptors.MolLogP(mol)
        hbd = Lipinski.NumHDonors(mol)
        hba = Lipinski.NumHAcceptors(mol)
        
        if not (mw >[C:1](=[O:2])[N:3]')
    
    products = []
    
    for acid_smi in acids:
        for amine_smi in amines:
            acid = Chem.MolFromSmiles(acid_smi)
            amine = Chem.MolFromSmiles(amine_smi)
            
            if acid is None or amine is None:
                continue
            
            # Run reaction
            products_tuple = rxn.RunReactants((acid, amine))
            
            if products_tuple:
                product = products_tuple[0][0]
                Chem.SanitizeMol(product)
                product_smi = Chem.MolToSmiles(product)
                products.append({
                    'acid': acid_smi,
                    'amine': amine_smi,
                    'product': product_smi
                })
    
    return products

# Usage
acids = ["CC(=O)O", "c1ccccc1C(=O)O"]
amines = ["CCN", "c1ccccc1N"]
products = enumerate_amide_coupling(acids, amines)
print(f"Generated {len(products)} amides")

Performance Optimization

Bulk Operations

from rdkit import Chem
from rdkit.Chem import AllChem
import pandas as pd

# Instead of loop
# ❌ SLOW
# fps = []
# for smi in smiles_list:
#     mol = Chem.MolFromSmiles(smi)
#     fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
#     fps.append(fp)

# ✅ FAST: Use PandasTools for bulk operations
from rdkit.Chem import PandasTools

df = pd.DataFrame({'SMILES': smiles_list})
PandasTools.AddMoleculeColumnToFrame(df, 'SMILES', 'Molecule')

# Calculate properties in bulk
df['MW'] = df['Molecule'].apply(lambda x: Descriptors.MolWt(x) if x else None)
df['LogP'] = df['Molecule'].apply(lambda x: Descriptors.MolLogP(x) if x else None)

Parallel Processing

from rdkit import Chem
from rdkit.Chem import AllChem
from multiprocessing import Pool
import pandas as pd

def process_molecule(smiles):
    """Process single molecule."""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None
    
    return {
        'smiles': smiles,
        'mw': Descriptors.MolWt(mol),
        'logp': Descriptors.MolLogP(mol),
        'fp': AllChem.GetMorganFingerprintAsBitVect(mol, 2)
    }

def process_library_parallel(smiles_list, n_jobs=4):
    """Process library in parallel."""
    with Pool(n_jobs) as pool:
        results = pool.map(process_molecule, smiles_list)
    
    # Filter None results
    results = [r for r in results if r is not None]
    return pd.DataFrame(results)

# Usage
# df = process_library_parallel(large_smiles_list, n_jobs=8)

Caching Calculations

from functools import lru_cache
from rdkit import Chem
from rdkit.Chem import Descriptors

@lru_cache(maxsize=10000)
def get_mol_properties(smiles):
    """Calculate properties with caching."""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None
    
    return {
        'mw': Descriptors.MolWt(mol),
        'logp': Descriptors.MolLogP(mol),
        'tpsa': Descriptors.TPSA(mol)
    }

# Repeated calls will use cache
props1 = get_mol_properties("CCO")  # Calculated
props2 = get_mol_properties("CCO")  # Cached (instant)

Common Pitfalls and Solutions

The "Invalid SMILES" Problem

Not all SMILES strings are valid.

# ❌ Problem: Assuming all SMILES are valid
smiles_list = ["CCO", "INVALID", "c1ccccc1"]
mols = [Chem.MolFromSmiles(s) for s in smiles_list]
# Contains None!

# ✅ Solution: Filter invalid molecules
mols = [Chem.MolFromSmiles(s) for s in smiles_list]
valid_mols = [m for m in mols if m is not None]

# ✅ Better: Track which failed
results = []
for smi in smiles_list:
    mol = Chem.MolFromSmiles(smi)
    if mol is not None:
        results.append({'smiles': smi, 'mol': mol, 'valid': True})
    else:
        results.append({'smiles': smi, 'mol': None, 'valid': False})

The "Stereochemistry Loss" Problem

SMILES generation can lose stereochemistry if not careful.

from rdkit import Chem

# Molecule with stereochemistry
chiral_smiles = "C[C@H](O)CC"  # (S)-2-butanol
mol = Chem.MolFromSmiles(chiral_smiles)

# ❌ BAD: Lose stereochemistry
non_iso = Chem.MolToSmiles(mol, isomericSmiles=False)
print(non_iso)  # "CC(O)CC" - lost chirality!

# ✅ GOOD: Preserve stereochemistry
iso = Chem.MolToSmiles(mol, isomericSmiles=True)
print(iso)  # "C[C@H](O)CC" - preserved!

The "3D Without Hydrogens" Problem

3D operations require explicit hydrogens.

from rdkit import Chem
from rdkit.Chem import AllChem

mol = Chem.MolFromSmiles("CCO")

# ❌ BAD: Generate 3D without hydrogens
result = AllChem.EmbedMolecule(mol)
# Poor quality or failure

# ✅ GOOD: Add hydrogens first
mol_h = Chem.AddHs(mol)
result = AllChem.EmbedMolecule(mol_h)
AllChem.MMFFOptimizeMolecule(mol_h)

The "Fingerprint Type Mismatch" Problem

Comparing different fingerprint types gives meaningless results.

from rdkit import Chem
from rdkit.Chem import AllChem, MACCSkeys
from rdkit import DataStructs

mol1 = Chem.MolFromSmiles("CCO")
mol2 = Chem.MolFromSmiles("CCCO")

# ❌ BAD: Comparing different fingerprint types
fp1_morgan = AllChem.GetMorganFingerprintAsBitVect(mol1, 2)
fp2_maccs = MACCSkeys.GenMACCSKeys(mol2)

# This will error or give nonsense!
# similarity = DataStructs.TanimotoSimilarity(fp1_morgan, fp2_maccs)

# ✅ GOOD: Use same fingerprint type
fp1 = AllChem.GetMorganFingerprintAsBitVect(mol1, 2)
fp2 = AllChem.GetMorganFingerprintAsBitVect(mol2, 2)
similarity = DataStructs.TanimotoSimilarity(fp1, fp2)

The "Memory Explosion" Problem

Processing millions of molecules can exhaust memory.

# ❌ BAD: Load entire library into memory
suppl = Chem.SDMolSupplier('huge_library.sdf')
mols = [mol for mol in suppl]  # Out of memory!

# ✅ GOOD: Process in batches
def process_in_batches(sdf_file, batch_size=10000):
    suppl = Chem.SDMolSupplier(sdf_file)
    batch = []
    
    for mol in suppl:
        if mol is not None:
            batch.append(mol)
        
        if len(batch) >= batch_size:
            # Process batch
            yield batch
            batch = []
    
    # Process remaining
    if batch:
        yield batch

# Usage
for batch in process_in_batches('huge_library.sdf'):
    # Process each batch
    pass

RDKit is the cornerstone of computational drug discovery and cheminformatics. Its comprehensive toolkit for molecular manipulation, descriptor calculation, and similarity searching makes it indispensable for pharmaceutical research, virtual screening, and chemical data analysis. Master RDKit, and you'll have the power to computationally explore vast chemical spaces and accelerate drug discovery.

Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

• Author: tondevrel
• Source: tondevrel/scientific-agent-skills
• License: MIT

Install and usage instructions live in the source repository linked above.