All energy terms are computed using solution phase geometries unless noted otherwise. DH f (X) is computed using either PM6-DH+, PM6, PM7, PM3, or AM1, while DG solv (X) is computed with the COSMO solvation method using MOPAC2016. A maximum of 200 optimization cycles were used for geometry optimizations. The conformers for each protonated isomer were generated from SMILES strings using RDKit. The number of conformers is min(1