protonate-utils 0.1.3

Add hydrogens to ligands and proteins at a target pH.

protonate_utils

A single utility for adding hydrogens to ligands and proteins at a target pH, for use in molecular modeling and structure-based drug design.

Why this exists

Most structures you download (a ligand from a database, a protein from the PDB) are missing hydrogens, or carry hydrogens that don't reflect the protonation state at physiological pH. Getting these right matters: a carboxylic acid is deprotonated (-COO⁻) at pH 7.4, a basic amine is protonated (-NH₃⁺), and a histidine side chain can go either way. Downstream tasks (docking, free-energy calculations, MD simulations, electrostatics) all depend on the correct charge and hydrogen placement.

Ligands and proteins need different tools for this. Small molecules are best handled with cheminformatics pKa models; proteins need residue-aware logic and geometry-based hydrogen placement. protonate_utils.py wraps the appropriate specialist tool for each case behind one consistent interface, so you don't have to remember two separate workflows:

• Ligands use Dimorphite-DL for pH-aware protonation states and the RDKit for structure handling. When the input has 3D coordinates, the heavy-atom geometry is preserved exactly; only the newly added hydrogens are given computed positions.
• Proteins use Hydride for geometry-based hydrogen addition and Biotite for PDB handling, with formal charges estimated per amino acid at the requested pH.

Everything is exposed both as a command-line tool and as an importable Python API.

Installation

Install the latest release from PyPI with pip:

pip install protonate-utils

Or install from a checkout for development:

git clone https://github.com/PatWalters/protonate_utils
cd protonate_utils
pip install -e .

Either way installs the dependencies for both modes (RDKit + Dimorphite-DL for ligands, Biotite + Hydride + NumPy for proteins), puts a protonate-utils command on your PATH, and makes import protonate_utils available.

Command-line usage

Once installed, use the protonate-utils command. The first argument selects the mode: ligand or protein. (You can also run it without installing via python protonate_utils.py … from a checkout.)

Ligands

# SDF in, SDF out (3D coordinates preserved, hydrogens placed from geometry)
protonate-utils ligand input.sdf output.sdf

# SMILES in, SMILES out, at a custom pH
protonate-utils ligand input.smi output.smi --ph 7.4

# Mixed: read SDF, write SMILES
protonate-utils ligand input.sdf output.smi

Input and output formats are inferred from the file extension: .smi/.smiles is treated as SMILES, anything else as SDF. SMILES files are read one molecule per line as SMILES [optional name]; an optional header line (e.g. SMILES Name) is detected and skipped automatically. When the output is SDF, explicit hydrogens are written even for SMILES input (positioned from the geometry when the input was 3D, otherwise added without coordinates).

Option	Default	Description
--ph	7.4	Target pH for protonation.

Molecules that fail to parse or protonate are skipped with a warning on stderr; the run reports how many were read, written, and skipped.

Proteins

# Remove a bound ligand by residue name, then add hydrogens
protonate-utils protein input.pdb AP5 output.pdb

# Remove several ligands at once with a comma-delimited list. The bundled
# 7axj_protein.pdb has two pocket ligands, EST and CL6:
protonate-utils protein 7axj_protein.pdb "EST,CL6" 7axj_out.pdb

# Keep everything (no ligand removal)
protonate-utils protein input.pdb none output.pdb --ph 7.0

The second positional argument is the residue name (3-letter CCD code) of a ligand to remove before protonation. Pass a comma-delimited list (e.g. "EST,CL6") to remove several residues at once — handy for clearing both a cofactor and a buffer/ion from a pocket — or none to keep all atoms. An error is raised naming any residue that isn't present. Output hydrogens are reordered so each one immediately follows the heavy atom it is bonded to.

Option	Default	Description
--ph	7.0	pH used to estimate amino-acid formal charges.
--no-relax	off	Skip dihedral relaxation of the added hydrogens.
--no-honor-protonation	off	Ignore force-field protonation names and (re)protonate every residue from --ph.

Force-field protonation/tautomer residue names

