How to Set Protonation States and Add Hydrogens to a Protein Before Docking or MD (PDB2PQR, PROPKA and H++)

X-ray structures in the PDB almost never contain hydrogens, so you must add them yourself before docking or MD. Choose a pH, predict side-chain pKa values with PROPKA or H++, let PDB2PQR build the hydrogens at that pH, then check the histidine tautomers and any pKa near your working pH by hand before running anything.

A student downloads 1HPX, runs it straight through prepare_receptor, docks a ligand, and gets a number. The number is meaningless, because the aspartate dyad at the centre of the HIV protease active site was silently assigned the wrong charge. This is the single most common preparation error we see, and it costs nothing to avoid. This guide covers why the hydrogens are missing, which pH to use, the exact PDB2PQR and PROPKA commands, how to decide between HID, HIE and HIP, when H++ is the better tool, and how to verify the output before it reaches Vina or pdb2gmx. Written by the StemSkills Lab team, from 10+ years in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling.

Why does a PDB file have no hydrogens in it?

Because X-ray crystallography measures electron density, and a hydrogen atom has one electron. At the resolution of a typical protein structure, that density is not distinguishable from noise, so the depositors never model the hydrogens at all. Only ultra-high-resolution X-ray structures and neutron structures resolve them.

The scale of the gap is easy to underestimate. We downloaded 1HPX (HIV-1 protease with the inhibitor KNI-272) from the RCSB PDB and counted: 1,516 ATOM records, and exactly zero of them are hydrogens. After protonating it at pH 7.4 with PDB2PQR, the same structure has 3,129 atoms. More than half of the atoms in the model you eventually simulate were never measured. They were predicted, by whichever tool you happened to run.

That matters because charge drives almost everything downstream. Hydrogen bonds in a docking score, salt bridges in an MD trajectory, and the total system charge you have to neutralise with counter-ions all depend on which side chains carry a proton. If you have not yet dealt with gaps in the backbone, do that first: our guide on fixing missing residues and loops in a PDB structure comes before this step, not after.

What pH should I use for protonating a protein?

Use the pH your experiment or your reference paper used. If there is no experimental context, pH 7.4 is the standard choice for a physiological cytosolic or serum target, and pH 7.0 is the tool default in PDB2PQR. Lysosomal and endosomal proteins sit closer to pH 5.0, and gastric enzymes such as pepsin are active near pH 2.

Two points that beginners get wrong. First, pH is not a property of the structure file, so nothing in the PDB tells you what to pick. The crystallisation pH in the header is the pH of the drop, not of the biology, and it is often chosen for crystal quality rather than function. Second, the choice only changes residues whose predicted pKa is close to your pH. A lysine with a pKa of 10.5 is protonated at pH 5 and at pH 7.4 alike. The residues that flip are the ones within roughly one pH unit of your working value, and those are exactly the ones worth inspecting by hand.

How do I predict pKa values with PROPKA?

PROPKA is an empirical predictor: it starts from a model pKa for each residue type and applies structure-based perturbations for desolvation, hydrogen bonding and charge interactions. It is fast, running in seconds on a normal protein, which is why it is the default engine inside PDB2PQR.

Install and run it from the jensengroup/propka package:

pip install propka
propka3 1HPX.pdb

This writes 1HPX.pka next to the input. The section headed SUMMARY OF THIS PREDICTION is the part you read. Here is the real output for the catalytic dyad and the single histidine in that structure:

Read the two aspartates carefully. They are the same residue in the two halves of a homodimer, and PROPKA shifts one by 1.3 units and the other by 5.5 units. That asymmetry is the well-known monoprotonated state of the aspartyl protease dyad, and a tool that assigned both copies the same charge would have got the active site wrong. This is why you look at the numbers instead of accepting a default.

How accurate is the prediction? The PROPKA3 paper (Olsson, Søndergaard, Rostkowski and Jensen, Journal of Chemical Theory and Computation, 2011) reports agreement with experiment of “rmsd = 0.79 (0.91) for Asp and Glu, 0.75 (0.97) for Tyr, 0.65 (0.72) for Lys, and 1.00 (1.37) for His residues”, with the bracketed figures being the earlier version. Roughly one pKa unit of error on histidine is the number to keep in mind: it is the reason a histidine predicted at 7.0 is a decision for you to make, not for the software.

How do I add the hydrogens with PDB2PQR?

