AutoDock Vina Docking Tutorial: From a PDB File to Publishable Results

Version changes

Current toolchain versus old tutorials

Most tutorials still rely on the AutoDockTools GUI (MGLTools 1.5.7) and Vina 1.1.2. The table lists what goes wrong when you copy them on Vina 1.2.x.

Installation

Install Vina 1.2.7, Meeko and supporting tools

The Vina pip package contains only the Python bindings, not the vina command. PyPI has prebuilt vina 1.2.7 wheels only for Linux x86_64 (Python 3.8–3.12); on macOS pip falls back to a source build and fails with Boost library location was not found! when Boost is missing. The conda-forge vina 1.2.7 package provides the vina and vina_split commands and the Python bindings on Linux and macOS, and the environment below uses it. Windows users run it in WSL or download the win.exe build from GitHub Releases. Tested on this page: part 1 of the commands below installs from scratch with micromamba on macOS arm64 (Apple M2).

Open Babel and RDKit play supporting roles in the current workflow. RDKit reads and writes SDF, assigns bond orders to PDB ligands and computes RMSD. Open Babel handles format conversion and is a PLIP dependency; it is no longer used to write PDBQT files for docking.

# 1) Environment: take vina from conda-forge; the package ships the vina and vina_split commands and the Python bindings
#    (linux-64, linux-aarch64, osx-64, osx-arm64; Python 3.10-3.14)
micromamba create -n dock -c conda-forge python=3.11 rdkit gemmi prody pdbfixer openmm vina=1.2.7 -y
micromamba activate dock
# Meeko 0.8.0 and molscrub 0.3.0 are only on PyPI; molscrub 0.3.0 imports joblib without declaring it, so install it by hand
pip install meeko==0.8.0 molscrub==0.3.0 propka joblib
vina --version        # AutoDock Vina v1.2.7

# 2) Without conda: download the binary for your platform from GitHub Releases
#    File suffixes: linux_x86_64, linux_aarch64, mac_aarch64, mac_x86_64, win.exe
curl -fL -o vina https://github.com/ccsb-scripps/AutoDock-Vina/releases/download/v1.2.7/vina_1.2.7_linux_x86_64
chmod +x vina && ./vina --version
#    pip install vina installs only the Python bindings; PyPI has prebuilt 1.2.7 wheels only for Linux x86_64,
#    so on macOS pip builds from source and fails with Boost library location was not found!

# 3) Optional: pocket prediction and interaction analysis
#    P2Rank 2.5.1: download from https://github.com/rdk/p2rank/releases and unpack (needs Java); the command is prank
#    fpocket: micromamba install -c conda-forge -c bioconda fpocket
#    PLIP: install openbabel from conda-forge, then pip install plip

Receptor preparation

Receptor preparation: waters, missing residues, protonation, metals and cofactors

Meeko matches each residue against templates. Heavy atoms must match the template exactly; hydrogens are optional. Receptor preparation therefore centers on completing heavy atoms and fixing the states of titratable residues.

Ligand preparation

Ligand preparation: 3D conformers, protonation and tautomers

mk_prepare_ligand.py requires a 3D structure with all hydrogens; SDF is recommended. PDB format has no bond-order information, and the official docs advise against using it for small molecules.

Docking straight from SMILES or a 2D structure fails for two reasons. Vina only rotates rotatable bonds during docking, so rings keep their input conformation and a wrong chair/boat cannot be corrected. A flat 2D input is also one of the failure causes listed in the official FAQ. molscrub embeds 3D conformers with ETKDGv3, minimizes them with a force field and converts six-membered boats to chairs.

Vina uses a united-atom scoring function, so output hydrogen positions are meaningless. Input hydrogens still decide which atoms are H-bond donors or acceptors, so the protonation state directly changes the score. For imatinib, molscrub gives two states over pH 5–9, with a neutral or a protonated pyridine, and the piperazine amines protonated in both.