Inputs from force fields often name residues by their protonation or tautomer state — HID/HIE/HIP (and CHARMM HSD/HSE/HSP), ASH, GLH, LYN, ARN, CYM, TYM, CYX. These are not CCD codes, so Biotite cannot assign their bonds (HID/HIE even collide with unrelated CCD entries), which previously left those residues unbonded and unprotonated. They are now normalized to their canonical CCD residue, and by default the exact state the name encodes is preserved, overriding the pH estimate:

Name(s)	Enforced state
HID / HSD	neutral histidine, proton on ND1 (δ)
HIE / HSE	neutral histidine, proton on NE2 (ε)
HIP / HSP	+1 imidazolium (both ring N protonated)
ASH / GLH	neutral (protonated) Asp / Glu
LYN / ARN	neutral Lys / Arg
CYM	thiolate (−1)
TYM	tyrosinate (−1)
CYX	disulfide cysteine — S–S bonded to its nearest partner, left unprotonated

Pass --no-honor-protonation to discard the input distinction and protonate every residue purely from --ph instead. Structures using only standard residue names are unaffected either way.

Python API

Import the functions directly from protonate_utils. There are symmetric in-memory and file-to-file entry points for both ligands and proteins.

	Ligands	Proteins
In-memory core	protonate_molecule(mol, ph)	protonate_structure(structure, …)
Convenience	protonate_smiles_string(smiles, ph)	N/A
File → file	protonate_ligands(in, out, ph)	prepare_structure(in, res, out, …)
I/O helpers	read_molecules(path), make_writer(path)	(Biotite PDBFile)

Ligands

Protonate a single SMILES string and get a SMILES string back:

from protonate_utils import protonate_smiles_string

protonate_smiles_string("CC(=O)O")             # 'CC(=O)[O-]'
protonate_smiles_string("OP(=O)(O)O", ph=7.4)  # 'O=P([O-])([O-])O'

protonate_smiles_string raises ValueError on an unparseable SMILES; other failures (e.g. Dimorphite-DL cannot handle the molecule) propagate as exceptions.

Protonate an RDKit Mol while preserving its 3D coordinates:

from rdkit import Chem
from protonate_utils import protonate_molecule, read_molecules

mol = next(read_molecules("ligand.sdf"))
protonated = protonate_molecule(mol, ph=7.4)   # Mol with explicit Hs + coords

Pass add_coord_hs=False to keep protonation implicit (no explicit hydrogen atoms added), appropriate when you intend to serialize to SMILES.

Batch-convert a whole file (the CLI ligand path):

from protonate_utils import protonate_ligands

protonate_ligands("input.sdf", "output.sdf", ph=7.4)

Proteins

Protonate an in-memory Biotite AtomArray and get a hydrogenated one back:

import biotite.structure.io.pdb as pdb
from protonate_utils import protonate_structure

structure = pdb.PDBFile.read("input.pdb").get_structure(model=1)
hydrogenated = protonate_structure(
    structure,
    ligand_res_name="AP5",     # or None / "none" to keep all atoms
    ph=7.0,
    relax=True,
    honor_protonation=True,    # preserve HID/HIE/HIP/CYX/... states (default)
)

protonate_structure raises ValueError if ligand_res_name is given but no atoms with that residue name exist. The returned AtomArray has hydrogens added and reordered to follow their bonded heavy atoms. Force-field protonation residue names are honored by default; pass honor_protonation=False to protonate from ph alone (see Force-field protonation/tautomer residue names).

Read a PDB, protonate, and write a PDB in one call (the CLI protein path):

from protonate_utils import prepare_structure

prepare_structure("input.pdb", "AP5", "output.pdb", ph=7.0, relax=True)

How it works

Ligand protonation

