pK_a Prediction

Knowing the pK_a of a molecule is key to understanding its structure and reactivity. Rowan's pK_a prediction workflow uses machine-learned interatomic potentials and semiempirical solvation energies to enable fast and accurate prediction of pK_a values with minimal empiricism.

Two quantum mechanics–style methods are available:

• Rowan pKa (Wagen 2026), a g-xTB-powered workflow (our recommended default), and
• Rowan pKa (Wagen 2024), our original AIMNet2-powered method.

At a high level, the Rowan pK_a prediction workflow works much like other quantum chemistry–based pK_a prediction workflows, but replaces the slow DFT calculations in these workflows with g-xTB or the AIMNet2 neural network potential (trained on ωB97M-D3BJ/def2-TZVPP data), leading to multiple orders-of-magnitude speedup over DFT. g-xTB performs well on thermochemical benchmarks and can handle organometallic systems that are out of reach for AIMNet2. Since AIMNet2 and g-xTB do not take solvation into account, ∆G_solv is computed for both from a single-point GFN2-xTB calculation with the CPCM-X implicit water model.

An example API script is available on GitHub. For fast, high-throughput macroscopic pK_a prediction, see Rowan's macroscopic pK_a workflow instead.

Protonation, deprotonation, and conformer generation

The workflow begins by iteratively adding or removing hydrogens to the molecule of interest using RDKit and quickly estimating the proton affinity for the conjugate acid/conjugate base pair using single-point energy calculations. If (de)protonation is predicted to yield a pK_a value close to the desired range, a conformer ensemble is generated using openconf.

The g-xTB-based method then reduces this set of conformers using the solvent-aware clustering and filtering from Rowan's solvent-dependent conformer search workflow. The retained structures are optimized with g-xTB in implicit solvent, receive g-xTB vibrational free-energy corrections from single-point Hessian calculations, and are scored with the GFN2-xTB/CPCM-X water solvation correction.

In the AIMNet2-based method, the retained structures are optimized with AIMNet2, receive AIMNet2 vibrational free-energy corrections, and are scored with the GFN2-xTB/CPCM-X water solvation correction.

In both, conformers transferred from the neutral parent are reused where possible so that related protonation states are compared consistently.

ΔGs to pK_as

Following Pracht and Grimme, we combine energies from each conformer via Boltzmann averaging to generate a single hybrid ∆G for the conjugate acid/conjugate base pair. The two methods diverge in how ∆G is converted into a pK_a value.

The AIMNet2-based method applies a quadratic free-energy relationship, plus an element- and valence-specific constant for each conjugate acid/conjugate base pair to correct for systematic errors in thermochemistry and solvation. (The correction for divalent oxygen, i.e. deprotonation of ROH, is defined as zero to eliminate extra degrees of freedom.) The final relationship contains 7 empirical parameters, shown here:

Because it relies on AIMNet2, this method is restricted to molecules within AIMNet2's supported elements, but remains a fast option for organic, drug-like molecules.

The g-xTB-based method instead uses per–functional-group linear free-energy relationships. Scaling factors were fit to a subset of the Dwar-iBond dataset, and functional groups were added automatically based on the Bayesian information criterion.

Modes

Rowan pK_a can be run in three modes: careful, rapid, or reckless, which tune the search breadth. For the AIMNet2 method, which uses openconf, the modes correspond to:

In general, we've found that the "careful" mode offers the best combination of speed and accuracy, but faster modes offer increased performance where high throughput is required.

Accuracy

On a variety of benchmarks (also described in the preprint), Rowan pK_a displays mean absolute errors of around 1 pK_a unit, and decent predictive accuracy as assessed by Kendall's τ and R² values. Here's the AIMNet2 method's performance on the SAMPL7 benchmark set:

For a full list of all assays surveyed, see the preprint.

We reinvestigated a variety of benchmark systems with the new g-xTB method after development was finished (without using the benchmarks to tune any scaling factors), and found improved accuracy across the board. On the Rombouts set of tricyclic amidine BACE1 inhibitors, for instance, g-xTB shrinks the MAE from 1.06 to 0.59 pK_a units.

We also see improvements on the Miller–Doukas–Seydel folate-inhibitor dataset:

And on the Müller dataset examining the effect of α- and β-substituents (e.g. oxetanes) on cyclic amine basicity: