is evaluated with the partial charges of λ-group k in the deprotonated state (q B i ) and the same interpolated charges for all other λ-groups. Thus, 2N sites computations are needed to evaluate the gradients for all titratable sites. The same arguments apply to pair-wise Lennard-Jones interactions, but because the contribution of Lennard-Jones interaction to pK a shift is minor, we neglected them in this work (see Supporting Information).