Following work by Ree et at., we present a fully automated QM-based workflow that for computing C-H pK a values. A given SMILES string undergoes modifications to produce a list of SMILES for each deprotonated C-H bond. We generate min(1+3 * n rot , 20) conformers for each SMILES using RDKit (v. 2022.09.4), where (n rot ) represents the number of rotatable bonds. Each conformer undergoes optimization in dimethylsulfoxide (DMSO, ϵ = 47.2) using GFN-FF-xTB and the analytical linearized Poisson-Boltzmann (ALPB) as the implicit solvation model. We then remove conformers with relative energies above 3 kcal/mol and select unique conformers by taking the centroids of a butina clustering using pairwise heavyatom root mean square deviation (RMSD) with a threshold of 0.5 Å. Subsequently, we re-optimize the remaining conformers in DMSO with GFN2-xTB and the ALPB implicit solvation model to identify the lowest energy conformer. We then conduct re-optimization in ORCA (v. 5.0.4) 20 , 21 using the dispersion D4 corrected DFT functional CAM-B3LYP, the Karlsruhe 24,25 triple-ζ basis set, def2-TZVPPD, and the Conductorlike Polarizable Continuum Model (CPCM) as the implicit solvation model. CAM-B3LYP is chosen as the optimal functional based on a benchmark study that evaluates the accuracy of different levels of theory, ranging from semiempirical methods (xTB), composite electronic structure methods (r 2 SCAN-3c) to DTF methods (CAM-B3LYP). All these methods are evaluated as single-point calculations or optimization and frequency calculations. For comprehensive details, refer to the SI Benchmark study -computational methods section. Hereafter, we check the geometries for imaginary frequencies and use the total thermal energy at 298.15 K. Following the approach of the Grzybowski group, we compute the heterolytic dissociation energy through the direct deprotonation reaction, AH solv ⇀ ↽ A - solv ; see Equation .