From the titration of protonatable group, one can read the so-called pK which is equal to the pH value where the group is half-protonated (i.e. when 50% such groups would be protonated). The pK is equal to the Henderson–Hasselbalch pK (pK) if the titration curve follows the Henderson–Hasselbalch equation. Most pK calculation methods silently assume that all titration curves are Henderson–Hasselbalch shaped, and pK values in pK calculation programs are therefore often determined in this way. In the general case of multiple interacting protonatable sites, the pK value is not thermodynamically meaningful. In contrast, the Henderson–Hasselbalch pK value can be computed from the protonation free energy via
and is thus in turn related to the protonation free energy of the site via
The protonation free energy can in principle be computed from the protonation probability of the group (pH) which can be read from its titration curve
Titration curves can be computed within a continuum electrostatics approach with formally exact but more elaborate analytical or Monte Carlo (MC) methods, or inexact but fast approximate methods. MC methods that have been used to compute titration curves are Metropolis MC or Wang–Landau MC. Approximate methods that use a mean-field approach for computing titration curves are the Tanford–Roxby method and hybrids of this method that combine an exact statistical mechanics treatment within clusters of strongly interacting sites with a mean-field treatment of intercluster interactions.
In practice, it can be difficult to obtain statistically converged and accurate protonation free energies from titration curves if  is close to a value of 1 or 0. In this case, one can use various free energy calculation methods to obtain the protonation free energy such as biased Metropolis MC, free-energy perturbation, thermodynamic integration, the non-equilibrium work method or the Bennett acceptance ratio method.
Note that the pK value does in general depend on the pH value.
This dependence is small for weakly interacting groups like well solvated amino acid side chains on the protein surface, but can be large for strongly interacting groups like those buried in enzyme active sites or integral membrane proteins.
While many protein pKa prediction methods are available, their accuracies often differ significantly due to subtle and often drastic differences in strategy.