While the linear scaling of the gradients for the pair-wise potentials in equation 15 in principle is a great improvement over the formal exponential scaling in equation 8, the requirement of performing 2N sites calculations per MD step still poses a computational bottleneck, in particular for larger systems. To overcome this bottleneck for electrostatic interactions, we follow the suggestion by Brooks and co-workers to interpolate charges rather than interaction functions. When interpolating the partial charges between the protonation states of N sites chemically uncoupled titratable sites, the λ-dependent Coulomb energy becomes: