Free energy peturbation mdp file settings for solvation free energy calculation

User discussions
1

Dear gromacs users,

I am setting up a solvation free energy calculation for a small molecule in water. I have an equilibrated solute-solvent box which I will use to setup the free energy perturbation calculations.

I am using the following mdp file setting: (for the 21st lambda index)
free_energy = yes
init_lambda_state = 20
delta_lambda = 0
calc_lambda_neighbors = 1 ; only immediate neighboring windows
couple-moltype = Merocyanine ; name of moleculetype to decouple
couple-lambda0 = none ; vdw and coulomb interactions are turned off for lamda=0
couple-lambda1 = vdw-q ; turn on vdw and coulomb interactions for lambda=1
couple-intramol = yes
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
vdw_lambdas = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
coul_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
; We are not transforming any bonded or restrained interactions
fep_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
sc-alpha = 0.5
sc-coul = no ; linear interpolation of Coulomb (none in this case)
sc-power = 1
sc-sigma = 0.3
nstdhdl = 10

I am using lamda0 state as solute-solvent fully decoupled (no vdw + no coulomb interactions) and lambda1 fully coupled (vdw and coulomb are fully switched on). In order to perturb the vdw and coulomb interactions only I am using an array of lamda values for the vdw_lambdas and coul_lambdas. I am setting fep_lambdas to an array of 0.00 as I don’t want any other lambda values or to perturb any other things like bonded interactions, mass, temperature, restraints, etc.

According to the mdp file details on gromacs website, it is mentioned that setting fep_lambdas to an array of values will use the same lambda values for vdw, coulomb, bond, temperature, restrains, mass lambda vectors. In the methanol solvation free energy tutorial on gromacs website, they used:

couple-lambda1 = vdwq
couple-lambda0 = none
fep-lambdas = 0.0 0.2 0.4 0.6 0.8 0.9 1.0 ;
Now it was mentioned that to avoid longer simulations, in this case the vdw and coulomb interactions are varied together (which might not be ideal always right?). Also if we set an array of lamda values to fep-lambdas, does it not mean that in addition to vdw and coulomb interactions, bonded, temperature, restrains—→ all these parameters are also going to get perturbed?

In my calculations I only want to vary vdw and coulomb interactions, setting fep_lambdas to null array would mean other *_lambdas are set to null arrays except for vdw_lambdas and coulomb_lambdas. Is this the standard way for calculating solvation free energy? Or it is mandatory to use same array of values for bonded_lambdas, temperature_lambdas, restraint_lambdas, mass_lambdas? Does it have anything to do with the:
couple-lambda1 = vdwq
couple-lambda0 = none

Is it that if we use couple_lambda1=vdwq, that means all other lambda values are irrelevant except for vdw and coulomb?

Also if we want to perturb the bonded interactions, mass, temperature, restraints, etc. what should we use for couple_lambda1?

Any help would be much appreciated, thank you

2

Hi! Your mdp parameters look alright to me for calculating SFEs.

  • Right. It’s usually better to switch of Coulomb interactions first, lest you get overlapping charges which can lead to your system exploding.
  • They would, if you had set alternative values (state A and B topology) for them.
  • Yes, that’s the standard. You don’t have to set all the other arrays explicitly. the couple-lambdaX indicates which non-bonded interactions you’re (de)coupling at the lambda end states. Your setup should be fine for what you want to do.
  • The couple-lambdaX are basically for convenience if you want to turn non-bondeds (vdw and/or q) completely on or off, so that you don’t have to alter the topology file only to add a bunch of zeros. For the other XXX_lambdas, you would need to set explicit A and B states, either in the .mdp or .top file. How you use these two options in conjunction depends entirely on what you want to achieve.

Hope this helps.

3

Thank you for detailed explanations. Is there any resources or documents that might point to the setting up of FEP simulations where one might want to perturb temperature, bonded interactions, mass, etc. The settings is not very clear from gromacs manual.

4

I’m not aware of any tutorials for more involved FEP simulations other than SFE. AFAIK, temperature-lambdas is only used for simulated tempering. The others are needed e.g. for absolute/relative binding free energy calculations. For these, you might want to have a look at pmx, for example.