Virtual Screening Library Preparation

This notebook demonstrates how to prepare compound libraries for virtual screening campaigns using Auto3D.

CLI Quick Start (Recommended)

For most virtual screening workflows, the CLI is the fastest approach:

# Basic virtual screening - generate 5 conformers per molecule
auto3d run compounds.smi --k=5 --gpu

# With energy window selection
auto3d run compounds.smi --k=5 --window=3.0 --gpu

# Fast screening with ANI2xt
auto3d run compounds.smi --k=3 --engine=ANI2xt --gpu

# Large library with multi-GPU
auto3d run large_library.smi --k=5 --gpu --gpu-idx="0,1,2,3"

# Validate input before processing
auto3d validate compounds.smi

# Use configuration preset for screening
auto3d config init -p quick -o screening.yaml
auto3d run compounds.smi -c screening.yaml --gpu

Example screening.yaml:

k: 5
window: 3.0
optimizing_engine: ANI2xt
threshold: 0.3
max_confs: 100
use_gpu: true

Python Workflow

This notebook covers:

  1. Library preprocessing - SMILES validation, standardization, filtering

  2. 3D conformer generation - optimal settings for docking

  3. Property-based filtering - Lipinski’s Rule of Five, lead-likeness

  4. Output formats - docking-ready SDF files with proper metadata

Chemistry Background

Virtual screening requires high-quality 3D conformers that represent the bioactive conformation. Key considerations:

  • Conformational sampling: Drug-like molecules typically have 1-10 rotatable bonds, generating 10s-1000s of conformers

  • Energy window: Bioactive conformations are usually within 3-5 kcal/mol of the global minimum

  • Protonation state: Most docking programs expect neutral or physiological pH forms

[ ]:

import os
import tempfile
from pathlib import Path

import Auto3D
from Auto3D import Auto3DOptions, main
from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski, AllChem, PandasTools
from rdkit.Chem.FilterCatalog import FilterCatalog, FilterCatalogParams
import pandas as pd

print(f"Auto3D version: {Auto3D.__version__}")

1. Library Preprocessing

Before 3D generation, we validate and standardize the input SMILES. This catches common issues:

  • Invalid SMILES syntax

  • Disconnected fragments (salts, counterions)

  • Non-drug-like molecules (too large, reactive groups)

[ ]:

# Example compound library - common drug scaffolds
# In practice, this would come from vendor catalogs or in-house databases
compound_library = {
    "aspirin": "CC(=O)OC1=CC=CC=C1C(=O)O",
    "ibuprofen": "CC(C)CC1=CC=C(C=C1)C(C)C(=O)O",
    "acetaminophen": "CC(=O)NC1=CC=C(C=C1)O",
    "caffeine": "CN1C=NC2=C1C(=O)N(C(=O)N2C)C",
    "naproxen": "COC1=CC2=CC(C(C)C(=O)O)=CC=C2C=C1",
    "celecoxib": "CC1=CC=C(C=C1)C2=CC(=NN2C3=CC=C(C=C3)S(=O)(=O)N)C(F)(F)F",
    "omeprazole": "COC1=CC2=NC(CS(=O)C3=NC4=CC=CC=C4N3)=NC2=C(OC)C1",
    "atorvastatin": "CC(C)C1=C(C(=C(N1CCC(CC(CC(=O)O)O)O)C2=CC=C(C=C2)F)C3=CC=CC=C3)C(=O)NC4=CC=CC=C4",
    # Invalid/problematic entries for demonstration
    "invalid_smiles": "not_a_smiles",
    "salt_form": "CC(=O)OC1=CC=CC=C1C(=O)[O-].[Na+]",  # Sodium salt
    "too_large": "C" * 100,  # Very long alkane
}

[ ]:

def validate_and_standardize(smiles_dict, max_mw=600, max_rotatable=12):
    """
    Validate and standardize a SMILES library for virtual screening.

    Filters applied:
    - Valid SMILES parsing
    - Remove salts/counterions (keep largest fragment)
    - Molecular weight cutoff
    - Rotatable bond limit (conformational complexity)

    Returns:
        dict: Validated {name: canonical_smiles}
        list: Rejected compounds with reasons
    """
    valid = {}
    rejected = []

    for name, smi in smiles_dict.items():
        # Try to parse
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            rejected.append((name, "Invalid SMILES"))
            continue

        # Handle salts - keep largest fragment
        fragments = Chem.GetMolFrags(mol, asMols=True, sanitizeFrags=True)
        if len(fragments) > 1:
            # Keep the largest fragment by heavy atom count
            mol = max(fragments, key=lambda m: m.GetNumHeavyAtoms())

        # Check molecular weight
        mw = Descriptors.MolWt(mol)
        if mw > max_mw:
            rejected.append((name, f"MW {mw:.1f} > {max_mw}"))
            continue

        # Check rotatable bonds (conformational complexity)
        n_rot = Lipinski.NumRotatableBonds(mol)
        if n_rot > max_rotatable:
            rejected.append((name, f"Rotatable bonds {n_rot} > {max_rotatable}"))
            continue

        # Store canonical SMILES
        valid[name] = Chem.MolToSmiles(mol, canonical=True)

    return valid, rejected


