Ligand Protonation States for Docking: Get pH 7.4 Right

To set a ligand’s protonation state for docking, start from SMILES and enumerate the states present at pH 7.4 with Dimorphite-DL (carboxylic acids become carboxylates, aliphatic amines become ammonium ions). Add RDKit tautomers if the scaffold has them, build 3D coordinates with explicit hydrogens, convert with Meeko, and check the PDBQT’s total charge before running Vina.

Most students fix the protein carefully, add hydrogens to histidines, and then dock the ligand exactly as PubChem drew it. PubChem’s records for aspirin (CID 2244) and imatinib (CID 5291) both show the neutral molecule with a formal charge of 0, yet at pH 7.4 aspirin’s carboxylic acid is almost entirely ionized and imatinib carries a positive charge on its piperazine ring. Dock the neutral forms and Vina sees a hydrogen bond donor that is not there, or misses one that is. This guide from the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modelling) walks through checking, setting and verifying ligand charge states with free tools, using the same imatinib and Abl kinase complex (PDB 1IEP) that the official AutoDock Vina tutorial uses. Every command and output below was run and checked with Open Babel 3.1.0, Dimorphite-DL 2.1.0, RDKit 2026.09.2 and Meeko 0.8.0.

Why does the ligand’s protonation state change a docking result?

Protonation decides which atoms can donate a hydrogen bond, which can accept one, and whether the ligand can form a salt bridge with an aspartate or a lysine. AutoDock Vina does not use the partial charges in your PDBQT file. Its official FAQ states that Vina ignores user-supplied partial charges and handles electrostatics through its hydrophobic and hydrogen-bonding terms (the scoring function is described in Trott and Olson, J. Comput. Chem. 31:455-461, 2010). That makes the hydrogens themselves the thing that matters: an OH written as a donor atom (type HD on the hydrogen) scores differently from a bare carboxylate oxygen.

You can see this directly. We built aspirin twice, neutral and as its anion, and converted both with Meeko. The neutral PDBQT has 14 atoms, including one polar hydrogen typed HD. The anion has 13 atoms and no HD at all, with four oxygen acceptors (type OA). Same molecule, different set of interactions available to the scoring function. The Vina documentation puts it plainly in its ligand preparation section: “Your success can sometimes hang by just an hydrogen atom.”

Which functional groups change charge at pH 7.4?

A group is mostly ionized when the pH is more than about one unit past its pKa. You do not need to memorise pKa tables, but you should recognise the usual suspects before trusting any tool’s output.

The groups with a pKa near 7.4 (some piperazines, imidazoles, sulfonamides with electron-withdrawing rings) are the ones where you genuinely need more than one state.

How do you set the ligand’s charge at pH 7.4 with Open Babel?

Open Babel’s -p option adds hydrogens appropriate for a given pH, using the transforms in its phmodel.txt file, as described in the obabel documentation. Write the pH explicitly:

obabel ligands.smi -osmi -p 7.4
obabel -:"CC(=O)Oc1ccccc1C(=O)O" -osdf --gen3d -p 7.4 -O aspirin_ph74.sdf

With -p 7.4, aspirin becomes CC(=O)Oc1ccccc1C(=O)[O-] and dopamine becomes [NH3+]CCc1ccc(O)c(O)c1, which is what you want. Two cautions from our own runs:

• Always type the number. In Open Babel 3.1.0, running -p with no value left aspirin’s carboxylic acid neutral, while -p 7.4 ionized it. Do not rely on the default.
• Open Babel works atom by atom. Its documentation warns that molecules with several ionizable centres “will have all centers ionized”. For imatinib, -p 7.4 puts a positive charge on both piperazine nitrogens, a total charge of +2. That is unlikely: once one nitrogen is protonated, the second becomes much harder to protonate. A single fixed answer is the weakness of this method.

For simple ligands with one ionizable group, Open Babel is fine. For anything with two or more nearby centres, enumerate states instead.

How do you enumerate protonation states with Dimorphite-DL?

Dimorphite-DL (Ropp, Kaminsky, Yablonski and Durrant, J. Cheminform. 11:14, 2019) uses substructure rules with pKa ranges and returns every state that is plausible inside a pH window. Install and run it as shown in the Dimorphite-DL README:

pip install dimorphite_dl
dimorphite_dl --ph_min 7.4 --ph_max 7.4 "CC(=O)Oc1ccccc1C(=O)O"

Aspirin returns one state, the carboxylate. Imatinib is more instructive:

