Molecular Docking Basics

Lesson 5 of 12 · 10 min

Preparing the Ligand

Docking scores a ligand against a receptor, but the structure you start with is never ready as-is: it may be a flat 2D graph, it lacks hydrogens, it has no partial charges, and the docking engine has no idea which bonds it is allowed to rotate. This lesson turns any starting ligand into a docking-ready PDBQT with correct 3D geometry, explicit hydrogens, Gasteiger charges, and a torsion tree. We do it with Open Babel, so the whole job is one command line.

Get the ligand

You obtain a ligand one of two ways: extract it from a co-crystal structure, where it sits as HETATM records under a three-letter residue code, or build it from a SMILES string or an SDF you fetched from PubChem or ChEMBL. A crystal ligand arrives with a real bound 3D pose but no hydrogens, while a SMILES string is only a flat connectivity graph with neither coordinates nor hydrogens. A quick grep on the residue name pulls the ligand out, though for messy files a structure-aware tool like PyMOL avoids grabbing the wrong atoms.

bash
# Route A: carve the ligand out of a co-crystal PDB (STI = imatinib in PDB 1IEP)
grep '^HETATM' complex.pdb | grep ' STI ' > ligand.pdb

# Route B: start from a SMILES string (acetylsalicylic acid, i.e. aspirin)
echo "CC(=O)Oc1ccccc1C(=O)O aspirin" > aspirin.smi

Add 3D, hydrogens, charges

In one pass, --gen3d builds a single low-energy 3D conformer from the flat SMILES graph, -h adds explicit hydrogens, --partialcharge gasteiger computes the Gasteiger partial charges, and the PDBQT writer assigns AutoDock atom types plus the torsion tree. That charge flag is not optional: leave it off and Open Babel still writes a valid-looking file, but every atom lands with a meaningless +0.000 charge. When the ligand came straight out of a crystal complex it already carries its bound pose, so drop --gen3d and run obabel ligand.pdb -O ligand.pdbqt -h --partialcharge gasteiger to keep that geometry. Confirm the protonation state suits your target: -h protonates for a neutral molecule, so use -p 7.4 instead when you need a physiological pH.

bash
obabel aspirin.smi -O aspirin.pdbqt --gen3d -h --partialcharge gasteiger

Rotatable bonds and torsions

A rotatable bond, or torsion, is a single acyclic bond that lets the fragment on one side swing relative to the other. Open Babel flags every such bond between two non-terminal heavy atoms, and it also flags the polar hydrogens on terminal -OH and -NH groups, since swinging a hydrogen-bond donor changes how the ligand binds; it skips only symmetric spins like a terminal methyl, which change nothing. It then roots the molecule on one rigid ROOT fragment, hangs the flexible pieces off it as BRANCH to ENDBRANCH blocks, and records the total as TORSDOF on the last line. Keeping these torsions active lets the docking engine sample the ligand's real conformations; freezing too many hides the bound pose, while inventing extra ones inflates the search space and wastes sampling.

bash
# Count the rotatable bonds Open Babel wrote into the torsion tree
grep -c '^BRANCH ' aspirin.pdbqt

# ...and read the torsional degrees of freedom on the last line
grep '^TORSDOF' aspirin.pdbqt
4
TORSDOF 4

Try it yourself: build dopamine with echo "NCCc1ccc(O)c(O)c1 dopamine" > dopamine.smi, convert it with obabel dopamine.smi -O dopamine.pdbqt --gen3d -h --partialcharge gasteiger, then run grep '^TORSDOF' dopamine.pdbqt. You get TORSDOF 5: three torsions down the aminoethyl chain (the ring-CH2, CH2-CH2, and CH2-NH2 bonds, the last one swinging the amine hydrogens) plus one for each catechol -OH, while only the aromatic ring stays rigid. Compare that to aspirin's 4 and predict which ligand hands the search engine more conformations to explore before you check.