Redocking with AutoDock Vina 1.2 and Meeko: a reproducible 1HSG walkthrough

Most docking tutorials end the moment Vina prints a table of scores. That is exactly where the real work begins. A docking score means nothing until you have shown that your setup can reproduce a known answer. This walkthrough runs the full modern AutoDock Vina pipeline (Meeko for preparation, Vina 1.2 for docking, RDKit for analysis) on one of the most-used validation systems there is: HIV-1 protease bound to indinavir (PDB 1HSG, ligand code MK1).

The result is more instructive than a clean success. With a reasonable box and default-ish settings, Vina’s top-ranked pose was 10.5 Å away from the crystal pose. A pose at 0.87 Å was in the output, but ranked fourth. Every number below comes from an actual run, and the commands are given in full so you can reproduce them.

What you need

Create a clean environment. Mixing old and new versions of Meeko is a common source of confusing import errors. For example, Meeko 0.5 fails with No module named 'rdkit.six' on current RDKit, because rdkit.six was removed years ago.

micromamba create -n dock -c conda-forge \
    python=3.12 "meeko=0.8" "vina=1.2.7" \
    rdkit gemmi openbabel scipy
micromamba activate dock
python -c "import meeko, vina, rdkit; \
print(meeko.__version__, vina.__version__, rdkit.__version__)"
# 0.8.0 1.2.7 2026.03.1

Download the structure:

wget https://files.rcsb.org/download/1HSG.pdb -O 1hsg.pdb

Step 1: split the complex

1HSG contains two protease chains (A and B, 99 residues each), the inhibitor MK1, and 127 crystallographic waters. For a first redocking we keep only the protein:

grep -E '^(ATOM|TER)' 1hsg.pdb > protein.pdb && echo END >> protein.pdb
grep '^HETATM.\{11\}MK1' 1hsg.pdb > mk1_xtal.pdb      # 45 heavy atoms

Removing every water is a choice, not a law. HIV protease has a famous structural “flap water” (HOH 308 in this entry, about 3 Å from the backbone N of both Ile50s) that bridges the inhibitor carbonyls to the flap residues Ile50/Ile50′. Many protocols keep it. We will come back to this when we interpret the results.

Step 2: build a correct ligand

PDB files store no bond orders and, in most entries, no hydrogens. If you feed mk1_xtal.pdb straight into a preparation tool, it has to guess the aromaticity, the amide bonds and the protonation from the geometry. Sometimes it guesses wrong without telling you. The robust approach is to assign bond orders from a trusted SMILES and then add hydrogens with 3D coordinates:

# make_ligand.py
from rdkit import Chem
from rdkit.Chem import AllChem

# Connectivity of indinavir (MK1). Stereo is then taken from the crystal coordinates.
template = Chem.MolFromSmiles(
    "CC(C)(C)NC(=O)C1CN(Cc2cccnc2)CCN1CC(O)CC(Cc1ccccc1)C(=O)NC1c2ccccc2CC1O")
xtal = Chem.MolFromPDBFile("mk1_xtal.pdb", removeHs=True)
xtal = AllChem.AssignBondOrdersFromTemplate(template, xtal)
Chem.AssignStereochemistryFrom3D(xtal)
xtal_h = Chem.AddHs(xtal, addCoords=True)
Chem.MolToMolFile(xtal_h, "mk1_xtal_h.sdf")          # reference pose, 92 atoms

# Redocking must NOT start from the answer: build a fresh 3D conformer.
smi = Chem.MolToSmiles(xtal)
mol = Chem.AddHs(Chem.MolFromSmiles(smi))
AllChem.EmbedMolecule(mol, randomSeed=42)
AllChem.MMFFOptimizeMolecule(mol)
Chem.MolToMolFile(mol, "mk1_start.sdf")
print(smi)

Two details matter here:

1. Stereochemistry comes from the crystal. The template SMILES has no stereo flags. AssignStereochemistryFrom3D reads the five stereocentres from the coordinates, so you cannot accidentally dock the wrong diastereomer. A SMILES typed from memory is the classic way to get this wrong. Check smi against the PDB Chemical Component Dictionary entry for MK1.
2. Protonation. Indinavir’s most basic nitrogen (the piperazine, pKa around 6) is mostly neutral at pH 7.4, so we dock the neutral form (formal charge 0). For your own ligands, decide the protonation state deliberately, for example with a pKa predictor or obabel -p 7.4, rather than accepting whatever the input file implies.