dimorphite_dl --ph_min 7.4 --ph_max 7.4 "CN1CCN(CC1)Cc1ccc(cc1)C(=O)Nc1ccc(C)c(Nc2nccc(n2)-c2cccnc2)c1"

It returns four states: neutral, protonated on the N-methyl nitrogen, protonated on the benzylic nitrogen, and doubly protonated. Dimorphite-DL is telling you that both piperazine nitrogens fall inside its uncertainty window, so it will not choose for you.

The pH window matters a lot. The defaults are --ph_min 6.4 and --ph_max 8.4, with --precision 1.0 (standard deviations of the pKa estimate). Run imatinib with those defaults and you get 64 states, because the aromatic nitrogens of the pyridine and pyrimidine rings and the anilino NH also land inside the wider window. The README lists weakly basic aromatic nitrogens such as pyridine as a known case that comes back in both states at the default range, even though the cation is a minor species at physiological pH. Other options you will use:

• --max_variants caps the output per molecule (default 128).
• --label_states tags each output as PROTONATED, DEPROTONATED or BOTH.
• --output_file states.smi writes the results to a file.

For docking, set both --ph_min and --ph_max to your assay or physiological pH and treat the result as a shortlist to reason about, not a list to dock blindly.

How do you enumerate tautomers with RDKit?

Tautomers move a hydrogen and a double bond together, so they change donors and acceptors without changing the total charge. 2-hydroxypyridine and 2-pyridone are the textbook pair. RDKit’s TautomerEnumerator in rdkit.Chem.MolStandardize.rdMolStandardize lists them:

from rdkit import Chem
from rdkit.Chem.MolStandardize import rdMolStandardize

mol = Chem.MolFromSmiles("Oc1ccccn1")   # 2-hydroxypyridine as drawn
te = rdMolStandardize.TautomerEnumerator()
for t in te.Enumerate(mol):
    print(Chem.MolToSmiles(t))
print("canonical:", Chem.MolToSmiles(te.Canonicalize(mol)))

This prints three tautomers (Oc1ccccn1, O=c1cccc[nH]1 and the non-aromatic O=C1CC=CC=N1) and picks the pyridone as canonical. On other molecules the counts climb quickly: warfarin gave 9 tautomers and guanine 15 in our run. The enumerator’s limit is reported by te.GetMaxTautomers() (1000 by default) and changed with te.SetMaxTautomers().

Two things to know before you dock any of them:

• The canonical tautomer is a label, not a prediction. The RDKit documentation notes that the canonical form is very likely not the most stable one under any given conditions; its purpose is to give the same answer whatever the input. For guanine, RDKit’s canonical form is not the 1H,9H tautomer found in DNA.
• Many enumerated forms are high in energy. Imatinib returns 6 tautomers, and four of them turn the amide into an imidic acid (N=C(O)), which you would not dock. Use chemistry to prune: keep amides as amides, keep aromatic rings aromatic unless the literature says otherwise, and keep the forms reported for your compound class.

How do you choose which states to dock?

Use evidence in this order:

1. Measured pKa values for your compound, from the published literature.
2. A co-crystal structure. If your ligand or a close analogue is in the PDB, look at what the ionizable atoms touch. In PDB 1IEP (mouse Abl kinase with imatinib, ligand code STI, 2.10 Å), the N-methyl nitrogen of the piperazine sits 2.67 Å from the backbone carbonyl oxygen of Ile360 and 3.19 Å from that of His361 in chain A. A nitrogen that close to two carbonyl oxygens is acting as a hydrogen bond donor, which only works if it is protonated. The benzylic nitrogen has no polar partner within 3.6 Å.
3. Rule-based enumeration (Dimorphite-DL, Molscrub) when you have neither.

The official AutoDock Vina basic docking tutorial agrees with the crystal: its prepared 1iep_ligand.sdf has a total charge of +1, with the hydrogen on the N-methyl nitrogen. That is one of the four states Dimorphite-DL returned, and the only one the structure supports.

When the evidence is weak, dock the two or three most plausible states as separate ligands and report each. Do not average their scores. Pick the state whose pose and interactions make chemical sense, and say in your methods section which states you tested and why.

Tap to explore in 3D

Abl kinase with imatinibThe binding pocketPiperazine and its carbonyl partners

How do you get from SMILES states to Vina-ready PDBQT files?

This script takes the Dimorphite-DL states, builds 3D coordinates with explicit hydrogens in RDKit, and writes one multi-molecule SDF:

from rdkit import Chem
from rdkit.Chem import AllChem
from dimorphite_dl import protonate_smiles