PROPKA predicts the pKa values. PDB2PQR uses them to build the actual hydrogens, optimise the hydrogen-bond network, and assign force-field charges and radii. Install it and run the two together:

pip install pdb2pqr

pdb2pqr30 --ff=AMBER --ffout=AMBER \
  --titration-state-method propka --with-ph 7.4 \
  --drop-water --keep-chain \
  --pdb-output 1HPX_H.pdb \
  1HPX.pdb 1HPX.pqr

What each flag does, quoting the tool’s own help text:

• --titration-state-method propka: “Method used to calculate titration states. If a titration state method is selected, titratable residue charge states will be set by the pH value supplied by –with_ph”. PROPKA is currently the only accepted value.
• --with-ph: “pH values to use when applying the results of the selected pH calculation method. (default: 7.0)”. Nothing happens with this flag unless you also set the titration state method.
• --ff: the force field whose charges and radii are written, one of AMBER, CHARMM, PARSE, TYL06, PEOEPB or SWANSON. The default is PARSE, which suits continuum electrostatics rather than MD.
• --ffout=AMBER: writes residue and atom names in that force field’s convention rather than the canonical PDB scheme. This is what turns a protonated aspartate into ASH and a delta-protonated histidine into HID.
• --pdb-output: “Create a PDB file based on input. This will be missing charges and radii”. Use this file for docking and for visual inspection; use the .pqr for electrostatics.
• --keep-chain and --drop-water: preserve chain IDs, and discard crystallographic waters before processing.

On our 1HPX run the log ends with Applying pKa values at a pH of 7.40, and the output contains exactly one ASH B 25 and two HID residues at position 69. The asymmetric dyad that PROPKA predicted was carried through into the structure automatically. There is also a browser version at the PDB2PQR web server if you cannot install anything locally, though the command line is easier to record in a methods section.

Which histidine tautomer should I choose: HID, HIE or HIP?

Histidine is the only common residue with a real tautomer decision, because its neutral form can carry the proton on either ring nitrogen. In AMBER naming the three states are:

• HID: neutral, proton on ND1 (the delta nitrogen). CHARMM calls this HSD.
• HIE: neutral, proton on NE2 (the epsilon nitrogen). CHARMM calls this HSE.
• HIP: both nitrogens protonated, net charge +1. CHARMM calls this HSP.

Pick HIP only when the predicted pKa is clearly above your working pH, which usually means the histidine sits in a negatively charged pocket or coordinates a phosphate. Between HID and HIE, the decision is made by the local hydrogen-bond network: whichever nitrogen has an acceptor pointing at it should keep the proton, and whichever has a donor nearby should stay bare. Metal coordination overrides everything, because a histidine ligating a zinc or a haem iron donates the lone pair on the coordinating nitrogen and must carry its proton on the other one.

GROMACS makes the same decision automatically inside pdb2gmx, and the documentation is explicit that “by default these selections are done automatically. For His, this is based on an optimal hydrogen bonding conformation.” You can take the decision back with -his, described as “Interactive histidine selection, instead of checking H-bonds”, and the sibling flags -asp, -glu, -lys, -arg and -gln do the same for the other titratable residues. Two warnings that catch people out. First, pdb2gmx checks geometry, not pKa, so it will never give you a HIP unless you ask. Second, if you have already protonated with PDB2PQR and you want those states preserved, do not pass -ignh, which is documented as “Ignore hydrogen atoms that are in the coordinate file” and will discard your work.

When is H++ a better choice than PDB2PQR?

Use H++ when you want a physics-based pKa rather than an empirical one, or when you want AMBER topology files produced in the same pass. H++ solves the Poisson-Boltzmann equation on the structure instead of applying fitted rules, which is slower but responds to the electrostatic environment directly.

The H++ 3.0 paper (Anandakrishnan, Aguilar and Onufriev, Nucleic Acids Research, 2012, W537 to W541) benchmarks the server on 23 high-quality structures containing 201 titratable groups. For groups with large pKa shifts it reports an RMS error of 1.44 pK units while still assigning the correct protonation state at pH 7 in 98% of cases. That distinction is the practical one: the state is usually right even when the number is not precise, and the state is what you feed to the simulation.

In practice these are complementary rather than competing. A workflow we use often is PROPKA for the numbers, PDB2PQR for the structure, then pdb2gmx with explicit residue choices for the topology. If you are still deciding which simulation you actually need, our comparison of molecular docking versus molecular dynamics is the place to start.

How do I check the protonated structure before docking or MD?

