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

Developers discussions
,
1

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.

2

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.

3

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.

4

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.

5

Coming back to this. It’s clear from the data shown that the free energy is cutoff-dependent for force-switch and potential-shift, and it’s only cutoff-independent for potential-switch and LJ-PME (data for that not shown, but will be in the next release of the preprint/submission). For anthracene, for standard cutoffs distances, the difference is 0.5 kJ or so, enough to matter for reproduciblity (see preprint cited above)

For potential-switch, all of the missing energy is at g(r) > rvdw-switch, and so handling the lambda-dependence of the switch correctly results in free energy differences that are independent of the cutoff. However with force-switch or potential shift, there is a lambda-dependent contribution to the free energy at all g(r). If you are at lambda=0.5 and are calculating the potential energy at lambda=0.6, then you need to take into account that the shift in the potential at all r, equal to the value of the unmodified potential at the cutoff at lambda=0.6, is different than the shift in the potential at lambda=0.5, equal to the value of the unmodified potential at the cutoff at lambda=0.5.

Or more algorithmically, when looping over the perturbed interactions at “foreign” lambdas to determine U(lambda+dlambda)-(lambda), as far as I can tell, it uses the cutoff-dependent shift in the potential at the CURRENT lambda, rather than the cutoff-dependent shift in the potential at that foreign lambda. This contribution is also neglected in dhdl.

I spent a bit of time in gmxlib/nonbonded/nb_free_energy.cpp today, but haven’t quite figured out how this could be implemented.

6

If I understand correctly, when looping over foreign lambdas, the gmx_nb_free_energy_kernel takes each foreign lambda as input. Therefore, it computes the potential shift at each foreign lambda rather than at the current lambda. The code I looked at is in nbnxm/freeenergydispatch.cpp – dispatchFreeEnergyKernel

7

Ah, annoying issues. This should be relatively straightforward to fix, but complicates the code even more.

We should ensure that our tests catch this.

8

Though, I think that’s what it’s supposed to do? The potential energy is different at the foreign lambda because it uses the cutoff at the foreign lambda. If one uses the potential shift at the current lambda to compute the foreign lambda energy, then one would not get that energy difference.

9

Thinking through this, what might actually be the case is that the free energy of the potential shifted and force switched molecules is DIFFERENT than the free energy of the potential-switched molecule. This is because the fully coupled state is actually a different Hamiltonian - one where there is a shift of U(r) at all r. One could argue that that is indeed the “right” free energy value, but then one finds that the free energy would depend on an arbitrary simulation parameter of the cutoff, which is not great - one really would prefer to have a cutoff-independent result! For example, we have checked the LJ-PME free energies as a function of cutoff, and it is also independent of cutoff - LJ-PME is probably the best reference for the “right” value of the energy.

We are finalizing submission of the paper that was first put out as [2410.14187] Force switching and potential shifting lead to errors in free energies of alchemical transformations (now including the LJ-PME results and also some results from LAMMPS) and and working to provide a definitive statement that either 1) this is a bug that should be fixed, 2) this is a bug that HAS BEEN fixed (which is preferable to 1), or 3) fundamentally, one should use something like LJ-PME or potential-switch for cutoff-independent free energies because of the subtraction of energy at U(r).

10

Of course the free-energy depends on the LJ cut-off. The LJ cut-off is a parameter of the force field, just like all other parameters. It is clear that this is undesirable. One solution is using LJ-PME, if that is deemed compatible with the given force field.

11

Sure, but it’s MUCH better if one can get them aligned; otherwise one often has a substantial cost tradeoff. potential-switch is relatively independent of cutoffs - almost completely independent of cutoffs for solvation in homogeoenous liquids, somewhat but not completely independent for protein-ligand interactions (depends on structure in g(r) outside the cutoff (https://doi.org/10.1021/jp0735987)), whereas potential-shift and force-switch free energies are significantly dependent on cutoff (depend on structure of g(r) INSIDE the cutoff).

So I think the question to resolve is whether the current cutoff-dependent free energies are “working as is mathematically intended, even if unfortunate” or “actually, there’s a cutoff-dependent correction term that is being omitted”). YiqiChen’s comment above suggest the former, that the lambda-dependence of the shift is being included.

12

I don’t understand YiqiChen’s comment. The potential shift is lambda independent. It is e.g. 1/r_c^6 for the dispersion.

13

Well, if YiQiChen’s comment is wrong, and the potential-shift and force-switch terms are not lambda independent, then that is a coding error, as they should be - the potential at r=r_c depends on lambda.

14

Sorry, I was wrong. The dispersion shift and repulsion shift terms are both constants and donot vary with lambda.

15

You mean that they depend on lambda at r=r_c because of the soft-core distance? That effect should be extremely small at the cut-off.

16

At r > 2*sigma or so, then U(r_c) varies roughly linearly with lambda, so it’s not really dependent on the soft-core potential itself.

17

OK, I dug into this a bit more. The question is, is the lambda dependent shift properly accounted for in GROMACS or not?

For potential-shift, the answer appears to be yes. We can do this by checking to make sure that the potential energy differences generated in the .dhdl.xvg file are equivalent to the potential energy differences generated by running two independent simulations of GROMACS at different lambda values, with one of them generating a trajectory and the other analyzing the trajectory with the --rereun flag set with the first trajectory.

In this case, there’s no question that the “foreign lambda” is carried out at the proper cutoff cutoffs, since it is done in a completely independent calculation. I did these calculations in double precision, as it was getting a little noisy in single precision.

What this showed was that yes, the internal Delta U calculation matched the Delta U calculation done by subtracting within 8 or so decimal places. So the lambda dependence of the cutoff (at least the vdW cutoff) appears to be implemented correctly.

This means (as I was coming around to anyway) that the free energy of a potential-shifted (and therefore also, to a lesser extent, force switching) is significantly cutoff-dependent, in a way that a potential switched, or PME-vdw calculation is not, because the potential at r < r_vdw_switch/r_vdw_shift is affected for potential shift/force-switch, but the potential-switched force field is only different from the no-cutoff force field when g(r) > r_vdw_switch - the close-range g(r) is essentially unchanged.

So my conclusion there is no associated bug here, but a “feature” of correct implementation of potential shifting.

My feeling is that there should be a warning with free energy calculations when using something other than vdw-modifier=potential switch or not using LJ-PME . . . would be interesting to see other people’s comments.

It is possible to calculate a correction to this effect, though I don’t like calculating things that are corrections to a “correct” (i.e. internally consistent) simulation. One would be producing a free energy that it would be if the simulation was not potential shifted, not the simulation that was performed.

I’ll test force-switch next.

18

Yeah, force-switch agrees as well, indicating that it’s correctly implemented, but that a force-switch potential has a different potential, dependent on the cutoff and thus a different free energy.

19

I don’t see any incorrectness issues here. If the user uses the force field with the intended cut-off treatmeant, the user gets free-energy differences consistent with that. That the results differ from a no cut-off,e.g. LJ-PME, result is kind off irrelevant and holds for any other property as well. The cut-off effects (should) have been parametrized into the force field when using a cut-off during parametrization. The extreme example of this is lipid force fields, most of which can not be used with LJ-PME as the area per lipid would then be off.

20

Yes, and no. There are a bigger and smaller cutoff effects. But a good point to raise. Let me draft up something complete that will go into the paper, this does need to be explained well.

I think we have reached the point we are pretty sure the code behaves as expected, we just have to be careful to define what “as expected” means.