smiles = "CN1CCN(CC1)Cc1ccc(cc1)C(=O)Nc1ccc(C)c(Nc2nccc(n2)-c2cccnc2)c1"  # imatinib
states = protonate_smiles(smiles, ph_min=7.4, ph_max=7.4)

writer = Chem.SDWriter("imatinib_states.sdf")
for i, smi in enumerate(states):
    mol = Chem.AddHs(Chem.MolFromSmiles(smi))
    AllChem.EmbedMolecule(mol, randomSeed=42)
    AllChem.MMFFOptimizeMolecule(mol)
    mol.SetProp("_Name", f"imatinib_state{i}")
    print(mol.GetProp("_Name"), Chem.GetFormalCharge(mol), smi)
    writer.write(mol)
writer.close()

Then convert every molecule with Meeko (our Meeko preparation guide covers the options in detail):

mk_prepare_ligand.py -i imatinib_states.sdf --multimol_outdir pdbqt_out

This writes imatinib_state0.pdbqt to imatinib_state3.pdbqt. Now verify the charge. The partial charge sits in columns 71 to 76 of each atom line, so this one-liner sums it:

for f in pdbqt_out/*.pdbqt; do
  awk '/^(ATOM|HETATM)/{s+=substr($0,71,6)} END{printf "%s %.3f\n", FILENAME, s}' "$f"
done

Our run gave sums of -0.001, 0.999, 0.999 and 1.999, matching formal charges of 0, +1, +1 and +2. A sum that does not round to the charge you intended means a hydrogen was lost or added on the way. If you want one tool to do enumeration and 3D generation together, the Forli lab’s Molscrub (scrub.py, with a --ph option) is the route the Vina documentation suggests for ligands without hydrogens.

What goes wrong most often, and how do you fix it?

• Meeko prints “RDKit molecule has implicit Hs. Need explicit Hs.” and writes no PDBQT. Your SDF has no hydrogens. Add them with Chem.AddHs() before embedding, or with obabel ... -p 7.4, and regenerate 3D coordinates.
• The carboxylic acid is still neutral after Open Babel. You ran -p without a number. Use -p 7.4 and check the SMILES output before building 3D.
• Every amine is charged and the ligand is +2 or +3. That is Open Babel’s atom-by-atom model. Enumerate with Dimorphite-DL and choose using pKa data or a crystal structure.
• Dimorphite-DL returns dozens of states. You left the default 6.4 to 8.4 window. Set --ph_min and --ph_max to the same value, or lower --precision.
• Bond orders and charges disappear after conversion. You prepared the ligand from a PDB file. The Vina documentation advises against PDB format for small molecules because it does not store bond information. Start from SMILES or SDF. Our guide to extracting a ligand from a PDB file shows how to recover bond orders.
• The PDBQT charge sum is off by 1. A hydrogen was added or stripped by a format conversion. Rebuild from SMILES and run the awk check again.

For problems that only show up once Vina runs, see common AutoDock Vina errors and their fixes.

Where does this fit in a docking workflow?

Ligand protonation is one half of preparation. The other half is the receptor, covered in our guide to protein protonation states, and the two meet in preparing a protein and ligand for docking. The full sequence, from choosing a target to validating poses, is on our molecular docking learning path, and the wider set of skills is mapped in the computational biology skills roadmap. If you plan to take the docked complex into MD, the ligand’s charge state also sets its total charge in the force field, so the choice you make here carries through to ligand parameterization.

Frequently asked questions

Does AutoDock Vina use the partial charges in my ligand PDBQT?

No. The Vina FAQ states that it ignores user-supplied partial charges. Protonation still matters, because hydrogens decide which atoms are typed as donors (HD) and acceptors (OA, NA), and those types drive the hydrogen-bond term.

Should I always dock the ligand at pH 7.4?

Use the pH of the biological compartment or the assay you are comparing against. Most cytosolic and plasma targets use 7.4, but lysosomal or stomach targets sit lower. Set Dimorphite-DL’s --ph_min and --ph_max to that value.

Is Open Babel’s -p option good enough?

For a ligand with one ionizable group, yes, as long as you write -p 7.4 explicitly. For two or more nearby groups it charges all of them, so enumerate with Dimorphite-DL and choose.

How many protonation states or tautomers should I dock?

Usually one to three per ligand. Prune with measured pKa values, co-crystal contacts and basic chemistry, dock the survivors as separate ligands, and report which ones you tested.

Can I use the RDKit canonical tautomer for docking?

Only after checking it. RDKit’s documentation says the canonical tautomer is chosen for consistency, not stability, so confirm it matches the form reported for your compound class.