# Implausible values in free energy calculations
**URL:** https://gromacs.bioexcel.eu/t/implausible-values-in-free-energy-calculations/7645
**Category:** User discussions
**Tags:** mdp-parameters, mdrun, free-energy
**Created:** [November 14, 2023, 2:29pm UTC](https://gromacs.bioexcel.eu/t/implausible-values-in-free-energy-calculations/7645 "2023-11-14T14:29:38Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@alkosh](https://gromacs.bioexcel.eu/u/alkosh)
#### Post date: [November 14, 2023, 2:29pm UTC](https://gromacs.bioexcel.eu/t/implausible-values-in-free-energy-calculations/7645/1 "2023-11-14T14:29:38Z")
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
---
### Author:  [@hess](https://gromacs.bioexcel.eu/u/hess)
#### Post date: [November 14, 2023, 6:22pm UTC](https://gromacs.bioexcel.eu/t/implausible-values-in-free-energy-calculations/7645/2 "2023-11-14T18:22:13Z")
How large is your box? Maybe this is due to issue 4665, which has been fixed in 2022.5.
---
### Author:  [@alkosh](https://gromacs.bioexcel.eu/u/alkosh)
#### Post date: [November 15, 2023, 6:14pm UTC](https://gromacs.bioexcel.eu/t/implausible-values-in-free-energy-calculations/7645/3 "2023-11-15T18:14:57Z")
I have appr. 15nmx15nmx10nm box. And the NP has about 13000 atoms (plus solvent). Today I tried to run the same calculations on 2022.5; unfortunately, it gives the same problem.
---
### Author:  [@hess](https://gromacs.bioexcel.eu/u/hess)
#### Post date: [November 15, 2023, 6:30pm UTC](https://gromacs.bioexcel.eu/t/implausible-values-in-free-energy-calculations/7645/4 "2023-11-15T18:30:28Z")
I have no clue then.
Have you tried turning off the bonded interactions along with the vdw interactions?
---
### Author:  [@cla254](https://gromacs.bioexcel.eu/u/cla254)
#### Post date: [February 16, 2024, 9:53pm UTC](https://gromacs.bioexcel.eu/t/implausible-values-in-free-energy-calculations/7645/5 "2024-02-16T21:53:43Z")
Hi Alkosh,
Were you able to fix the issue you were experiencing? I also have a similar issue.