Here, we assess the ability of NEQ-based free energy calculations to resolve the pK a values of 144 residues across 13 proteins. Although largescale studies on the application of NEQ alchemical calculations for predicting mutagenic folding-free energy changes and relative and absolute ligandbinding affinities already exist, such an extension to protein pK a values was heretofore absent from the literature. A seamless free energy workflow that can probe the role of protonation on ligand binding -particularly relevant at an enzymatic active site -and resolve the underlying pK a values of both individual residues and bound molecules is highly desirable. Here, we take a step towards that goal. Although (de)protonation is the smallest topological change that a residue can undergo, it nevertheless results in a significant charge shift. We find that such perturbations and the corresponding free energies can be readily resolved using our pmx-based approach (i.e., AUE: 0.68 ± 0.05 pK), with accuracy comparable to FEP+ and demonstrate the ability of this approach to resolve the pK a of coupled residues.