valid_compounds, rejected = validate_and_standardize(compound_library)

print(f"Valid compounds: {len(valid_compounds)}")
print(f"Rejected: {len(rejected)}")
print("\nRejection reasons:")
for name, reason in rejected:
    print(f"  {name}: {reason}")

2. Lipinski’s Rule of Five Filtering

For oral drug candidates, Lipinski’s Rule of Five provides a useful filter:

  • MW ≤ 500 - Molecular weight

  • LogP ≤ 5 - Lipophilicity

  • HBD ≤ 5 - Hydrogen bond donors

  • HBA ≤ 10 - Hydrogen bond acceptors

Compounds violating ≥2 rules have poor oral bioavailability.

[ ]:

def calculate_lipinski_properties(smiles_dict):
    """
    Calculate Lipinski Rule of Five properties.

    Returns DataFrame with:
    - MW, LogP, HBD, HBA
    - Number of violations
    - Pass/Fail status
    """
    data = []

    for name, smi in smiles_dict.items():
        mol = Chem.MolFromSmiles(smi)

        mw = Descriptors.MolWt(mol)
        logp = Descriptors.MolLogP(mol)
        hbd = Lipinski.NumHDonors(mol)
        hba = Lipinski.NumHAcceptors(mol)

        # Count violations
        violations = sum([
            mw > 500,
            logp > 5,
            hbd > 5,
            hba > 10
        ])

        data.append({
            "Name": name,
            "SMILES": smi,
            "MW": round(mw, 1),
            "LogP": round(logp, 2),
            "HBD": hbd,
            "HBA": hba,
            "Violations": violations,
            "Ro5_Pass": violations < 2
        })

    return pd.DataFrame(data)


df_props = calculate_lipinski_properties(valid_compounds)
print(df_props.to_string(index=False))

[ ]:

# Filter to Ro5-compliant compounds
df_druglike = df_props[df_props["Ro5_Pass"]].copy()
print(f"Drug-like compounds (Ro5 pass): {len(df_druglike)}")

3. PAINS Filtering

Pan-Assay Interference Compounds (PAINS) are molecules that give false positives in many assays due to:

  • Aggregation

  • Redox cycling

  • Fluorescence interference

  • Covalent modification

RDKit includes PAINS filters that should be applied before screening.

[ ]:

def filter_pains(smiles_dict):
    """
    Filter out PAINS (Pan-Assay Interference Compounds).

    These compounds cause false positives in biochemical assays.
    """
    # Initialize PAINS filter catalog
    params = FilterCatalogParams()
    params.AddCatalog(FilterCatalogParams.FilterCatalogs.PAINS)
    catalog = FilterCatalog(params)

    clean = {}
    pains_hits = []

    for name, smi in smiles_dict.items():
        mol = Chem.MolFromSmiles(smi)
        entry = catalog.GetFirstMatch(mol)

        if entry is None:
            clean[name] = smi
        else:
            pains_hits.append((name, entry.GetDescription()))

    return clean, pains_hits


# Apply PAINS filter
druglike_smiles = dict(zip(df_druglike["Name"], df_druglike["SMILES"]))
clean_compounds, pains = filter_pains(druglike_smiles)

print(f"Compounds passing PAINS filter: {len(clean_compounds)}")
if pains:
    print(f"\nPAINS hits removed:")
    for name, desc in pains:
        print(f"  {name}: {desc}")

4. 3D Conformer Generation for Docking

For virtual screening, we need:

  • Multiple conformers per molecule (k=3-10) to sample conformational space

  • Energy window (typically 3-5 kcal/mol) to include bioactive conformations

  • Reasonable optimization - not too tight (slow) but well-minimized

Auto3D’s neural network potentials provide accurate geometries much faster than force fields.

CLI (Recommended)

# Basic virtual screening setup
auto3d run filtered_compounds.smi --k=5 --window=3.0 --gpu

# Fast screening with ANI2xt
auto3d run filtered_compounds.smi --k=3 --engine=ANI2xt --gpu

