3D Molecular Descriptors for ML/QSAR

This notebook demonstrates calculating 3D molecular descriptors from Auto3D-generated conformers for:

• QSAR modeling - Quantitative Structure-Activity Relationships

• Machine learning - Features for property prediction

• Molecular similarity - 3D pharmacophore fingerprints

• Drug design - Shape and property optimization

Descriptor Categories

Type	Description	Examples
Constitutional	Atom/bond counts	MW, nAtoms, nRotBonds
Topological	2D connectivity	TPSA, Chi indices
Geometrical	3D shape/size	RadiusOfGyration, Asphericity
Electronic	Charge distribution	Dipole, partial charges
Quantum	From NNP energies	HOMO-LUMO gap, atomization energy

The key advantage of 3D descriptors is capturing conformationally-dependent properties.

import os
import tempfile
from pathlib import Path

import numpy as np
import Auto3D
from Auto3D import Auto3DOptions, main
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors, Descriptors3D
from rdkit.Chem import rdMolDescriptors, Lipinski
from rdkit.Chem import rdFreeSASA
from rdkit.ML.Descriptors import MoleculeDescriptors
import pandas as pd

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

1. Generate 3D Structures

First, generate optimized 3D structures using Auto3D.

# Drug molecules for descriptor calculation
drug_dataset = {
    "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",
    "diazepam": "CN1C(=O)CN=C(C2=CC=CC=C2)C3=CC=CC(Cl)=C31",
    "fluoxetine": "CNCCC(C1=CC=CC=C1)OC2=CC=C(C=C2)C(F)(F)F",
    "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",
    "metformin": "CN(C)C(=N)NC(=N)N",
    "omeprazole": "COC1=CC2=NC(CS(=O)C3=NC4=CC=CC=C4N3)=NC2=CC1OC",
}

with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:
    for name, smi in drug_dataset.items():
        f.write(f"{smi} {name}\n")
    input_file = f.name

print(f"Dataset: {len(drug_dataset)} molecules")

# Generate 3D structures
if __name__ == "__main__":
    config = Auto3DOptions(
        path=input_file,
        k=1,                        # Lowest energy conformer
        optimizing_engine="AIMNET",
        use_gpu=True,
    )

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

2. 2D Descriptors (Constitutional & Topological)

These are calculated from the molecular graph (don’t require 3D).

def calculate_2d_descriptors(mol):
    """
    Calculate key 2D molecular descriptors.
    """
    return {
        # Constitutional
        "MW": Descriptors.MolWt(mol),
        "nHeavyAtoms": mol.GetNumHeavyAtoms(),
        "nRotBonds": Lipinski.NumRotatableBonds(mol),
        "nRings": rdMolDescriptors.CalcNumRings(mol),
        "nAromaticRings": rdMolDescriptors.CalcNumAromaticRings(mol),
        "nHeteroatoms": rdMolDescriptors.CalcNumHeteroatoms(mol),

        # Lipinski
        "HBD": Lipinski.NumHDonors(mol),
        "HBA": Lipinski.NumHAcceptors(mol),
        "LogP": Descriptors.MolLogP(mol),
        "TPSA": Descriptors.TPSA(mol),

        # Electronic
        "FormalCharge": Chem.GetFormalCharge(mol),

        # Complexity
        "FractionCSP3": rdMolDescriptors.CalcFractionCSP3(mol),
        "NumAliphaticRings": rdMolDescriptors.CalcNumAliphaticRings(mol),
    }


# Calculate for all molecules
if 'output_sdf' in dir():
    mols = list(Chem.SDMolSupplier(output_sdf, removeHs=False))

    data_2d = []
    for mol in mols:
        if mol is None:
            continue
        name = mol.GetProp("_Name")
        desc = calculate_2d_descriptors(mol)
        desc["Name"] = name
        data_2d.append(desc)

    df_2d = pd.DataFrame(data_2d)
    df_2d = df_2d.set_index("Name")

    print("2D Descriptors:")
    print(df_2d.round(2).to_string())