Four checks, in order, and none of them takes more than a minute.

1. Count the atoms. The protonated file should have roughly twice the atoms of the input. If the count barely moved, the hydrogens were not added.
2. Grep the residue names. With AMBER output naming, look for ASH, GLH, HID, HIE, HIP, LYN and CYM. Each one is a residue where the tool departed from the standard charge state, and each deserves a look: grep -E " (ASH|GLH|HID|HIE|HIP|LYN|CYM) " protein_H.pdb | cut -c18-27 | sort -u
3. Inspect the active site visually. Open the file in PyMOL or ChimeraX, show hydrogens on the binding-site residues, and confirm the donors and acceptors point where the chemistry says they should. Our computational biology skills roadmap covers the visualisation step in more depth.
4. Check the total charge. Sum the charge column of the .pqr file. It should be close to an integer, and it should match what the sequence implies. A non-integer total means a residue was parameterised badly.

What goes wrong, and how do I fix it?

Your ligand disappeared from the output. PDB2PQR does not know your ligand. On our 1HPX run the log printed Unable to find amino or nucleic acid definition for KNI. Parsing as new residue, and the inhibitor appears in the --pdb-output file but is absent from the .pqr entirely. Pass the ligand separately with --ligand ligand.mol2, or protonate the protein and the ligand independently and recombine them afterwards.

PROPKA reports a residue it could not identify. The same cause. Modified residues, cofactors and non-standard amino acids are outside the parameter set, so they are skipped rather than guessed at. Handle them explicitly instead of hoping.

Your histidine sits at exactly your working pH. HIS 69 in our run came out at 6.98 and 7.11 in the two chains, both within one unit of pH 7.4 and within the tool’s own stated error. When this happens, decide from the local hydrogen-bond geometry, state the choice in your methods, and if the residue is inside the binding site, run both tautomers and report the difference.

Docking fails or gives implausible scores. Check that the receptor really has hydrogens on it. The AutoDock Vina documentation is direct about this: the “receptor coordinate file must contain all hydrogen atoms”. The PDBQT conversion then keeps only the polar hydrogens and merges the rest, but it cannot merge hydrogens that were never there.

GROMACS rejects your protonated file. The residue names from --ffout=AMBER are AMBER names. If your MD force field is a CHARMM variant it expects HSD, HSE and HSP instead. Either write the output in the matching convention or let pdb2gmx build the hydrogens itself with explicit interactive choices. The full sequence from raw PDB to a running simulation is in our walkthrough on preparing a protein structure for molecular dynamics.

Frequently asked questions

Do I need to protonate the protein if I am only docking, not running MD?

Yes. Docking scoring functions count hydrogen bonds, and a hydrogen bond needs a hydrogen. AutoDock Vina requires a receptor file containing all hydrogen atoms before conversion to PDBQT, so skipping this step does not save you anything, it just moves the guess into a tool that guesses worse.

Can I just let AutoDock Tools or ChimeraX add hydrogens and skip PROPKA?

You can, and for a solvent-exposed surface it is usually fine. The risk is in the binding site, where buried or coupled residues have the largest pKa shifts, which are precisely the cases a geometry-only method cannot see. Running PROPKA first costs seconds and tells you whether the shortcut is safe here.

What is the difference between the .pqr and .pdb output of PDB2PQR?

The .pqr replaces the occupancy and B-factor columns with per-atom charge and radius, which is what Poisson-Boltzmann solvers such as APBS need. The .pdb written by --pdb-output keeps standard PDB formatting and is what you pass to docking or visualisation tools.

Should I keep crystallographic waters when I protonate?

Drop them by default with --drop-water, then add back only the waters you have a reason to keep, such as one bridging the ligand and the protein in the reference structure. Keeping all of them adds noise to docking without adding information.

How do I report this step in a thesis or paper?

State the structure and its PDB ID, the tool and version, the pKa method, and the pH. For example: hydrogens were added with PDB2PQR 3.7.1 using PROPKA 3.5.1 titration states at pH 7.4 with the AMBER force field, and the histidine tautomers were verified against the local hydrogen-bond network. A reviewer can reproduce that in one command.

The short version

The structure you download is incomplete in a way that is invisible until it changes your results. Pick a pH deliberately, predict the pKa values instead of assuming them, look at every residue the tool assigned a non-standard state, and pay particular attention to histidines and to any titratable residue inside the binding site. The whole procedure takes a few minutes and is the difference between a docking score you can defend and one you cannot.