Our 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. If (de)protonation is predicted to yield a pKa value close to the desired range, conformers are generated using the ETKDG algorithm and optimized using the MMFF forcefield. After removing redundant geometries, low-energy conformers are then optimized using the GFN2-xTB semiempirical method, single-point energies are calculated using the AIMNet2 model trained on ωB97M-D3/def2-TZVPP (henceforth abbreviated as AN2-ωB97MD3), and the lowest N undergo full optimization and frequency calculations at the AN2-ωB97MD3 level of theory to generate a free energy for each conformer, including all standard entropy-and symmetry-based corrections. (N is set to 10 by default but can be varied by the end user; vide infra.) Since AIMNet2 does not take solvation into account, ∆Gsolv is computed from a single-point GFN2-xtb calculation with the CPCM-X implicit water model. Following Pracht and Grimme, we combine energies from each conformer to generate a single hybrid ∆G for the conjugate acid/conjugate base pair and apply a quadratic free-energy relationship to convert this into pKa values. Initial application of this strategy provided disappointing results, with considerable systematic error observed between different classes of acids. Functional group-dependent systematic error is common in ab initio pKa prediction, and many different corrections have been employed: Schrödinger's Jaguar pKa contains a different linear free-energy relationship for roughly a hundred different functional groups, while Jensen and co-workers employed different reference molecules for each class of compound under study. To minimize the number of empirical parameters employed in our workflow, we instead add an element-and valence-specific constant to ∆G for each conjugate acid/conjugate base pair. (The correction for divalent oxygen, i.e. deprotonation of ROH, is defined as zero to eliminate extra degrees of freedom.) These constants are determined by Levenberg-Marquardt least-squares minimization, not by comparison to any individual reference compound, and can be thought of as a correction for systematic errors in both thermochemistry and solvation. Similar corrections are often used to account for dispersion in mean-field electronic structure theory methods without substantially compromising the generality of the method. Workflow development and parameter optimization was conducted on a slightly modified version of the TR224 dataset, also used by Pracht and Grimme. The final Rowan pKa workflow contains 7 tunable parameters, shown in Table .