A new method is proposed for constant pH molecular dynamics (MD), employing generalized Born (GB) electrostatics. Protonation states are modeled with different charge sets, and titrating residues sample a Boltzmann distribution of protonation states as the simulation progresses, using Monte Carlo sampling based on GB‐derived energies. The method is applied to four different crystal structures of hen egg‐white lysozyme (HEWL). pKa predictions derived from the simulations have root‐mean‐square (RMS) error of 0.82 relative to experimental values. Similarity of results between the four crystal structures shows the method to be independent of starting crystal structure; this is in contrast to most electrostatics‐only models. A strong correlation between conformation and protonation state is noted and quantitatively analyzed, emphasizing the importance of sampling protonation states in conjunction with dynamics. © 2004 Wiley Periodicals, Inc. J Comput Chem 25: 2038–2048, 2004