3. 3D Descriptors (Geometrical)

These require 3D coordinates and describe molecular shape.

def calculate_3d_descriptors(mol):
    """
    Calculate 3D geometrical descriptors.

    These depend on the conformer and capture shape information.
    """
    descriptors = {}

    try:
        # Shape descriptors
        descriptors["RadiusOfGyration"] = rdMolDescriptors.CalcRadiusOfGyration(mol)
        descriptors["Asphericity"] = rdMolDescriptors.CalcAsphericity(mol)
        descriptors["Eccentricity"] = rdMolDescriptors.CalcEccentricity(mol)
        descriptors["InertialShapeFactor"] = rdMolDescriptors.CalcInertialShapeFactor(mol)
        descriptors["SpherocityIndex"] = rdMolDescriptors.CalcSpherocityIndex(mol)

        # Principal moments of inertia
        pmi = rdMolDescriptors.CalcPMI1(mol), rdMolDescriptors.CalcPMI2(mol), rdMolDescriptors.CalcPMI3(mol)
        descriptors["PMI1"] = pmi[0]
        descriptors["PMI2"] = pmi[1]
        descriptors["PMI3"] = pmi[2]

        # Normalized PMI ratios (useful for shape classification)
        if pmi[2] > 0:
            descriptors["NPR1"] = pmi[0] / pmi[2]  # Rod-likeness
            descriptors["NPR2"] = pmi[1] / pmi[2]  # Disk-likeness

        # Plane of best fit
        descriptors["PBF"] = rdMolDescriptors.CalcPBF(mol)

    except Exception as e:
        print(f"Error calculating 3D descriptors: {e}")

    return descriptors


if 'output_sdf' in dir():
    data_3d = []
    for mol in mols:
        if mol is None:
            continue
        name = mol.GetProp("_Name")
        desc = calculate_3d_descriptors(mol)
        desc["Name"] = name
        data_3d.append(desc)

    df_3d = pd.DataFrame(data_3d)
    df_3d = df_3d.set_index("Name")

    print("\n3D Shape Descriptors:")
    print(df_3d.round(3).to_string())

4. Surface Area Descriptors

Solvent-accessible surface area (SASA) is important for solubility and binding.

def calculate_surface_descriptors(mol):
    """
    Calculate surface area descriptors.
    """
    descriptors = {}

    try:
        # Need to compute radii first
        radii = rdFreeSASA.classifyAtoms(mol)

        # Total SASA
        sasa = rdFreeSASA.CalcSASA(mol, radii)
        descriptors["SASA_total"] = sasa

        # LabuteASA (approximate)
        descriptors["LabuteASA"] = rdMolDescriptors.CalcLabuteASA(mol)

        # PEOE VSA descriptors (partial charge weighted)
        peoe_vsa = rdMolDescriptors.PEOE_VSA_(mol)
        descriptors["PEOE_VSA_sum"] = sum(peoe_vsa)

        # SMR VSA descriptors (MR weighted)
        smr_vsa = rdMolDescriptors.SMR_VSA_(mol)
        descriptors["SMR_VSA_sum"] = sum(smr_vsa)

    except Exception as e:
        print(f"Error calculating surface descriptors: {e}")

    return descriptors


if 'output_sdf' in dir():
    data_surf = []
    for mol in mols:
        if mol is None:
            continue
        name = mol.GetProp("_Name")
        desc = calculate_surface_descriptors(mol)
        desc["Name"] = name
        data_surf.append(desc)

    df_surf = pd.DataFrame(data_surf)
    df_surf = df_surf.set_index("Name")

    print("\nSurface Area Descriptors:")
    print(df_surf.round(1).to_string())

5. Energy-Based Descriptors from Auto3D

The neural network potentials provide quantum-mechanical-like properties.