• Structures downloaded from PubChem often include salts and counterions; scrub.py keeps only the largest fragment by default.
• When the protonation state is uncertain, enumerate states with --ph_low and --ph_high, dock each, and report the best one with its state.
• Resolve undefined stereocenters in the SMILES first. molscrub 0.3.0 has no stereoisomer enumeration option, and ETKDG assigns undefined stereocenters arbitrarily; when enumeration is needed, expand the SMILES with RDKit EnumerateStereoisomers first and pass the results to scrub.py.
• Macrocycles: Meeko makes macrocycles flexible by default, and Vina computes the conformational entropy penalty from the number of BRANCH records, which can shift scores by about 2 kcal/mol. Use the same setting for all molecules you compare, and add --rigid_macrocycles when needed.
• Ligands containing elements without Vina parameters or metals fail with Atom type ... is not a valid AutoDock type and need custom atom type parameters.

Grid box

How to set the docking box: center and size

A smaller box lets the search converge more easily; positions outside the box are never sampled. The Vina FAQ recommends a box "as small as possible, but not smaller" and suggests raising exhaustiveness when the box exceeds 30×30×30 Å.

Box too small: part of the ligand is pushed out, and Vina returns truncated poses or clearly worse scores. Box too large: at the same exhaustiveness the search is insufficient, and results vary strongly between seeds. Feinstein and Brylinski found on 3659 complexes that Vina poses were most accurate when the cube edge was 2.857 times the ligand radius of gyration; a generic box of "at least 22.5 Å per dimension" gave a mean RMSD of 4.9 Å versus 4.0 Å for the optimized box.

Blind docking (a box enclosing the whole protein) is common in network pharmacology papers. Che and Zhang reviewed 35 such papers from 2024 and found 17 using or likely using blind docking. In the CASF-2016 data they cite, Vina blind docking placed 34–47% of ligands within 2 Å RMSD, versus 90.2% when docking to the specified site.

Tested on this page (P2Rank 2.5.1, fpocket 4.2.3, input: chain A of 1IEP after PDBFixer as in the workflow below): the center of the top P2Rank pocket (probability 0.937) lies 3.0 Å from the crystal ligand center, the center of the top fpocket pocket lies 3.4 Å from it, and the pocket residues include Thr315 and Met318. --box_enveloping xtal_lig.sdf --padding 5 produced an 18.7×26.7×23.5 Å box; 2.857 × radius of gyration gives 19.1 Å for the crystal pose and 18.7 Å for the scrub.py conformer. P2Rank takes about 1 minute per structure on an 8-core machine; the column names in its CSV have leading spaces, so read it in pandas with skipinitialspace=True.

# Crystal structure
prank predict -f protA_H.pdb -o p2rank_out
# AlphaFold, NMR and cryo-EM models: the default model uses B-factors, which hold pLDDT in AlphaFold files
prank predict -c alphafold -f AF-model.pdb -o p2rank_out
# The first row is the top-ranked pocket; columns center_x/y/z give the box center
head -3 p2rank_out/protA_H.pdb_predictions.csv

# Cross-check: fpocket writes to protA_H_out/
fpocket -f protA_H.pdb

# Box from pocket center and ligand size (Feinstein & Brylinski 2015: edge = 2.857 x radius of gyration)
python - <<'EOF'
from rdkit import Chem
from rdkit.Chem import Descriptors3D
m = Chem.MolFromMolFile("lig_scrub.sdf")  # Hs removed by default: heavy atoms only
print("cube edge (A):", round(2.857 * Descriptors3D.RadiusOfGyration(m), 1))
EOF
mk_prepare_receptor.py --read_pdb protA_H.pdb -o rec -p -v \
    --box_center 12.3 45.6 7.8 --box_size 22 22 22

Parameters

exhaustiveness, num_modes, energy_range and random seeds

The same seed does not guarantee the same result as someone else: the random number generator comes from Boost, and builds against different Boost versions follow different trajectories. The maintainers suggest running about 10 times each and comparing the score distributions. In the data cited by Che and Zhang, raising exhaustiveness from 8 to 64 lowered the median RMSD from 3.37 Å to 2.21 Å at about 8 times the runtime.

