We describe the development and application of a computational method for the prediction and rationalization of pK(a) values of ionizable residues in proteins, based on ab initio quantum mechanics (QM) and the effective fragment potential (EFPs) method (a hybrid QM/MM method). The theoretical developments include (1) a covalent boundary method based on frozen localized orbitals, (2) divide-and-conquer methods for the ab initio computation of protein EFPs consisting of multipoles up to octupoles and dipole polarizability tensors, (3) a method for computing vibrational free energies for a localized molecular region, and (4) solutions of the polarized continuum model of bulk solvation equations for protein-sized systems. The QM-based pK(a) prediction method is one of the most accurate methods currently available and can be used in cases where other pK(a) prediction methods fail. Preliminary analysis of the computed results indicate that many pK(a) values (1) are primarily determined by hydrogen bonds rather than long-range charge-charge interactions and (2) are relatively insensitive to large-scale dynamical fluctuations of the protein structure.