Bio Applied Docking

Molecular Docking with AutoDock Vina

When to Use

• Predicting the binding pose and affinity of a small molecule against a protein target
• Virtual screening a compound library (e.g. ZINC subset) against one receptor
• Redocking a co-crystallized ligand to validate a docking protocol before trusting new scores
• Preparing receptor/ligand PDBQT files from a PDB structure or SMILES strings
• Parsing and ranking vina log output across many docking runs

Version Compatibility

• AutoDock Vina ≥ 1.2 (output table format also matches 1.1.2)
• Open Babel ≥ 3.1 (obabel CLI)
• RDKit ≥ 2023.09
• Python ≥ 3.10, pandas ≥ 2.0

Prerequisites

• vina binary on PATH (conda install -c conda-forge vina or pip install vina), openbabel, rdkit, pandas
• A receptor structure (PDB, e.g. from RCSB) and ligand SMILES/SDF files
• Prior: bio-structural-biology-structure-io for fetching/cleaning the receptor PDB

Goal: Convert a receptor PDB and ligand SMILES into the PDBQT format Vina requires, filtering ligands for drug-likeness first.

Approach: Use RDKit to embed 3D conformers and compute Lipinski descriptors, then hand off charge assignment/PDBQT writing to Open Babel.

import subprocess
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors


def check_lipinski(mol):
    """Check Lipinski's Rule of Five. A molecule is 'drug-like' only if
    FEWER THAN 2 of the four rules are violated (a single violation is fine)."""
    mw = Descriptors.MolWt(mol)
    logp = Descriptors.MolLogP(mol)
    hbd = Descriptors.NumHDonors(mol)
    hba = Descriptors.NumHAcceptors(mol)
    violations = sum([mw > 500, logp > 5, hbd > 5, hba > 10])
    return {"MW": round(mw, 1), "LogP": round(logp, 2), "HBD": hbd, "HBA": hba,
            "violations": violations, "drug_like": violations < 2}