Tested on this page (1IEP redocking, an 18.7×26.7×23.5 Å box from --box_enveloping with 5 Å padding; other conditions in the next section): exhaustiveness 8, 16, 32 and 64 with 3 seeds each (10 seeds for 8 and 32) all gave best-pose RMSDs of 0.77–0.85 Å and scores between −12.74 and −12.82 kcal/mol. The median CPU time (user) per run was 109, 281, 362 and 821 s, and the median wall time 70, 175, 94 and 280 s (heavily affected by machine load; use only as an order of magnitude). This box fits the crystal ligand tightly, and the default of 8 already reached 2 Å in 10 of 10 runs. The official docs note that the default 8 occasionally misses the correct pose, so large boxes or more flexible ligands should still start at 32. Repeated runs on the same machine with the same vina binary and the same seed give identical scores and poses.

Complete commands

Complete workflow: redocking imatinib into 1IEP

Using the official Vina example system 1IEP (c-Abl kinase + imatinib), these commands go from the PDB file to a table of scores and RMSDs. The ligand starts from SMILES, so the test checks whether the whole preparation workflow reproduces the crystal pose. In the official example the best pose scores about −13 kcal/mol with the Vina scoring function.

Tested on this page (AutoDock Vina 1.2.7 (conda-forge), Meeko 0.8.0, molscrub 0.3.0, PDBFixer 1.12, RDKit 2025.09 and 2026.03; macOS arm64, Apple M2 with 8 cores and 16 GB; other jobs were running on the machine during the tests, with a system load of about 25–70): the four code blocks above run end to end. All chain A residues from the PDBFixer output matched Meeko templates; scrub.py gives a single state at pH 7.4 (both piperazine N protonated, net charge +2). mk_export.py writes one SDF record per pose, not one molecule with several conformers.

Redocking results: starting from SMILES with exhaustiveness 32 and seeds 1–10, the best pose scored −12.74 to −12.82 kcal/mol with an in-place RMSD of 0.80–0.85 Å to the crystal pose, below 2 Å in all 10 runs. The second-ranked pose is the end-to-end flipped binding mode (RMSD about 13 Å, about 1.4 kcal/mol worse). PLIP shows that the top pose reproduces the hydrogen bonds to the Met318 backbone and the Thr315 side chain. Scoring the crystal pose directly (--score_only) gives −12.51 kcal/mol, and local optimization in place (--local_only) gives −13.24 kcal/mol at an RMSD of 0.21 Å.

The starting conformer changes the score: rerunning the same script in an RDKit 2026.03 environment, scrub.py produced a different starting conformer and the best score of all 3 seeds was −12.41 to −12.42 kcal/mol, with RMSD 0.84–0.85 Å. Cross-docking confirmed that the difference comes from the ligand starting conformer, not from receptor hydrogens. Vina keeps bond lengths, angles and ring conformations fixed during docking, so fix the ligand preparation environment or keep the starting conformer files when comparing scores.

The Vina terminal table lists up to num_modes poses, but the output PDBQT only contains poses within energy_range of the best score: here the terminal listed 9 poses and the file held only 2–7 in each run. rmsd_table.py reads scores from the REMARK VINA RESULT lines in the PDBQT, which correspond one to one with the SDF records.

set -euo pipefail
SMI='Cc1ccc(NC(=O)c2ccc(CN3CCN(C)CC3)cc2)cc1Nc1nccc(-c2cccnc2)n1'   # imatinib

# 1. Download 1IEP; keep chain A protein (no altloc or altloc A) and the chain A ligand STI
curl -fsSLO https://files.rcsb.org/download/1IEP.pdb
awk 'substr($0,1,4)=="ATOM" && substr($0,22,1)=="A" && (substr($0,17,1)==" " || substr($0,17,1)=="A")' 1IEP.pdb > protA.pdb
echo END >> protA.pdb
awk 'substr($0,1,6)=="HETATM" && substr($0,18,3)=="STI" && substr($0,22,1)=="A"' 1IEP.pdb > xtal_lig.pdb

# 2. Add missing heavy atoms, add hydrogens at pH 7.4, drop water and all heterogens
pdbfixer protA.pdb --output=protA_H.pdb --add-atoms=all --keep-heterogens=none --ph=7.4

# 3. Check pKa values of titratable residues in the pocket (written to protA_H.pka)
propka3 protA_H.pdb

