The electrostatic calculations outlined above provide (free) energies of each of the 2N protonation microstates (10 (link)) in the system, where N is the number of ionizable sites. To make the subsequent calculation of the partition functions (and pK1/2) manageable, a fast variant of a clustering approach is used (34 (link)). The approach subdivides the interacting sites into independent clusters based upon the strength of electrostatic site–site interactions between them. All electrostatic interactions for each ionizable site are sorted from highest to lowest; the top Cmax sites are then selected to contribute to the calculation, and all others are ignored. The partition function for the site is then factored into computationally manageable components of maximum size Cmax. Here, Cmax = 17 is used: in tests on 600 representative proteins (35 (link)), we found (J. Myers, G. Grothaus and A. Onufriev, manuscript submitted) that Cmax = 17 resulted in average errors of 0.2 pK units, compared with a standard treatment based upon a Monte Carlo approach (16 (link)).
The probability of protonation is computed for every site over a range of pH values equally spaced by 0.1 pH units apart. Individual curves can be displayed for user-selected residues, and the total protonation curve is generated, showing the computed isoelectric point of the molecule. A diagram showing the 10 lowest protonation states and their relative free energies is also generated. These diagrams were found useful (10 (link)) for analysis of proton transfer events in biomolecular systems.