For each set of λ values, called a grid point, a 10 ns MD simulation was performed. ∂V /∂λ i were saved every ps, which is approximately equal to the auto-correlation times for the λ-coordinates. The total charge of the system was kept neutral by simultaneously changing the charge of a single buffer particle, discussed below. The ∂V /∂λ i values were averaged over the last 9 ns of the trajectories. To obtain an analytical expression for V MM a fifth order polynomial was fitted to these averages for Asp and Glu, while an eighth order polynomial was fitted for His, taking into account possible linear dependencies of the coefficients (see Supplementary Information). Fitting errors were below 0.5 kJ/mol for Asp and Glu and below 1 kJ/mol for His, which are of similar magnitude as the statistical accuracy of the derivatives.