Quantum Chemistry Features for ML in Drug Discovery

Classical ML drug discovery models use ECFP4 fingerprints or 2D graph representations. These encode topology but lack electronic information. Quantum chemistry descriptors capture reactivity, metabolic vulnerability, membrane permeability drivers, and binding electrostatic complementarity with higher physical grounding.

The practical challenge: reference-level DFT (B3LYP/6-31G*) takes minutes per molecule. GFN2-xTB (Grimme group, ORCA/xtb) is a semi-empirical tight-binding DFT method accurate to ~1 kcal/mol for drug-sized molecules and runs in ~1–3 seconds — practical for library-scale feature extraction.

QM Descriptors from xTB (GFN2-xTB)

Extracting xTB Descriptors

import subprocess
import json
import tempfile
import os
from pathlib import Path
from typing import Optional
from rdkit import Chem
from rdkit.Chem import AllChem


def smiles_to_xyz(smiles: str, charge: int = 0) -> Optional[str]:
    """
    Generate a 3D XYZ block from SMILES using RDKit ETKDG embedding.

    Args:
        smiles: Input SMILES string
        charge: Molecular charge for embedding

    Returns:
        XYZ format string, or None if embedding fails.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    mol = Chem.AddHs(mol)
    result = AllChem.EmbedMolecule(mol, AllChem.ETKDGv3())
    if result != 0:
        return None

    AllChem.MMFFOptimizeMolecule(mol, maxIters=2000)

    conf = mol.GetConformer()
    atoms = [mol.GetAtomWithIdx(i).GetSymbol() for i in range(mol.GetNumAtoms())]
    n_atoms = len(atoms)

    xyz_lines = [str(n_atoms), f"Generated from SMILES: {smiles}"]
    for atom, (x, y, z) in zip(atoms, [conf.GetAtomPosition(i) for i in range(n_atoms)]):
        xyz_lines.append(f"{atom:2s}  {x:12.6f}  {y:12.6f}  {z:12.6f}")

    return "\n".join(xyz_lines)


def run_xtb(xyz_content: str, charge: int = 0, multiplicity: int = 1) -> Optional[dict]:
    """
    Run xTB GFN2 semi-empirical calculation and extract electronic properties.

    Requires xtb binary in PATH (https://github.com/grimme-lab/xtb).

    Args:
        xyz_content: XYZ format molecular geometry
        charge: Molecular charge
        multiplicity: Spin multiplicity (1 = singlet)

    Returns:
        Dict of QM descriptors, or None if xTB fails.
    """
    with tempfile.TemporaryDirectory() as tmpdir:
        xyz_path = os.path.join(tmpdir, "molecule.xyz")
        with open(xyz_path, "w") as f:
            f.write(xyz_content)

        cmd = [
            "xtb", xyz_path,
            "--gfn", "2",
            "--chrg", str(charge),
            "--uhf", str(multiplicity - 1),
            "--json",
            "--alpb", "water",  # ALPB implicit solvation
            "--opt",            # Geometry optimisation
        ]

        result = subprocess.run(
            cmd, capture_output=True, text=True, cwd=tmpdir, timeout=120
        )

        if result.returncode != 0:
            return None

        # xtb writes xtbout.json
        json_path = os.path.join(tmpdir, "xtbout.json")
        if not os.path.exists(json_path):
            return None

        with open(json_path) as f:
            xtb_data = json.load(f)

    return {
        "total_energy_hartree": xtb_data.get("total energy"),
        "homo_energy_ev": xtb_data.get("HOMO") * 27.2114 if xtb_data.get("HOMO") else None,  # Hartree → eV
        "lumo_energy_ev": xtb_data.get("LUMO") * 27.2114 if xtb_data.get("LUMO") else None,
        "homo_lumo_gap_ev": xtb_data.get("HOMO-LUMO gap") * 27.2114 if xtb_data.get("HOMO-LUMO gap") else None,
        "dipole_debye": xtb_data.get("dipole"),
        "polarisability_au": xtb_data.get("polarizability"),
        "partial_charges": xtb_data.get("partial charges"),
        "n_atoms": xtb_data.get("number of atoms"),
    }

Delta-ML: Correcting Fast Methods with ML

import numpy as np
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import cross_val_score

def build_delta_ml_model(
    smiles_list: list[str],
    y_dft: np.ndarray,
    y_xtb: np.ndarray,
    fingerprint_matrix: np.ndarray,
    xtb_features: np.ndarray,
) -> dict:
    """
    Train a delta-ML model to correct xTB predictions towards DFT quality.

    Delta-ML paradigm:
        y_pred = y_xtb + delta_model(features)
        where delta = y_dft - y_xtb

    This is significantly more data-efficient than direct DFT prediction:
    the model only needs to learn the residual, which is smoother than
    the full property surface.

    Args:
        smiles_list: List of compound SMILES
        y_dft: High-level DFT target values (reference, expensive)
        y_xtb: xTB baseline predictions for same set
        fingerprint_matrix: ECFP4 fingerprints [n_compounds, 2048]
        xtb_features: xTB electronic descriptors [n_compounds, n_xtb_feats]

    Returns:
        Dict with trained model, CV scores, and feature importance.
    """
    delta = y_dft - y_xtb  # What the ML model must learn

    # Combine fingerprints and xTB features
    X = np.hstack([fingerprint_matrix, xtb_features])

    model = GradientBoostingRegressor(
        n_estimators=500,
        max_depth=4,
        learning_rate=0.05,
        subsample=0.8,
        min_samples_leaf=5,
        random_state=42,
    )

    cv_scores = cross_val_score(
        model, X, delta,
        cv=5,
        scoring="neg_mean_absolute_error",
        n_jobs=-1,
    )

    model.fit(X, delta)

    return {
        "delta_model": model,
        "cv_mae_mean": float(-cv_scores.mean()),
        "cv_mae_std": float(cv_scores.std()),
        "description": "Predicts Δ = DFT - xTB; add to xTB prediction for final value",
    }


def predict_delta_ml(
    xtb_value: float,
    features: np.ndarray,
    delta_model,
) -> float:
    """Apply delta-ML correction: y_pred = y_xtb + delta_model(features)."""
    delta = delta_model.predict(features.reshape(1, -1))[0]
    return float(xtb_value + delta)

QM9 Benchmark and Model Comparison

MAE values are approximate; benchmark results vary by train/test split. QM9 dataset: 134k drug-like molecules with DFT-computed properties.

Key Takeaways

• GFN2-xTB provides physically meaningful electronic descriptors (HOMO/LUMO, charges, dipole) in ~1–3 seconds per molecule — the practical sweet spot for lead optimisation campaigns of 10k–100k compounds.
• Delta-ML (ML correcting xTB towards DFT) typically requires only 500–2000 DFT reference calculations to reach sub-0.05 eV MAE on electronic properties — 50–100× fewer than direct ML-from-DFT.
• ANI neural network potentials (TorchANI) are best for conformer energetics and geometry optimisation at near-DFT accuracy. They are not suitable for properties requiring orbital information (HOMO/LUMO).
• SchNet/DimeNet++ require explicit 3D coordinates. Generating conformers adds computational overhead but enables distance and angle features unavailable to 2D graph models. Use ETKDG + MMFF for conformer generation.
• QM features add the most value for: metabolic site prediction (CYP450 SOM), covalent fragment screening (electrophilicity), CNS drug design (dipole-driven permeability), and reactive metabolite flags.

Try it live

Explore a related model hands-on: open Quantum Chemical Reaction Visualization and change its parameters while it runs. Everything stays in your browser.