Dear gromacs users,
I am calculating the solvation free energy of an organic molecule (merocyanine) in water and in ethylene glycol. I prepared two simulation systems by solvating one MC molecule in a waterbox (4x4x4 nm) and in a box containing ethylene glycol as solvent. The energy minimized and equilibrated the system for a few ns and used that equilibrated configuration for alchemical free energy calculation.
For alchemical free energy calculation I have used energy minimization, NVT equilibration and velocity generation, NPT equilibration and then production run, at each lambda values. Following are my mdp file setting for free-energy calculation: (I am using gromacs 2024, sd integrator and C-rescale barostat)
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
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
; Options for the decoupling
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 have kept couple-intramol = yes as setting it to “no” was generating many errors related to pair interactions and setting this to yes is not throwing any errors, I ran the NVT and NPT equilibrations (with alchemical setup) for 500 ps each and the production MD for 1 ns. After that I copied all md*.xvg files in a separate directory and ran: gmx bar -f md*xvg -o -oi and the output was something like this:
lam_A lam_B DG +/- s_A +/- s_B +/- stdev +/-
0 1 1.92 0.10 1.42 0.11 1.67 0.12 1.84 0.07
1 2 3.64 0.22 5.93 0.22 4.93 0.28 4.98 0.31
2 3 -0.64 0.24 6.43 0.22 3.43 0.20 4.04 0.25
3 4 -2.86 0.09 3.11 0.07 1.85 0.08 2.42 0.05
4 5 -3.43 0.07 2.22 0.06 1.42 0.05 2.00 0.03
5 6 -3.57 0.04 1.82 0.09 1.19 0.05 1.80 0.04
6 7 -3.41 0.03 1.80 0.05 1.18 0.03 1.75 0.02
7 8 -3.06 0.02 1.65 0.02 1.09 0.01 1.69 0.01
8 9 -2.49 0.02 1.85 0.02 1.19 0.01 1.77 0.01
9 10 -1.84 0.02 2.07 0.01 1.30 0.01 1.86 0.01
10 11 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
11 12 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
12 13 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
13 14 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
14 15 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
15 16 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
16 17 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
17 18 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
18 19 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
19 20 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
Final results in kJ/mol:
point 0 - 1, DG 4.77 +/- 0.25
point 1 - 2, DG 9.01 +/- 0.53
point 2 - 3, DG -1.59 +/- 0.58
point 3 - 4, DG -7.10 +/- 0.22
point 4 - 5, DG -8.51 +/- 0.16
point 5 - 6, DG -8.84 +/- 0.11
point 6 - 7, DG -8.45 +/- 0.07
point 7 - 8, DG -7.58 +/- 0.04
point 8 - 9, DG -6.17 +/- 0.04
point 9 - 10, DG -4.57 +/- 0.06
point 10 - 11, DG 0.00 +/- 0.00
point 11 - 12, DG 0.00 +/- 0.00
point 12 - 13, DG 0.00 +/- 0.00
point 13 - 14, DG 0.00 +/- 0.00
point 14 - 15, DG 0.00 +/- 0.00
point 15 - 16, DG 0.00 +/- 0.00
point 16 - 17, DG 0.00 +/- 0.00
point 17 - 18, DG 0.00 +/- 0.00
point 18 - 19, DG 0.00 +/- 0.00
point 19 - 20, DG 0.00 +/- 0.00
total 0 - 20, DG -39.05 +/- 1.28
Now what I see that once the coulomb interactions are switched on (11th lambda point onwards i.e. lambda vector 10 onwards) the DG values are all zero and. Is there any error(s) in the mdp file setup? Is there any correlation with the output and the mdp setting that I am using?
Any help would be much appreciated, thank you.