# 4. Assign bond orders to the crystal ligand to get a reference SDF (used for the box and RMSD)
python make_ref.py "$SMI" xtal_lig.pdb xtal_lig.sdf

# 5. Receptor PDBQT + box: centered on the crystal ligand, padded by 5 A on each side
#    To change a protonation state add -n, e.g. -n A:<resnum>=HIP
mk_prepare_receptor.py --read_pdb protA_H.pdb -o rec -p -v \
    --box_enveloping xtal_lig.sdf --padding 5
#    Writes rec.pdbqt, rec.box.txt (Vina config) and rec.box.pdb (view the box in PyMOL)

# 6. Ligand: protonation state at pH 7.4 and 3D conformer from SMILES, then PDBQT
scrub.py "$SMI" -o lig_scrub.sdf --ph 7.4 --skip_tautomers
#    scrub.py writes unnamed molecules; without --multimol_prefix a single molecule becomes the hidden file lig_pdbqt/.pdbqt
mk_prepare_ligand.py -i lig_scrub.sdf --multimol_outdir lig_pdbqt --multimol_prefix lig

# 7. Redocking: 3 random seeds x exhaustiveness 32
mkdir -p out
for lig in lig_pdbqt/*.pdbqt; do
  name=$(basename "$lig" .pdbqt)
  for seed in 1 2 3; do
    vina --receptor rec.pdbqt --ligand "$lig" --config rec.box.txt \
      --exhaustiveness 32 --num_modes 9 --energy_range 3 --seed "$seed" \
      --out "out/${name}_s${seed}.pdbqt" | tee "out/${name}_s${seed}.log"
    mk_export.py "out/${name}_s${seed}.pdbqt" -s "out/${name}_s${seed}.sdf"
  done
done

# 8. Score and RMSD to the crystal pose for every pose
python rmsd_table.py xtal_lig.sdf out > redock_rmsd.tsv
sort -t$'\t' -k3,3g redock_rmsd.tsv | head

# make_ref.py: assign bond orders to the PDB ligand from a SMILES template
import sys
from rdkit import Chem
from rdkit.Chem import AllChem

smi, pdb_in, sdf_out = sys.argv[1:4]
template = Chem.MolFromSmiles(smi)
lig = Chem.MolFromPDBFile(pdb_in, removeHs=True)
ref = AllChem.AssignBondOrdersFromTemplate(template, lig)
Chem.MolToMolFile(ref, sdf_out)
print("heavy atoms:", ref.GetNumHeavyAtoms())

# rmsd_table.py: in-place (no superposition), symmetry-aware heavy-atom RMSD
import glob, os, sys
from rdkit import Chem
from rdkit.Chem import rdMolAlign

def heavy_neutral(m):
    # Remove H and formal charges so different protonation states still match
    m = Chem.RemoveHs(m)
    for a in m.GetAtoms():
        if a.GetFormalCharge() != 0:
            # A protonated N keeps one explicit H after RemoveHs; clearing only the charge raises a valence error
            a.SetFormalCharge(0)
            a.SetNumExplicitHs(0)
            a.SetNoImplicit(False)
    Chem.SanitizeMol(m)
    return m

ref = heavy_neutral(Chem.MolFromMolFile(sys.argv[1]))
print("file\tpose\tvina_score\trmsd_to_xtal")
for sdf in sorted(glob.glob(os.path.join(sys.argv[2], "*.sdf"))):
    pdbqt = sdf[:-4] + ".pdbqt"
    scores = [float(l.split()[3]) for l in open(pdbqt) if l.startswith("REMARK VINA RESULT")]
    poses = []
    for m in Chem.SDMolSupplier(sdf, removeHs=False):
        if m is None:
            continue
        for conf in m.GetConformers():
            p = Chem.Mol(m, confId=conf.GetId())
            poses.append(p)
    for i, p in enumerate(poses):
        rmsd = rdMolAlign.CalcRMS(heavy_neutral(p), ref)
        print(f"{os.path.basename(sdf)}\t{i+1}\t{scores[i]:.2f}\t{rmsd:.2f}")

# Batch docking: compute maps once and reuse them for every ligand
from vina import Vina
import glob, os

v = Vina(sf_name="vina", cpu=0, seed=42)
v.set_receptor("rec.pdbqt")
os.makedirs("out_batch", exist_ok=True)
# Without a ligand loaded, compute_vina_maps builds maps for every atom type in the force field, so they are reused across ligands
v.compute_vina_maps(center=[15.19, 53.90, 16.92], box_size=[22, 22, 22])  # replace with the values in rec.box.txt

for lig in sorted(glob.glob("lib_pdbqt/*.pdbqt")):
    v.set_ligand_from_file(lig)
    v.dock(exhaustiveness=16, n_poses=20)
    v.write_poses(os.path.join("out_batch", os.path.basename(lig)), n_poses=9, energy_range=3.0, overwrite=True)
    print(lig, v.energies(n_poses=1)[0][0])

Validation

Can you trust the result: redocking, positive controls and scoring limits

The −4.25, −5.0 and −7.0 kcal/mol thresholds common in network pharmacology papers have no single source, and scores from different docking programs cannot share a threshold. When comparing molecules of different size, also report ligand efficiency LE = −score / heavy-atom count; 0.3 kcal/mol per heavy atom is a common reference, and Vina scores become more negative as molecules grow.

1. 01

   Redocking RMSD < 2 Å

   Dock the co-crystal ligand back into its own structure. The symmetry-aware heavy-atom RMSD between the best pose and the crystal pose should be below 2 Å. This threshold appears in the original Vina paper (Trott & Olson 2010) and in the PoseBusters benchmark, where Vina met it for 58% of the 85 Astex Diverse systems. Compute RMSD in place in the receptor frame: RDKit CalcRMS does not superimpose, while GetBestRMS first aligns the ligand onto the reference and underestimates the value.

2. 02

   If redocking fails, check in this order

   Box units and position; ligand and receptor protonation states; higher exhaustiveness or other seeds; ligand ring conformations; receptor structure quality. If the crystal pose is not a minimum of the scoring function (compare vina --score_only and --local_only; both modes require the ligand coordinates to be inside the box, so prepare a PDBQT from the crystal ligand with hydrogens), a larger search will not help; change the scoring function or the method. Tested on 1IEP on this page: the crystal pose scores −12.51 with --score_only and −13.24 kcal/mol after --local_only, and the best docked pose scores −12.8 kcal/mol, with all three within 1 Å RMSD.

3. 03

   Positive control

   Dock a known active (ideally with an experimental Kd or IC50) with the same receptor, box and parameters, and compare new molecules against its score. In the network pharmacology papers reviewed by Che and Zhang, none used a positive control.

4. 04

   Agreement across seeds

   When the top poses from 3 seeds lie within 2 Å RMSD of each other, the search has converged. If changing only the seed changes the score a lot, shrink the box or raise exhaustiveness first.

Result analysis

Result analysis: choosing poses, PyMOL figures and 2D interaction diagrams

# Run in the PyMOL command line
load protA_H.pdb, rec
load out/lig-1_s1.sdf, poses
load xtal_lig.sdf, xtal
split_states poses, prefix=pose
delete poses
hide everything
show cartoon, rec
set cartoon_transparency, 0.5
select pocket, byres (rec within 4 of pose0001)
show sticks, pocket and not name N+C+O
show sticks, pose0001 or xtal
color grey70, xtal and elem C
color green, pose0001 and elem C
distance hb, pose0001, pocket, mode=2
label pocket and name CA, "%s%s" % (resn, resi)
orient pose0001
set ray_opaque_background, 0
png pose1.png, width=2400, height=1800, dpi=300, ray=1

# Export the complex for PLIP / LigPlot+: mark the ligand as HETATM with one residue name
alter pose0001, resn="LIG"
alter pose0001, chain="L"
alter pose0001, resi="900"
alter pose0001, type="HETATM"
save complex_pose1.pdb, rec or pose0001

# Interaction tables (txt + xml) and a PyMOL session; -p also renders images
plip -f complex_pose1.pdb -t -x -y -o plip_pose1

1. 01

   Read the output table

   affinity is the Vina score (kcal/mol). rmsd l.b. and rmsd u.b. are distances from the top pose of the same run and have nothing to do with the crystal structure; u.b. matches atoms one-to-one without symmetry, l.b. matches each atom to the nearest atom of the same element. The idea that "a smaller docking RMSD is better" confuses these two kinds of RMSD.

2. 02

   Choose a pose

   Do not just take rank 1. Check whether rank 1 recurs across seeds. When the top scores differ by much less than the Vina scoring error (2.85 kcal/mol), prefer the pose that reproduces known key interactions (in 1IEP, imatinib hydrogen-bonds to the hinge Met318 backbone and the gatekeeper Thr315). Discard poses that are mostly solvent-exposed or pressed against the box edge.

3. 03

   Export SDF

   mk_export.py restores bond orders, formal charges and all hydrogens from the SMILES stored in the PDBQT header. The SDF can go directly to PyMOL, PLIP, RDKit and downstream molecular dynamics.

4. 04

   3D figure

   Use the PyMOL commands below to draw the pocket, hydrogen bonds and an overlay with the crystal pose. If Vina 1.2.x output displays incompletely in PyMOL, the file usually contains NUL characters; view the SDF instead.

5. 05

   2D interaction diagrams

   PLIP reports residues and distances for hydrogen bonds, hydrophobic contacts, π stacking, salt bridges and halogen bonds. LigPlot+ and PoseView (proteins.plus) draw 2D diagrams, and Discovery Studio Visualizer can as well. All of them expect the receptor and ligand in one PDB file with the ligand as HETATM.

Writing the paper

Reporting docking results in a paper

Reviewers need the following to reproduce a docking study. A missing box center and size is the most common gap in network pharmacology papers.

• Software and versions: AutoDock Vina 1.2.7, Meeko 0.8.0, molscrub 0.3.0, PDBFixer, etc.
• Receptor: PDB ID, chain, resolution; which waters, ions and ligands were removed, which metals and cofactors were kept; which residues were rebuilt.
• Protonation: target pH, tools used (PDBFixer, PROPKA), residue states set manually.
• Ligands: source (e.g. PubChem CID), protonation and 3D generation method, number of states enumerated.
• Box: center coordinates and dimensions (Å), and the basis for the center (co-crystal ligand, literature residues, P2Rank).
• Parameters: scoring function, exhaustiveness, num_modes, energy_range, random seeds and number of repeats.
• Validation: redocking RMSD of the co-crystal ligand; positive control molecule and its score.
• Results: best score per molecule across seeds (mean ± SD), the basis for pose selection, interacting residue table and 3D figures.

Troubleshooting

Common errors and fixes

Tools common in China

Practical problems with tools common in China: CB-Dock2, metalloproteins, Discovery Studio and MD after docking

These lessons come from the Keinsci computational chemistry forum, Sobereva's blog and Chinese Q&A sites, and were checked against tool documentation or the original papers.

Hand off to an agent

Hand this workflow to an agent

Example instruction: "Use AutoDock Vina to dock these 5 compounds (SMILES attached) into the ATP pocket of PDB 1IEP. First validate by redocking the co-crystal imatinib, then use imatinib as the positive control, with 3 random seeds per molecule. Deliver a score table, an RMSD table, PyMOL figures and PLIP interaction tables."

The agent installs Vina 1.2.7, Meeko and molscrub in the workspace (RDKit, OpenMM and AmberTools are preinstalled), prepares the receptor and ligands as described here, redocks and checks the RMSD, runs the batch docking, exports SDF files, renders figures and compiles the interaction tables. Deliverables include scripts such as run_redock.sh, rec.pdbqt and rec.box.txt, docking outputs and logs for each seed, redock_rmsd.tsv, a score summary, PNG figures and PLIP reports. The agent reviews the results adversarially, for example checking that redocking really falls within 2 Å and that seeds agree.

You still need to check whether pocket residue protonation states match the experimental conditions, whether the box matches the binding site reported in the literature, whether the key interactions of the chosen pose have experimental support, and whether the positive control is appropriate. For molecular dynamics follow-up, the preinstalled GROMACS can run in the same workspace.