Lambdas in FEP calculations

Dear all,

I am slightly confused about lambda selection in FEP calculations.

When it comes to absolute FE calculations, the electrostatic interactions are switched off first, and then vdw interactions are switched off (just an example):

; lambda = 0 1 2 3 4
coul-lambdas = 0.0 0.5 1.0 1.0 1.0
vdw-lambdas = 0.0 0.0 0.0 0.5 1.0

How one should choose the lambda values in relative (perturbation) calculations?

If I use the same approach, then I would end up with situation where I have switched off coulombic interactions for several atoms, and switched them on for several other atoms; all this while vdw interactions still have not been changed. That is, at lambda 2 I would already have coulombic interactions from state B, whereas vdw interactions would still be from state A.

Could it really be as straight forward as:

; lambda = 0 1 2 3 4
fep-lambdas = 0.0 0.25 0.5 0.75 1.0

Or am I missing something?


Your wording is confusing.

lambda always moves between state A and B, independently of what you want to calculate in the end. You don’t specify what the A and B states of your topology are for your second case, so I can’t say what is right or wrong.

I am sorry for confusing problem formulation.
I will try to explain it better:

In many of the free energy calculation tutorials, the focus is usually on the solvation free energy or absolute binding free energy. In these cases the state A is molecule of interest in the solution or bound to the protein, whereas the state B is annihilated molecule where all interactions are switched OFF (basically, molecule is “switched off”). Here, moving from state A to B switches OFF various interaction types subsequently, and it makes sense to me, as there are no intermediate states where some neighboring atoms will have various/different interaction types switched ON/OFF.

In free energy perturbation calculations, lets say, I want to transform methyl group in a molecule (state A) to methoxy group (state B). At the state lambda 2 in previously mentioned lambda array, I would have a situation where charges of the metyl group (state A) would be changed to the charges of the methoxy group (state B); while the vdw interactions would still correspond to initial state A (since dw-lambdas is still 0.0 at this point). Thus, there would be a situation where methoxy group would have coulombic charges ON while vdw interactions would still be OFF (which would, probably, lead to system blowing up).

Besides, the functional group that we are changing is connected to main molecule where all the interactions are ON throughout perturbation process.

Topology example:

[ atoms ]
; nr type resno resnm atom cgnr charge mass typeB chargeB massB
… atoms that does not change …
40 c 1 LIG C01 40 0.100000 1.0080 du 0.000000 1.0000
41 h 1 LIG H01 41 0.130000 1.0080 du 0.000000 1.0000
42 h 1 LIG H02 42 0.130000 1.0080 du 0.000000 1.0000
43 h 1 LIG H03 43 0.130000 1.0080 du 0.000000 1.0000
44 du 1 LIG DU44 44 0.000000 1.0000 o -0.300000 16.0000
45 du 1 LIG DU45 45 0.000000 1.0000 c 0.300000 12.0100
46 du 1 LIG DU46 46 0.000000 1.0000 h 0.130000 1.0080
47 du 1 LIG DU47 47 0.000000 1.0000 h 0.130000 1.0080
48 du 1 LIG DU48 48 0.000000 1.0000 h 0.130000 1.0080

I hope this helps to explain my concern a little bit better.