def calculate_energy_descriptors(mol):
    """
    Extract energy-based descriptors from Auto3D output.
    """
    descriptors = {}

    # Electronic energy from NNP
    if mol.HasProp("E_hartree"):
        e_hartree = float(mol.GetProp("E_hartree"))
        descriptors["E_hartree"] = e_hartree

        # Atomization energy (E / n_atoms) - rough measure of stability
        n_atoms = mol.GetNumAtoms()
        descriptors["E_per_atom"] = e_hartree * 627.509 / n_atoms  # kcal/mol per atom

    # Relative energy (if multiple conformers)
    if mol.HasProp("E_rel"):
        descriptors["E_rel"] = float(mol.GetProp("E_rel"))

    # Thermodynamic properties (if calculated)
    if mol.HasProp("G_hartree"):
        descriptors["G_hartree"] = float(mol.GetProp("G_hartree"))
    if mol.HasProp("H_hartree"):
        descriptors["H_hartree"] = float(mol.GetProp("H_hartree"))

    return descriptors


if 'output_sdf' in dir():
    data_energy = []
    for mol in mols:
        if mol is None:
            continue
        name = mol.GetProp("_Name")
        desc = calculate_energy_descriptors(mol)
        desc["Name"] = name
        data_energy.append(desc)

    df_energy = pd.DataFrame(data_energy)
    df_energy = df_energy.set_index("Name")

    print("\nEnergy-Based Descriptors:")
    print(df_energy.round(4).to_string())

6. Distance-Based Descriptors

Pharmacophore-relevant distances from 3D structures.

def calculate_distance_descriptors(mol):
    """
    Calculate distance-based descriptors.
    """
    conf = mol.GetConformer()
    n_atoms = mol.GetNumAtoms()

    # Get all pairwise distances
    distances = []
    for i in range(n_atoms):
        for j in range(i+1, n_atoms):
            pos_i = conf.GetAtomPosition(i)
            pos_j = conf.GetAtomPosition(j)
            dist = pos_i.Distance(pos_j)
            distances.append(dist)

    distances = np.array(distances)

    return {
        "MaxDist": distances.max(),           # Molecular diameter
        "MinDist": distances.min(),           # Closest atoms
        "MeanDist": distances.mean(),         # Average distance
        "StdDist": distances.std(),           # Distance variability
    }


if 'output_sdf' in dir():
    data_dist = []
    for mol in mols:
        if mol is None:
            continue
        name = mol.GetProp("_Name")
        desc = calculate_distance_descriptors(mol)
        desc["Name"] = name
        data_dist.append(desc)

    df_dist = pd.DataFrame(data_dist)
    df_dist = df_dist.set_index("Name")

    print("\nDistance-Based Descriptors (Å):")
    print(df_dist.round(2).to_string())

7. Combined Descriptor Table

Merge all descriptors for ML modeling.

if 'df_2d' in dir():
    # Combine all descriptors
    df_all = df_2d.copy()

    if 'df_3d' in dir():
        df_all = df_all.join(df_3d, how='left')
    if 'df_surf' in dir():
        df_all = df_all.join(df_surf, how='left')
    if 'df_dist' in dir():
        df_all = df_all.join(df_dist, how='left')
    if 'df_energy' in dir():
        df_all = df_all.join(df_energy, how='left')

    print(f"\nComplete Descriptor Matrix: {df_all.shape[0]} molecules × {df_all.shape[1]} descriptors")
    print("\nDescriptor columns:")
    for i, col in enumerate(df_all.columns):
        print(f"  {i+1:2d}. {col}")

# Export to CSV for ML
if 'df_all' in dir():
    output_csv = "/tmp/molecular_descriptors.csv"
    df_all.to_csv(output_csv)
    print(f"\nExported descriptors to: {output_csv}")

    # Show summary statistics
    print("\nDescriptor Statistics:")
    print(df_all.describe().round(2).T[["mean", "std", "min", "max"]].to_string())

8. Molecular Fingerprints

3D pharmacophore fingerprints capture spatial arrangement of features.

from rdkit.Chem import AllChem
from rdkit.Chem.Pharm2D import Gobbi_Pharm2D, Generate

