Predicting how aqueous solvent modulates the conformational transitions and influences the pKa values that regulate the biological functions of biomolecules remains an unsolved challenge. To address this problem, we developed FDPB MF, a rotamer repacking method that exhaustively samples side chain conformational space and rigorously calculates multibody protein-solvent interactions. FDPB MF predicts the effects on pKa values of various solvent exposures, large ionic strength variations, strong energetic couplings, structural reorganizations and sequence mutations. The method achieves high accuracy, with root mean square deviations within 0.3 pH unit of the experimental values measured for turkey ovomucoid third domain, hen lysozyme, Bacillus circulans xylanase, and human and Escherichia coli thioredoxins. FDPB MF provides a faithful, quantitative assessment of electrostatic interactions in biological macromolecules.conformational flexibility ͉ electrostatics ͉ pKa prediction ͉ protein modeling B y strongly interacting with charges, the solvent significantly modulates the electrostatic properties of biomolecular systems. Although variations in pH and ionic strength are responsible for physiologically important conformational transitions (1), quantitatively assessing their role in modulating electrostatic energies is particularly challenging (2). Approaches in which solvent is modeled as a continuum polarizable medium while biomolecules are modeled in explicit molecular detail have become widely used in recent years (3, 4). However, no current method combines broad conformational sampling with a rigorous solvation model that can predict quantitatively and efficiently the effects of solvent on protein energetics and conformations.In principle, molecular dynamics and free-energy simulations monitoring protonation events can model accurately the energetic effects of solvation and conformational relaxation of biomolecules. However, despite recent progress, these techniques often fail to converge in a reasonable amount of computer time and therefore do not provide reliable predictions of thermodynamic properties (5, 6).Rotamer repacking methods sample more efficiently and exhaustively the conformational space of biomolecules. However, these methods require the energy function be decomposed into self energies and interaction energies between pairs of residues (7-10). Consequently, repacking methods cannot model the effects of conformational relaxations involving more than two positions at a time. Many important properties of biomolecular surfaces (i.e., their interface with the solvent and the distribution of bulk ions around them) are not pairwise factorable and depend on the simultaneous knowledge of the conformations of all residues. Approximating these effects by relaxing the protein-solvent interface only at the vicinity of solventexposed charged residues does not improve the prediction of pKa values in proteins (11). Current rotamer repacking methods have also difficulties in assessing accurately and reliably the ...