MDP Templates: Why use fep-lambdas?

[a-ws-m]

I've been adapting my MDPOW example for the specific molecule I'm interested in, which is quite a bit larger, and instabilities in the FEP simulations have caused me to review the way the MDP templates work. I've noticed that the fep-lambdas key is used. From the docs:

fep-lambdas is different from the other -lambdas keywords because all components of the lambda vector that are not specified will use fep-lambdas.

Therefore, using this parameter couples everything, not just VDW or Coulomb components.

I understand that the couple-lambda parameters are used to turn off charges in the case of VDW coupling, so the two cases are distinct. However, is it necessary to use fep-lambdas, rather than vdw-lambdas or coul-lambdas?

[VOD555]

In short, when decoupling the interactions of a specific type of moelcules, using the couple-moltype option with couple-lambda0 and couple-lambda1, there's no difference between using fep-lambdas and vdw-lambdas or coul-lambdas.

For example, given couple-lambda0="vdw-q" and couple-lambda1=vdw, changing fep-lambdas from 0 to 1 will only decouple the colomb interactions on the molecule and won't change the VDW interactions. For more details about these options, please look at https://manual.gromacs.org/documentation/current/user-guide/mdp-options.html#mdp-couple-moltype .

One advantage of using couple-moltype and fep-lambdas is that one doesn't need to prepare two sets of parameters for the molecule before and after decoupling. In MDPOW, we modify the template mdp file with these options (see here).

[a-ws-m]

Thank you for the swift reply. I had noticed the difference in the couple-lambda parameters and I believed what you just described to be the case until recently. I can see that setting couple-moltype allows VDW/Coulomb interactions to be changed, and that using the correct values of couple-lambda will determine which components are affected by lambda.

However, the documentation does not specify that setting couple-moltype will disable mass, temperature, restraint or bonded parameter coupling. Furthermore, I have found that just by changing fep-lambdas to vdw-lambdas or coul-lambdas, my simulations become stable. This suggests a change in behaviour.

I will investigate a little bit more myself and see if changing the parameters in a simulation that is stable either way leads to different free energy calculations.

[a-ws-m]

I did a quick experiment: I changed the code to use vdw-lambdas or coul-lambdas, rather than fep-lambdas. The calculated solvation energies and resulting partition coefficient for my Martini benzene example notebook changed significantly: both solvation energies differed by more than the error estimate and there was a final difference of 0.14 in the final partition coefficient. You can find these results in this branch.

[VOD555]

I've run the example notebook and got the logP of benzene to be 2.25, while the original one is 2.34, with a difference of 0.4 kJ/mol in the solvation free energy in octanol. I'd say though the system is simple, it is still not easy to get a converged estimate for logP from a single run of fep simulations. Can you run another 4 runs for using fep-lambdas and vdw-lambdas or coul-lambdas to see some statistics?

It is obvious that these options don't affect the performance when decoupling coulomb interactions, as in both cases the coulomb part of the solvation free energies is 0 in either octanol or water. Can you try another run in which the vdw and coulomb interactions are turned off at the same time (set vdw-lambdas = coul-lambdas, and couple_lambda0="vdw", couple_lambda1="none")?

Any opinions? @orbeckst

[VOD555]

And it's always good to post the question on gromacs forum, there are people who are much more familiar with gromacs than Oliver and I.

[orbeckst]

I don't know if there is a functional difference between fep-lambdas and could-lambdas/vdw-lambdas. I don't think that it should affect stability but I don't know, the GROMACS free energy code is a complicated beast. What do the GROMACS docs say?

Otherwise I agree with @VOD555 that converged FEP calculations are necessary before drawing any conclusions as to the actual effect (assuming that it doesn't trigger a bug of some kind).

[orbeckst]

(Octanol can take forever to converge!!)

[a-ws-m]

Thank both for your input. I may have a chance to look at running some longer simulations but I don't have access to my usual machine at the moment. After looking around, I think I might have identified the reason my simulations were unstable.

My topology file contains these lines:

[constraints]
  1 2   1      0.261     17000

The final term is an artifact from some other software that I used; usually, you can copy>paste this line between bonds and constraints without any change in functionality, seeing as the force constant is ignored. However, when free-energy is turned on, I believe this acts as an end state. I think the path between constraints may be defined by the bonded-lambdas, which is why using vdw- or coul- explicitly negates the issue. Hopefully I can check in more detail next week and get back to you.

In any case, maybe it's better to use vdw- or coul- explicitly? The docs suggest you should "fill in all arrays explicitly, just to make sure things are properly assigned."