One of the most direct benchmarks for electrostatic models of macromolecules is provided by the pK(a)’s of ionizable groups in proteins. Obtaining accurate results for such a benchmark presents a major challenge. Microscopic models involve very large opposing contributions and suffer from convergence problems. Continuum models that consider the protein permanent dipoles as a part of the dielectric constant cannot reproduce the correct self-energy. Continuum models that treat the local environment in a semi-microscopic way do not take into account consistently the protein relaxation during the charging process. This work describes calculations of pK(a)’ s in protein in an accurate yet consistent way, using the semi-microscopic version of the protein dipoles Langevin dipoles (PDLD) model, which treats the protein relaxation in the microscopic framework of the linear response approximation. This approach allows one to take into account the protein structural reorganization during formation of charges, thus reducing the problems with the use of the so-called "protein dielectric constant", epsilon(p).