# Multi-GPU for large libraries
auto3d run library.smi --k=5 --gpu --gpu-idx="0,1,2,3"

# With configuration file
auto3d config init -p quick -o screening.yaml
auto3d run filtered_compounds.smi -c screening.yaml

Example screening.yaml for virtual screening:

k: 5
window: 3.0
optimizing_engine: ANI2xt
threshold: 0.3
max_confs: 100
use_gpu: true

Python API

Below is an example using the Python API for programmatic workflows:

[ ]:

# Configure Auto3D for docking preparation
# Key settings for virtual screening:
# - k=5: Generate up to 5 diverse conformers per molecule
# - window=3.0: Include conformers within 3 kcal/mol of minimum
# - threshold=0.3: RMSD threshold for conformer diversity

if __name__ == "__main__":
    config = Auto3DOptions(
        path=input_file,
        k=5,                        # Top 5 conformers
        window=3.0,                 # Within 3 kcal/mol of minimum
        optimizing_engine="AIMNET", # Fast and accurate
        threshold=0.3,              # RMSD diversity threshold (Angstroms)
        max_confs=100,              # Initial conformer pool
        use_gpu=True,
        verbose=True,
    )

    output_sdf = main(config)
    print(f"\nOutput: {output_sdf}")

5. Post-processing for Docking

Most docking programs require specific preparation:

  • AutoDock Vina: PDBQT format with Gasteiger charges

  • Glide: Maestro format with proper protonation

  • GOLD: MOL2 format with proper atom types

Here we add properties useful for post-docking analysis.

[ ]:

def add_docking_properties(sdf_path, output_path=None):
    """
    Add properties useful for docking analysis:
    - Formal charge
    - Number of rotatable bonds
    - TPSA (topological polar surface area)
    - Aromatic rings
    """
    if output_path is None:
        output_path = sdf_path.replace(".sdf", "_docking_ready.sdf")

    mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))

    with Chem.SDWriter(output_path) as writer:
        for mol in mols:
            if mol is None:
                continue

            # Calculate useful properties
            mol.SetProp("FormalCharge", str(Chem.GetFormalCharge(mol)))
            mol.SetProp("NumRotatableBonds", str(Lipinski.NumRotatableBonds(mol)))
            mol.SetProp("TPSA", f"{Descriptors.TPSA(mol):.1f}")
            mol.SetProp("NumAromaticRings", str(Descriptors.NumAromaticRings(mol)))
            mol.SetProp("NumHeavyAtoms", str(mol.GetNumHeavyAtoms()))

            writer.write(mol)

    print(f"Wrote {len(mols)} conformers to {output_path}")
    return output_path


# Add docking-relevant properties
if 'output_sdf' in dir():
    docking_ready = add_docking_properties(output_sdf)
    print(f"Docking-ready file: {docking_ready}")

6. Conformer Analysis

Before docking, inspect the generated conformers to ensure quality.

[ ]:

def analyze_conformers(sdf_path):
    """
    Analyze the generated conformer ensemble.
    """
    mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))

    # Group by molecule name
    from collections import defaultdict
    conformers_by_mol = defaultdict(list)

    for mol in mols:
        if mol is None:
            continue
        name = mol.GetProp("_Name").split("_")[0]  # Remove conformer suffix
        energy = float(mol.GetProp("E_rel")) if mol.HasProp("E_rel") else 0.0
        conformers_by_mol[name].append(energy)

    print("Conformer Summary:")
    print("-" * 50)
    print(f"{'Molecule':<20} {'Count':>6} {'E_range (kcal/mol)':>20}")
    print("-" * 50)

    for name, energies in conformers_by_mol.items():
        e_range = max(energies) - min(energies) if len(energies) > 1 else 0.0
        print(f"{name:<20} {len(energies):>6} {e_range:>20.2f}")

    print("-" * 50)
    print(f"Total conformers: {len(mols)}")


if 'output_sdf' in dir():
    analyze_conformers(output_sdf)

Summary

This workflow prepared a compound library for virtual screening:

  1. Validated SMILES and removed problematic entries

  2. Filtered by Lipinski’s Rule of Five for oral bioavailability

  3. Removed PAINS compounds that cause assay interference

  4. Generated diverse 3D conformers with Auto3D

  5. Annotated with docking-relevant properties

The output SDF file is ready for docking with programs like:

  • AutoDock Vina (convert to PDBQT with prepare_ligand4.py)

  • Glide (import into Maestro)

  • GOLD (convert to MOL2)

[ ]:

# Cleanup
if 'input_file' in dir():
    os.unlink(input_file)