GROMACS version: gromacs-2022.4
GROMACS modification: Yes/No
Hello! I am trying to calculate the free energy of a ligand on the nanoparticle. While calculating NP-Ligand - NP-Dummy free energy I faced with a strange behaviour of gromacs:
Detailed results in kT (see help for explanation):
lam_A lam_B DG +/- s_A +/- s_B +/- stdev +/-
0 1 32.23 0.19 1.39 0.09 1.43 0.08 1.77 0.04
1 2 27.24 0.20 3.56 0.39 3.82 0.12 4.08 0.50
2 3 22.33 0.15 1.09 0.17 0.84 0.10 1.36 0.07
3 4 20.90 0.08 0.59 0.12 0.56 0.12 1.11 0.04
4 5 19.51 0.08 0.82 0.18 0.85 0.18 1.20 0.12
5 6 -697.74 2.42 -107.13 2.33 -697.90 2.42 117213966.31 6690079.73
6 7 3.74 0.16 2.93 0.17 4.56 0.13 3.29 0.16
7 8 3.34 0.19 3.04 0.12 4.24 0.19 3.28 0.12
8 9 2.56 0.18 3.34 0.12 4.78 0.18 3.67 0.11
9 10 2.94 0.16 2.67 0.16 4.07 0.22 3.04 0.15
10 11 2.84 0.35 3.65 0.21 6.05 1.10 3.96 0.23
11 12 -0.64 0.36 5.81 0.66 7.07 1.11 5.82 1.47
12 13 -6.54 1.36 9.87 1.08 12.12 1.48 21.97 2.29
13 14 -8.86 0.79 7.16 0.51 8.18 0.86 9.32 1.27
14 15 -2.44 0.23 3.13 0.29 3.41 0.33 3.23 0.43
WARNING: Some of these results violate the Second Law of Thermodynamics:
This is can be the result of severe undersampling, or (more likely)
there is something wrong with the simulations.
The dG is incorrect on the step where coulomb interactions are switched off and the next lambda starts to switch off the vdw interactions. Do you know what could be the reason for that fluctuation? There are no any warnings during production but it appears in “gmx bar”.
Thank you in advance. Here I added md.mdp file that I use.
; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 500000 ; 1 ns
nstcomm = 100
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
nstxout-compressed = 0
cutoff-scheme = verlet
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.2
coulombtype = PME
rcoulomb = 1.2
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
DispCorr = EnerPres
fourierspacing = 0.12
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
tc-grps = NP_Mol ETH_Water ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 333 333 ; reference temperature, one for each group, in K
; Pressure coupling is on for NPT
Pcoupl = Parrinello-Rahman
tau_p = 1.0
compressibility = 4.5e-05
ref_p = 1.0
refcoord_scaling = com
; Free energy control stuff
free_energy = yes
init_lambda_state = 0
delta_lambda = 0
calc_lambda_neighbors = 1
couple-moltype = PDM
couple-lambda0 = vdw-q
couple-lambda1 = none
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
vdw_lambdas = 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
coul_lambdas = 0.00 0.20 0.40 0.60 0.80 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
bonded_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
restraint_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
sc-alpha = 0.5
sc-coul = no
sc-power = 1
sc-sigma = 0.3
nstdhdl = 10
gen_vel = no
constraints = h-bonds
constraint-algorithm = lincs
continuation = yes
lincs-order = 12