We have developed an efficient and reliable protocol for the calculation of pK(a) values in aqueous solution from density functional calculations. We establish a standard linear regression fit using only calculated energies of deprotonation and experimental pK(a) values; all other factors, including most entropic effects, are absorbed into the fitting constants. In this article we fit a small training set of 34 experimentally well-characterized molecules to determine the best level of theory among those tested (i.e., the optimum compromise between efficiency and accuracy for the basis set, the exchange-correlation functional, the (continuum) solvation model and the level of geometry optimization). Our main findings are that a relatively modest basis set (6-311+G**) suffices for the calculation of the energy differences, with an even small basis set (3-21G*) sufficient for the preceding geometry optimization. Using a solvation model (COSMO in our case) throughout is essential to achieve reliable results. The exchange-correlation functional plays only a modest role; in particular, pure DFT functionals that allow the efficient calculation of the Coulomb term are perfectly adequate. The final protocol will be applied subsequently to data sets much larger than commonly used in such studies.