# Bio Applied Docking

> Dock ligands into a receptor with AutoDock Vina: build PDBQT files (Open Babel/RDKit), set the grid box, run vina, parse/rank poses by affinity and RMSD. Use for molecular docking, virtual screening, redocking a co-crystal ligand, or PDB/SMILES to PDBQT conversion.

- Skill: `pavel-kravchenko/bio-applied-docking` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-docking`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-docking/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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.

```python
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.

```python
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.

```python
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