A small molecule nestled into a glowing protein binding pocket

Docking asks a narrow, useful question: given a protein structure and a small molecule, where does the molecule sit, and how favourable does that placement look? AutoDock Vina answers it with a search plus a scoring function, in seconds per compound. It is genuinely useful for triage — and routinely over-interpreted as a binding affinity prediction, which it is not.

This tutorial takes one real example end to end: imatinib into the ABL kinase domain (PDB 1IEP), including the preparation steps where almost every beginner’s result goes wrong, and the redocking validation that tells you whether to trust the run at all.

What Vina does and does not tell you

Vina gives youVina does not give you
A ranked set of poses in a defined boxA measured binding free energy
A score in kcal/mol, comparable within one runA score comparable across different proteins
Implicit flexibility of the ligandReceptor side-chain flexibility
Fast triage of hundreds of compoundsSolvent, entropy or membrane effects

The score is a heuristic. Two compounds differing by 0.3 kcal/mol are indistinguishable by this method, and a compound that scores −10 in one protein cannot be compared to one that scores −8 in another. Used as a ranking filter inside a consistent protocol, though, it is a reasonable tool.

Install the tools

conda create -n dock -c conda-forge -c bioconda \
  autodock-vina openbabel meeko pymol-open-source
conda activate dock

vina --version
obabel -V

Note the toolchain shift: Vina 1.2 replaced the old AutoDockTools ligand scripts with meeko, while receptor preparation still most commonly uses AutoDockTools or ADFR Suite. Understanding which script belongs to which version prevents a lot of confusion when following older tutorials.

Prepare the receptor

Preparation is where docking succeeds or fails, and it is boring work. A raw PDB file contains water, sometimes a co-crystallised ligand, alternate conformations for some residues, and no hydrogens. Docking into that structure gives plausible-looking nonsense.

# 1. Download the structure
wget https://files.rcsb.org/download/1IEP.pdb

# 2. Inspect: what is in the file besides protein?
grep '^HETATM' 1IEP.pdb | awk '{print $4}' | sort | uniq -c
#   1 STI   <- imatinib, the ligand we want to remove and later redock
#  many HOH  <- crystallographic waters

Remove waters and the co-crystallised ligand, keep any structurally important metal ions, resolve alternate locations, and add hydrogens at pH 7.4. PyMOL can do the removal cleanly:

pymol -cq <<'EOF'
load 1IEP.pdb, cplx
remove solvent
remove resn STI
remove alt B
h_add
save receptor_clean.pdb
EOF

Then convert to PDBQT, which is PDB plus partial charges and AutoDock atom types. Vina cannot read a PDB directly.

# AutoDockTools route (script ships with MGLTools / ADFR Suite)
prepare_receptor4.py -r receptor_clean.pdb -o receptor.pdbqt -A hydrogens

# If you only have Vina 1.2+, meeko can prepare receptors from a PDB too
mk_prepare_receptor.py --read_pdb receptor_clean.pdb -o receptor_prepared -p

Two checks before moving on. Confirm the file has HD or hydrogen atoms and AD4-style types, and confirm there are no remaining HETATM records you did not intend. A single stray water in the pocket changes everything.

Prepare the ligand

Starting from a SMILES string is the most reproducible route, because it documents exactly which molecule you docked. Open Babel generates a 3D conformer; meeko then produces the PDBQT with rotatable bonds identified.

# Imatinib, from its canonical SMILES
SMILES='CC1=C(C=C(C=C1)NC(=O)C2=CC=C(C=C2)CN3CCN(CC3)C)NC4=NC=CC(=N4)C5=CC=CN=C5'

obabel -:"$SMILES" -O ligand.sdf --gen3d --best

# Vina 1.2 toolchain
mk_prepare_ligand.py -i ligand.sdf -o ligand.pdbqt

# Legacy toolchain (Vina 1.1.x), for comparison
# prepare_ligand4.py -l ligand.sdf -o ligand.pdbqt

Check the protonation state. Imatinib’s piperazine nitrogen is basic and will be protonated at physiological pH; docking the neutral form can change both the pose and the score. If the protonation matters for your question, generate it explicitly rather than accepting the default.

How many rotatable bonds?

grep -c '^BRANCH' ligand.pdbqt
# 7 for imatinib

Every rotatable bond multiplies the search space. Seven is manageable; twenty is a compound that needs a much larger --exhaustiveness and should be treated with more scepticism when it reports a confident pose.

The box is not a technical detail, it is your hypothesis. Too small and the molecule cannot fit; too large and you waste search effort on solvent while the reported score describes a pose outside the site you care about.

# Where is the co-crystallised ligand's centre? Use it as the box centre.
python - <<'EOF'
from Bio.PDB import PDBParser
import numpy as np
s = PDBParser(QUIET=True).get_structure("x", "1IEP.pdb")
atoms = [a.coord for a in s[0]["A"] if a.get_parent().id == "STI"]
print(np.round(np.mean(atoms, axis=0), 2))
EOF
# [ 15.19  15.31  27.36]

Centring on the known ligand is the honest approach for redocking validation, because it defines the site by experimental evidence rather than by guesswork. A box of 20 Å per side comfortably accommodates a drug-sized molecule at this site.

