Cutoff-dependent deviations in free energy for force switching and potential shifting

My student Lindsey Whitmore and I have been looking at the different vdw-modifiers and (after finding a bug in the force-switch modifier that Berk impressively rapidly patched), found that there is still fundamentally a cutoff dependent term in the force-switch and potential-shift modifiers for free energies, whereas potential-switch gives essentially cutoff-independent results.

We suspect that this is occurring in most if not all MD codes - the lammps potential-shift van der Waals potential appears to have similar behavior (word of mouth from another user).

We wrote this up here: [2410.14187] Force switching and potential shifting lead to errors in free energies of alchemical transformations, see Fig 1 and 2.

The basic issue is that the constant energy subtracted out at the cutoff (force-switch has some energy subtracted at the cutoff as well) is lambda-dependent, but this is not currently included anywhere in the dH/dl or the dU(lambda)-dU(Delta lambda). I think it’s possible but ugly to include terms to calculate it, but I have not finished the entire calculations (for our free energy calculations in general, we have just always been using potential-switch).

Let me know what would make the most sense here as next steps. There probably should be warning with non-potential-switch vdw-modifiers if they are not tweaked to give correct results.

I don’t understand. The lambda dependence is handled fully by GROMACS.

If you shift the potential, then the free-energy difference is going to be dependent on the cut-off. This will not be the case with potential-switch. Is that the effect that you are seeing? This is an expected, but very small, effect.

If you shift the potential, then the free-energy difference is going to be dependent on the cut-off. This will not be the case with potential-switch. Is that the effect that you are seeing?

That is my understanding. It’s more of a medium size or somewhat small effect, though, not very small. Because you are shifting the energy in a way that affects all pairs of particles that are closer together, so it ends up getting multiplied by a fair number of interactions. But one should be able to include this as well, though I will need to sit down and get the exact terms.

But that is exactly what the force field specifies. If you want to take these terms into account, you need to add dispersion correct. GROMACS handles the lambda dependence of that correctly as well.