1. Pre-existing hydrogens are stripped; any 3D conformer on the heavy atoms is kept.
2. Dimorphite-DL enumerates candidate microstate(s) within a ±0.5 pH window. One is chosen deterministically by a site-by-site plausibility check rather than by net charge (see Correcting Dimorphite-DL microstates below), and any residual implausible ionization is repaired against the input. The SMILES string is a final tiebreak, so re-runs are stable.
3. The chosen template's formal charges and total hydrogen counts are mapped back onto the original atoms via a charge-insensitive substructure match (so -COOH still matches -COO⁻). Carrying the H count, not just the charge, keeps the RDKit's kekulization correct on aromatic heterocycles.
4. For SDF output, Chem.AddHs adds explicit hydrogens; with 3D input they are positioned from the existing geometry (addCoords=True) and heavy-atom coordinates are never moved, while SMILES input (no coordinates) still gets explicit hydrogens, just without positions. For SMILES output, protonation stays implicit so the writer renders it cleanly.

Correcting Dimorphite-DL microstates

Dimorphite-DL enumerates every microstate whose modeled pKa falls anywhere near the pH window, including many that are negligibly populated at pH 7.4. Left to a "most ionized" or "closest net charge" rule, the selector picks chemically wrong states: it deprotonates amides and phenols and protonates anilines. We add a per-atom legitimacy check (_charge_change_is_legitimate) that compares each candidate to the input atom-by-atom and accepts a formal-charge change only when that group genuinely ionizes near physiological pH:

Group	Typical pKa	At pH 7.4	Dimorphite enumerates	We
Aliphatic amine	pKaH ~10	cation	both	protonate
Amidine / guanidine	pKaH ~12–13	cation	both	protonate
Carboxylic acid	~4	anion	anion	deprotonate
Sulfonic / sulfinic / phosphate / phosphonate	<2–7	anion	anion	deprotonate
Sulfonamide / acylsulfonamide / tetrazole	~3–10	anion	both	deprotonate
Carboxamide N–H	~17–22	neutral	both → [N⁻] or [NH⁺]	keep neutral
Aniline / amino-heteroarene	pKaH ~3–5	neutral	both → [NH⁺]	keep neutral
Cyanamide (N–C≡N)	pKaH ~0	neutral	both → [NH⁺]	keep neutral
Imidazole / pyrazole / indazole / indole / triazole N–H	~10–17	neutral	both → [n⁻]	keep neutral
Phenol / alcohol	~10–16	neutral	both → [O⁻]	keep neutral
Plain thiol / thione	~7–10	neutral	both → [S⁻]	keep neutral

Two further safeguards:

• Repair fallback. When Dimorphite offers only an implausibly-ionized microstate (e.g. it returns just the [N⁻] form of an O-alkyl hydroxamate or imide, with no neutral alternative to select), the offending site is reverted to the input's protonation rather than emitted as-is.
• Input charges preserved. A change is only judged relative to the input, so charges already present in the SMILES (quaternary ammonium salts, N-oxides, mesoionic zwitterions) are never altered.

Borderline acids/bases whose pKa sits right at 7.4 (e.g. p-nitrophenol ~7.15, mercaptoazoles ~7) are deliberately defaulted to neutral; they are ~50/50 at physiological pH, so this is at least as defensible as ionizing them and avoids mis-ionizing the far more common ordinary phenols and amides. Validated across the 2,173-molecule Biogen logS set: no skips, no heavy-atom changes, and the selection is deterministic.

Protein protonation

1. Optionally remove one or more ligands by residue name (a comma-delimited list removes several at once), then strip any existing hydrogens.
2. Normalize force-field protonation/tautomer residue names (HID/HIE/HIP, CYX, ASH, …) to their canonical CCD codes so their bonds can be assigned.
3. Assign covalent bonds from CCD residue templates (connect_via_residue_names).
4. Estimate per-residue formal charges for canonical amino acids at the requested pH (hydride.estimate_amino_acid_charges). When honor_protonation is set (the default), each force-field-named residue is then pinned to the exact state its name encodes — overriding charges, and, for the HID δ-tautomer, swapping the imidazole ring double bond so the proton lands on ND1, and, for CYX, adding the S–S bond.
5. Add hydrogens with Hydride and, by default, relax their geometry.
6. Reorder atoms so each hydrogen immediately follows the heavy atom it is bonded to.

Acknowledgments

Thanks to Manish Sud for helpful suggestions.