def smiles_to_ligand_pdbqt(smiles, out_path):
    """Embed a SMILES string in 3D and write a Vina-ready ligand PDBQT.
    Requires the `obabel` CLI on PATH for Gasteiger charges + PDBQT torsion tree.
    Returns the Lipinski check so callers can filter before docking.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        raise ValueError(f"Could not parse SMILES: {smiles}")
    lipinski = check_lipinski(mol)
    mol = Chem.AddHs(mol)
    AllChem.EmbedMolecule(mol, randomSeed=42)
    AllChem.MMFFOptimizeMolecule(mol)
    sdf_path = out_path.replace(".pdbqt", ".sdf")
    with Chem.SDWriter(sdf_path) as writer:
        writer.write(mol)
    subprocess.run(["obabel", sdf_path, "-O", out_path, "--partialcharge", "gasteiger"],
                    check=True)
    return lipinski


def prepare_receptor(pdb_path, out_pdbqt, strip_hetero=True):
    """Strip waters/heteroatoms and convert a receptor PDB to a rigid PDBQT.
    For production runs prefer AutoDockTools' prepare_receptor4.py, which assigns
    AD4 atom types more carefully than obabel's defaults; obabel is fine for quick work.
    """
    cmd = ["obabel", pdb_path, "-O", out_pdbqt, "-xr", "--partialcharge", "gasteiger"]
    if strip_hetero:
        cmd += ["--delete", "HOH"]
    subprocess.run(cmd, check=True)

Goal: Run AutoDock Vina docking with a grid box centered on the known (or predicted) binding site.

Approach: Call the vina CLI via subprocess, capturing stdout so it can be parsed downstream; keep the box just large enough to cover the pocket plus a few Å margin.

import subprocess


def run_vina(receptor_pdbqt, ligand_pdbqt, center, box_size=(20, 20, 20),
             out_path="docked.pdbqt", exhaustiveness=8, num_modes=9):
    """Dock one ligand into one receptor with AutoDock Vina.
    center/box_size are Angstrom coordinates/lengths of the search grid box.
    Default docking is rigid-receptor / flexible-ligand (Vina's default mode).
    Returns Vina's stdout table (mode | affinity | rmsd l.b. | rmsd u.b.).
    """
    cx, cy, cz = center
    sx, sy, sz = box_size
    cmd = ["vina", "--receptor", receptor_pdbqt, "--ligand", ligand_pdbqt,
           "--center_x", str(cx), "--center_y", str(cy), "--center_z", str(cz),
           "--size_x", str(sx), "--size_y", str(sy), "--size_z", str(sz),
           "--exhaustiveness", str(exhaustiveness), "--num_modes", str(num_modes),
           "--out", out_path]
    result = subprocess.run(cmd, check=True, capture_output=True, text=True)
    return result.stdout

Goal: Parse Vina's output table and rank many docked ligands (a virtual screen) by best-mode affinity.

Approach: Regex-free line parsing keyed on the leading mode number, then a pandas ranking with a configurable hit threshold.

import pandas as pd


def parse_vina_output(text):
    """Parse AutoDock Vina's stdout/log table into a tidy DataFrame.
    Columns: mode, affinity_kcal_mol, rmsd_lb, rmsd_ub (rmsd is vs. the top mode).
    """
    rows = []
    for line in text.strip().split("\n"):
        line = line.strip()
        if line and line[0].isdigit():
            parts = line.split()
            rows.append({
                "mode": int(parts[0]),
                "affinity_kcal_mol": float(parts[1]),
                "rmsd_lb": float(parts[2]),
                "rmsd_ub": float(parts[3]),
            })
    return pd.DataFrame(rows)


def rank_hits(ligand_best_scores, threshold_kcal_mol=-7.0):
    """Rank a batch of docked ligands by their best-mode affinity and flag hits.
    ligand_best_scores: dict[name] -> best-mode affinity (kcal/mol), e.g. one
    entry per ligand from parse_vina_output(text)['affinity_kcal_mol'].min().
    """
    hits = pd.DataFrame(
        [{"ligand": name, "affinity_kcal_mol": score}
         for name, score in ligand_best_scores.items()]
    )
    hits = hits.sort_values("affinity_kcal_mol").reset_index(drop=True)
    hits["hit"] = hits["affinity_kcal_mol"] <= threshold_kcal_mol
    return hits


def validate_by_redocking(rmsd_to_crystal_pose, cutoff_angstrom=2.0):
    """Sanity-check a docking protocol by redocking a co-crystallized ligand
    into its own site before trusting scores on new compounds.
    Convention: RMSD < 2 Å from the crystal pose = 'reproduced the binding mode'.
    """
    return rmsd_to_crystal_pose < cutoff_angstrom

Pitfalls

• Scoring-function reliability — Vina's empirical score ranks poses of the same ligand well (good for pose selection) but is unreliable for ranking different ligands against each other; large score gaps between compounds are often overestimation, not real potency.
• Always validate first — redock the native/co-crystal ligand before screening new compounds; if it doesn't reproduce the crystal pose (RMSD < 2 Å convention), fix receptor prep/grid box before trusting any hit list.
• Rigid receptor by default — Vina's default holds the receptor fixed and only samples ligand torsions, missing induced-fit effects; use flexible side chains or ensemble docking if the pocket is known to move.
• Covalent ligands — a covalently bound reference ligand (e.g. a Cys-Michael acceptor) is treated as non-covalent by standard Vina, so its reported affinity understates true potency; don't use it as-is for score calibration.
• Lipinski's Rule of Five — thresholds are MW > 500, LogP > 5, HBD > 5, HBA > 10; a compound only fails drug-likeness with 2 or more violations, not any single one.
• Grid box sizing — too small clips the true pocket and biases poses to whatever fits; too large wastes exhaustiveness on irrelevant space. Center on the known/predicted site with a few Å margin.
• Charge/protonation order — add hydrogens and assign protonation state before computing partial charges; charging an unprotonated structure gives wrong electrostatics in the PDBQT.

See Also

• bio-applied-molecular-modeling — force fields, energy minimization, homology modeling that precede docking
• bio-chemoinformatics-virtual-screening — scaling ligand libraries and filtering hit lists beyond Lipinski
• bio-chemoinformatics-molecular-descriptors — deeper descriptor/fingerprint computation for ligand triage
• bio-structural-biology-structure-io — fetching and cleaning receptor structures from the PDB