def calculate_fingerprints(mol):
    """
    Calculate various molecular fingerprints.
    """
    fps = {}

    # Morgan (ECFP-like) - 2D
    fp_morgan = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=1024)
    fps["Morgan_OnBits"] = fp_morgan.GetNumOnBits()

    # RDKit fingerprint - 2D
    fp_rdkit = Chem.RDKFingerprint(mol)
    fps["RDKit_OnBits"] = fp_rdkit.GetNumOnBits()

    # 2D Pharmacophore
    try:
        factory = Gobbi_Pharm2D.factory
        fp_pharm = Generate.Gen2DFingerprint(mol, factory)
        fps["Pharm2D_OnBits"] = fp_pharm.GetNumOnBits()
    except:
        fps["Pharm2D_OnBits"] = np.nan

    return fps


if 'output_sdf' in dir():
    data_fp = []
    for mol in mols:
        if mol is None:
            continue
        name = mol.GetProp("_Name")
        fps = calculate_fingerprints(mol)
        fps["Name"] = name
        data_fp.append(fps)

    df_fp = pd.DataFrame(data_fp).set_index("Name")

    print("\nFingerprint Summary:")
    print(df_fp.to_string())

9. Descriptor Selection for QSAR

Best practices for selecting descriptors for ML models.

def select_descriptors(df, variance_threshold=0.01, correlation_threshold=0.95):
    """
    Select descriptors for ML modeling.

    Removes:
    - Constant or near-constant descriptors
    - Highly correlated descriptors
    """
    # Remove columns with NaN
    df_clean = df.dropna(axis=1)

    # Remove low variance
    variances = df_clean.var()
    high_var = variances[variances > variance_threshold].index
    df_clean = df_clean[high_var]

    # Remove highly correlated (keep first of each pair)
    corr_matrix = df_clean.corr().abs()
    upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(bool))
    to_drop = [col for col in upper.columns if any(upper[col] > correlation_threshold)]
    df_selected = df_clean.drop(columns=to_drop)

    return df_selected, {
        "original": len(df.columns),
        "after_nan": len(df_clean.columns) + len(to_drop),
        "after_variance": len(high_var),
        "after_correlation": len(df_selected.columns),
        "removed_correlated": to_drop
    }


if 'df_all' in dir():
    df_selected, info = select_descriptors(df_all)

    print("\nDescriptor Selection:")
    print(f"  Original: {info['original']} descriptors")
    print(f"  After removing NaN columns: {info['after_nan']}")
    print(f"  After variance filter: {info['after_variance']}")
    print(f"  After correlation filter: {info['after_correlation']}")
    print(f"\nRemoved (highly correlated): {info['removed_correlated']}")
    print(f"\nSelected descriptors: {list(df_selected.columns)}")

10. Best Practices

For QSAR/ML

1. Use 3D descriptors when conformation matters (binding, permeability)

2. Boltzmann-average for flexible molecules

3. Standardize descriptors before modeling (z-score)

4. Remove collinear descriptors to avoid overfitting

Descriptor Selection

• Start with comprehensive set (~200 RDKit descriptors)

• Remove constant/near-constant

• Remove highly correlated (keep more interpretable)

• Use feature importance from random forest for final selection

3D-Specific Considerations

• Shape descriptors (PMI, asphericity) are conformer-dependent

• Use lowest-energy conformer or Boltzmann average

• Energy-based descriptors from NNPs can capture electronic effects

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

Summary

This tutorial covered:

1. 2D descriptors - MW, LogP, TPSA, counts

2. 3D shape descriptors - RadiusOfGyration, Asphericity, PMI

3. Surface descriptors - SASA, VSA-based descriptors

4. Energy descriptors - from Auto3D NNP calculations

5. Distance descriptors - molecular dimensions

6. Fingerprints - Morgan, RDKit, pharmacophore

7. Descriptor selection - removing correlated features

Key workflow:

1. Generate 3D with Auto3D

2. Calculate descriptors with RDKit

3. Filter and select features

4. Use for QSAR/ML modeling