Geometries of all structures were optimized using the dispersion-corrected 36,37 B3LYP-D3(BJ) functional with the 6-31+G(d) basis set using the Gaussian 16 program 35 (Fig. ). Vibrational frequency calculations were performed at the same level of theory of the geometry optimization to confirm that each structure is a local minimum (i.e., with no imaginary frequencies). Single point energy calculations were carried out using the M06-2X functional with the 6-31+G(d) basis set. Solvation energy corrections were calculated using the SMD solvation model 41 in single point energy calculations with water as the solvent. Carboxylic acids (ArHet-CO2H) and carboxylate anions (ArHet-CO2 -) may have several conformers depending on whether the carboxylic acid or carboxylate group is co-planar with the heteroaromatic ring. The "SetDihedralDeg" function in RDKit was used to generate conformers of carboxylic acids and carboxylate anions by rotating about the Cipso-Ccarbonyl bond. Only the lowest energy conformer of each structure was used to compute the reported properties. The Automated Quantum Mechanical Environments (AQME) software was used in post-processing to check for self-consistent field (SCF) and geometry optimization convergence errors and imaginary frequencies. Calculations with convergence errors were resubmitted by using the intermediate structure during the previous geometry optimization with the lowest root-mean-square gradient as the input geometry. In cases where imaginary frequencies were present, the calculations were adjusted by slightly perturbing the geometry and resubmitted with the keyword "opt=(calcfc,maxstep=5)". This automated process was repeated twice, and any calculations still showing errors after the attempted re-calculations were not included in the final database.