cat > box.txt <<'EOF'
center_x = 15.19
center_y = 15.31
center_z = 27.36
size_x = 20
size_y = 20
size_z = 20
EOF

Run the docking

vina \
  --receptor receptor.pdbqt \
  --ligand ligand.pdbqt \
  --config box.txt \
  --out docked.pdbqt \
  --log vina.log \
  --exhaustiveness 32 \
  --num_modes 9 \
  --seed 42 \
  --cpu 8
OptionMeaningPractical advice
--exhaustiveness 32Search effort per ligand (default 8)32 for careful work, 64+ for flexible ligands
--num_modes 9How many poses to keepLook at all of them, not just mode 1
--seed 42Random seed for reproducibilityAlways set it, always record it
--cpu 8Parallel threadsMatch to available cores, not to cores you wish you had

The log reports each mode’s score and, crucially, the RMSD lower and upper bounds relative to the best mode. A first mode at −11.2 with nine near-identical poses is a converged result. Nine poses scattered by 4 Å mean the search did not converge and the top score is luck.

cat vina.log | sed -n '/mode/,$p' | head -12

Validate before you trust anything

Redocking validation is the difference between docking and pretending. Take the co-crystallised ligand, dock it back into its own structure, and compare the predicted pose with the experimental one. If your protocol cannot reproduce a pose you already know, its predictions for unknown compounds carry no weight.

# Compute RMSD between the top docked pose and the crystal pose
python - <<'EOF'
from pymol import cmd
cmd.load("docked.pdbqt", "docked")
cmd.load("1IEP.pdb", "crystal")
cmd.remove("crystal and not resn STI")
cmd.extract("crystal_pose", "crystal and resn STI")
# align on the crystal ligand, then measure
r = cmd.align("docked", "crystal_pose")
print("RMSD:", round(r[0], 2), "Angstrom")
EOF

Below about 2 Å is a pass. Between 2 and 3 Å is ambiguous and needs inspection. Above 3 Å, fix the preparation before running a single compound — the usual causes are missing hydrogens, a wrong protonation state, an oversized box, or a receptor still containing the original ligand.

The failures you will actually hit

SymptomLikely cause
“Unknown atom type” on receptorPDB converted without charges, not prepared
Every pose in the same trivial spotBox off-centre or far too large
Scores identical to two decimals for unrelated compoundsSearch not converging; raise exhaustiveness
Top pose outside the intended pocketBox includes solvent or a second cavity
Redocking RMSD above 3 ÅPreparation error, not a scoring failure

Docking many compounds without fooling yourself

Once the protocol validates, batch docking becomes a triage exercise. Keep the receptor, box, exhaustiveness and seed fixed across all compounds, because any change makes the scores incomparable. Rank, take the top few percent, and then interrogate the poses visually: contact residues, hydrogen bonds, and whether the molecule is even in the site you meant.

A docking score is a hypothesis generator, not a result. Turning it into a claim means experimental testing or at minimum a second, independent method such as molecular dynamics or free-energy calculation. If you want to see how docking fits into a wider structure-based analysis — from sequence to model to validated interaction — Ampersand Academy covers that progression.

Practice

  1. Reproduce the redocking validation for 1IEP and report the RMSD to two decimals.
  2. Re-run with the box enlarged to 40 Å and describe what changes in the pose distribution.
  3. Dock the neutral and protonated forms of imatinib and compare scores — then decide which is defensible.
  4. Run the same ligand twice with different seeds and check whether the top score is stable.

Is the Vina score a binding affinity?

No. It is a heuristic scoring function value in kcal/mol, useful for ranking compounds within one protocol. It omits entropy, explicit solvent and receptor flexibility, and it is not comparable between different proteins or even reliably between different box definitions.

What RMSD counts as a successful redocking?

Below 2 Å between the predicted and experimental pose is the usual threshold for success. Between 2 and 3 Å is borderline and should be inspected manually; above 3 Å indicates a preparation or search problem rather than a biological finding.

Why must I remove the co-crystallised ligand before docking?

Leaving it in the structure blocks the binding site and forces Vina to search elsewhere, producing confident poses in the wrong place. Remove it, then use its coordinates to define the box centre and, separately, to validate the protocol by redocking.

Should I use a rigid or flexible receptor?

Standard Vina treats the receptor as rigid, which is a real limitation for pockets that close on binding or for resistance mutations. If side-chain motion matters, consider flexible side chains with supported Vina variants or follow up with molecular dynamics rather than trusting a single rigid run.

How do I choose exhaustiveness?

Start at 32 for routine work and 64 or higher for flexible ligands with many rotatable bonds. Judge the setting by convergence: if the top nine poses cluster tightly and the score is stable across seeds, the search found the minimum. If poses scatter, increase effort before interpreting anything.

Can I dock into a homology model?

Yes, but the uncertainty compounds. A model built at 30 percent identity can have mispositioned side chains in the pocket, and docking will happily produce precise-looking poses in an imprecise site. Validate the model first and treat docking results as much weaker evidence.

Does Vina consider water molecules in the pocket?

Not by default — water is removed during preparation and treated as bulk solvent implicitly at best. Bridging waters can be essential for binding, so for well-characterised sites consider retaining specific structural waters or using a method that models them explicitly.

Last updated on · Written by