[gmx-users] Binding free energy calculation by thermodynamic integration

[Moutusi Manna]

Hi,
I want to use thermodynamic integration for calculating free energy of binding of ligand to its receptor protein. For that I would like to follow the two step approach: first decoupling Coulombic and then vdW interactions. 

1. Shall I start with the same coordinate and topology files for both decoupling process and only vary the lambda values? Or before turning off vdW terms all charges in the topology file should set to zero?  

2. The MDP files I am going to use are:
MD1.mdp
free_energy = yes
init_lambda = 0.2 ; different values of lambda 0.0,0.1,0.2,0.3,....,1
delta_lambda = 0
foreign_lambda = 0.1 0.3
sc-alpha = 0
sc-sigma = 0
couple-moltype = Ligand
couple-intramol = no
couple-lambda0 = vdwq
couple-lambda1 = vdw

MD2.mdp
free_energy = yes
init_lambda = 0.2 
delta_lambda = 0
foreign_lambda = 0.1 0.3
sc-alpha = 0.5
sc-power = 1.0
sc-sigma = 0.3
couple-moltype = Ligand
couple-intramol = no
couple-lambda0 = vdw
couple-lambda1 = none

foreign_lambda: I put only neighboring lambda points. Is that OK? Or I should put all 11 lambda points 0 0.1 0.2 .... 0.9 1.0 ? If at any stage I have to increase the number of intermediate lambda states, say at 0.05 interval- then the closest lambda values will be different (like, for init_lambda = 0.2 foreign_lambda = 0.15 0.25). In that case do I have to rerun all simulations with their new foreign_lambda values?

Soft-core parameters are set to zero during charge removal. Is it OK? 

3. During decoupling process I want to specify the position of ligand relative to protein. For which I use pull code. But how to correct its contribution to free energy of binding?
pull =  umbrella
pull_group0 = Ligand
pull_group1 = Protein
pull_rate1 =  0.0000
pull_k1 = 750

Thanks and best regards,
Moutusi