The starting conformer is 2.65 Å from the crystal geometry after optimal superposition, so the docking has to do real conformational work.

Step 3: prepare the ligand with Meeko

mk_prepare_ligand.py -i mk1_start.sdf -o mk1.pdbqt

Meeko assigns AutoDock atom types and Gasteiger charges, merges non-polar hydrogens, and builds the torsion tree. Indinavir has 14 active torsions. That is high: the search space grows quickly with every rotatable bond, and anything above about 10 deserves extra exhaustiveness. Meeko also writes the SMILES and atom mapping into the PDBQT REMARK lines, which is what later lets mk_export.py turn docked poses back into proper SDF molecules with bond orders.

Why not Open Babel (obabel -xr) or the old MGLTools prepare_ligand4.py? Both still work, but they cannot round-trip bond orders. MGLTools is also Python 2 software that no longer installs cleanly. Use Open Babel for format conversion (obabel in.mol2 -O out.sdf), not for the final PDBQT.

Step 4: prepare the receptor and define the box

The box is where most docking runs go wrong. Here we used the Vina box generator on the crystal ligand with 5 Å padding on each side, which gives:

That is a volume of about 9,360 ų: big enough for the ligand to rotate freely, small enough to keep the search focused. Meeko can write the receptor PDBQT and the box file in one step:

mk_prepare_receptor.py --read_pdb protein.pdb -o 1hsg_rec -p -v \
  --box_center 11.634 22.470 5.859 --box_size 24.854 20.214 18.635
# -> 1hsg_rec.pdbqt, 1hsg_rec.box.txt (Vina config), 1hsg_rec.box.pdb (for viewing)

Meeko matches each residue against chemical templates and adds the polar hydrogens Vina needs (328 HD atoms here). If it stops with a template error, the cause is almost always a structural problem in the PDB: a missing side chain, a non-standard residue, or an alternate location. Fix the structure (PDBFixer, Modeller), or use -d to delete the residue on purpose. Do not work around it by stripping atoms.

Always open 1hsg_rec.box.pdb together with the protein in PyMOL or ChimeraX before you dock. Ten seconds of looking will save you days.

Step 5: dock

vina --receptor 1hsg_rec.pdbqt --ligand mk1.pdbqt --config 1hsg_rec.box.txt \
     --exhaustiveness 16 --num_modes 9 --seed 42 --cpu 8 --out mk1_out.pdbqt

This took 91 seconds on 8 cores. Always set --seed so that a run can be repeated exactly. The output:

mode |   affinity | dist from best mode
     | (kcal/mol) | rmsd l.b.| rmsd u.b.
-----+------------+----------+----------
   1       -10.42          0          0
   2       -10.23      2.131      10.62
   3       -10.12      2.203      10.65
   4       -10.06      2.102      10.39
   5       -9.882       1.41      3.869
   ...

Note that the “rmsd” columns compare each pose with the best mode, not with the crystal structure. They tell you nothing about accuracy.

Step 6: measure the RMSD to the crystal pose

Convert the poses back to SDF and compute a symmetry-aware RMSD without superposition. The protein frame is fixed, so the poses must not be re-aligned. Re-aligning hides placement errors.

mk_export.py mk1_out.pdbqt -s mk1_out.sdf

# rmsd.py
from rdkit import Chem
from rdkit.Chem import rdMolAlign
ref = Chem.RemoveHs(Chem.MolFromMolFile("mk1_xtal_h.sdf"))
for i, pose in enumerate(Chem.SDMolSupplier("mk1_out.sdf"), 1):
    rms = rdMolAlign.CalcRMS(Chem.RemoveHs(pose), ref)   # symmetry-aware, no fitting
    print(f"pose {i}: {rms:.2f} Å")

The usual success criterion is a top-ranked pose within 2.0 Å. By that criterion this run failed, even though a near-perfect pose (0.87 Å) was found.

Step 7: diagnose the result (scoring or sampling?)

When redocking fails, there are two possible culprits. Either the scoring function prefers a wrong pose (a scoring failure), or the search never found the right one well enough (a sampling failure). Vina can tell you which by scoring the crystal pose directly:

mk_prepare_ligand.py -i mk1_xtal_h.sdf -o mk1_xtal.pdbqt
vina --receptor 1hsg_rec.pdbqt --ligand mk1_xtal.pdbqt \
     --config 1hsg_rec.box.txt --score_only
#   Estimated Free Energy of Binding : -10.950 kcal/mol
vina --receptor 1hsg_rec.pdbqt --ligand mk1_xtal.pdbqt --config 1hsg_rec.box.txt \
     --local_only --out mk1_xtal_min.pdbqt
#   Estimated Free Energy of Binding : -11.391 kcal/mol   (0.29 Å from the crystal)

The crystal pose, after a local optimisation that moves it only 0.29 Å, scores −11.39 kcal/mol. That is a full kcal/mol better than anything the global search reported. The scoring function is fine. The search did not converge. The 10.5 Å top pose is the inhibitor bound backwards in this C2-symmetric active site. Its centroid is only 1.3 Å from the crystal centroid, but the tert-butyl group sits where the aminoindanol should be, and the aminoindanol sits in the pyridyl pocket. Pseudo-symmetric sites and pseudo-symmetric peptidomimetics make this a well-known trap.

More exhaustiveness alone did not fix it. With --exhaustiveness 32 and three different seeds, the top pose was the backwards pose (−10.43, 10.55 Å) every time. The 0.86 Å pose showed up within the top 4 in only one of the three runs. Practical options, roughly in order of effort:

• Raise the exhaustiveness much further, or run several seeds and pool the poses. With 14 torsions this is expected, not a sign of a broken setup.
• Tighten the box. Use 3–4 Å of padding instead of 5, but only after you have confirmed that the native pose fits.
• Keep structural waters such as the flap water, either as part of the receptor or with Vina’s hydrated docking protocol (Forli et al., 2016 / Eberhardt et al., 2021).
• Rescore the top N poses with a second method, or check the binding modes visually. For example, see whether the central hydroxyl sits between the catalytic Asp25/Asp25′, as it does in every hydroxyethylamine inhibitor structure.

The general lesson: always redock a known ligand first, and always compare the crystal pose score with the docked scores. If the crystal pose scores better than your top pose, increase sampling. If it scores worse, no amount of sampling will help, and you need a different scoring approach or receptor setup.

Interpreting Vina scores

Vina reports kcal/mol, but the score is an empirical function fitted to reproduce affinities on average. Its typical error is around 2 kcal/mol, which is more than a factor of 20 in K_d. The score difference between ranks 1 and 4 here, 0.36 kcal/mol, is far inside that noise. Use scores to rank poses of one ligand and, with care, to triage many ligands. Do not use them to claim that compound A “binds 3× better” than compound B. For a feel for what kcal/mol differences mean in affinity terms, see the K_d ↔ ΔG calculator. For a more physics-based (but still approximate) estimate, run MD on the complex and then MM-PBSA with gmx_MMPBSA.

Reproducibility checklist

• Record the exact versions (vina --version, meeko.__version__), the box, the exhaustiveness and the seed.
• Keep the prepared PDBQT files. Small preparation differences change the results.
• Report the redocking RMSD of the top pose and of the best pose, together with the crystal-pose score.
• Do not report more decimal places than the method supports. Two is plenty.

References

1. Eberhardt J, Santos-Martins D, Tillack AF, Forli S. AutoDock Vina 1.2.0: New docking methods, expanded force field, and Python bindings. J Chem Inf Model 2021;61:3891–3898. doi:10.1021/acs.jcim.1c00203
2. Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem 2010;31:455–461. doi:10.1002/jcc.21334
3. Forli S, Huey R, Pique ME, Sanner MF, Goodsell DS, Olson AJ. Computational protein–ligand docking and virtual drug screening with the AutoDock suite. Nat Protoc 2016;11:905–919. doi:10.1038/nprot.2016.051
4. Chen Z, Li Y, Chen E, et al. Crystal structure at 1.9-Å resolution of human immunodeficiency virus (HIV) II protease complexed with L-735,524, an orally bioavailable inhibitor of the HIV proteases. J Biol Chem 1994;269:26344–26348. (PDB 1HSG)
5. Meeko documentation: meeko.readthedocs.io. AutoDock Vina documentation: autodock-vina.readthedocs.io.
6. O’Boyle NM, Banck M, James CA, et al. Open Babel: An open chemical toolbox. J Cheminform 2011;3:33. doi:10.1186/1758-2946-3-33