Free energy calc gives "WARNING: Some of these results violate the Second Law of Thermodynamics"

User discussions
, , ,
1

GROMACS version: 2020
GROMACS modification: No

Hello,

I am trying to run a free energy perturbation of a ligand bound to a protein using alchemical free energy calculations. However, after running bar on xvg files I get “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.” In order to increase sampling, I added more lambdas; currently I have 40 lambdas spaced out linearly from turning off electrostatics to Lennard-Jones. When I run solvation energy calculation of the ligand using same mdp file, I do not get the same warning. Below you can find the mdp file for the production run and below that a log file with the warning. Let me know if you need anything else!

Thank you!

mdp file:

;===========================================
; SD Run
;===========================================

; Run Control
;-------------------------------------------

integrator = sd
tinit = 0
dt = 0.002
nsteps = 5000000 ; 10ns
nstcomm = 100

; Output Control
;-------------------------------------------

nstlog = 500
nstenergy = 500
nstxout-compressed = 500
compressed-x-grps = Protein_BIT

; Neighbor Searching
;-------------------------------------------

cutoff-scheme = Verlet
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.2

; Electrostatics
;-------------------------------------------

coulombtype = PME
rcoulomb = 1.2

; van der Waals
;-------------------------------------------

vdwtype = Cut-off
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2

; Long Range Dispersion Correction
;-------------------------------------------

DispCorr = EnerPres

; Spacing
;-------------------------------------------

fourierspacing = 0.10

; EWALD/PME/PPPM parameters
;-------------------------------------------

pme_order = 4
ewald_rtol = 1e-06
epsilon_surface = 0

; Temperature Coupling
;-------------------------------------------

tc_grps = Protein_BIT Water_and_Ions ;tcoupl is implicitly handled by the sd integrator
tau_t = 2.0 2.0
ref_t = 298 298

; Pressure Coupling
;-------------------------------------------

Pcoupl = Parrinello-Rahman
tau_p = 1.0
compressibility = 4.5e-05
ref_p = 1.0

; Velocities
; ------------------------------------------

gen_vel = no

; Bonds
;-------------------------------------------

constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 12
lincs-iter = 2
; Continuation
;-------------------------------------------

continuation = yes

; Free Energy Parameters
; ------------------------------------------

free_energy = yes
init_lambda_state = 0
delta_lambda = 0
calc_lambda_neighbors = 1
couple-moltype = BIT
couple-lambda0 = vdw-q
couple-lambda1 = none
couple-intramol = no
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
vdw_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 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

; Soft-core Parameters
;-------------------------------------------

sc-alpha = 0.5
sc-coul = no
sc-power = 1
sc-sigma = 0.3
nstdhdl = 100

BAR log file:

Command line:
gmx bar -f md0.xvg md10.xvg md11.xvg md12.xvg md13.xvg md14.xvg md15.xvg md16.xvg md17.xvg md18.xvg md19.xvg md1.xvg md20.xvg md21.xvg md22.xvg md23.xvg md24.xvg md25.xvg md26.xvg md27.xvg md28.xvg md29.xvg md2.xvg md30.xvg md31.xvg md32.xvg md33.xvg md34.xvg md35.xvg md36.xvg md37.xvg md38.xvg md39.xvg md3.xvg md40.xvg md4.xvg md5.xvg md6.xvg md7.xvg md8.xvg md9.xvg -b 1000 -o bar.xvg -oi barint.xvg

md0.xvg: Ignoring set ‘pV (kJ/mol)’.
md0.xvg: 0.0 - 10000.0; lambda = (0, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0, 0) (50001 pts)
delta H to (0.05, 0) (50001 pts)

md10.xvg: Ignoring set ‘pV (kJ/mol)’.
md10.xvg: 0.0 - 10000.0; lambda = (0.5, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.45, 0) (50001 pts)
delta H to (0.5, 0) (50001 pts)
delta H to (0.55, 0) (50001 pts)

md11.xvg: Ignoring set ‘pV (kJ/mol)’.
md11.xvg: 0.0 - 10000.0; lambda = (0.55, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.5, 0) (50001 pts)
delta H to (0.55, 0) (50001 pts)
delta H to (0.6, 0) (50001 pts)

md12.xvg: Ignoring set ‘pV (kJ/mol)’.
md12.xvg: 0.0 - 10000.0; lambda = (0.6, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.55, 0) (50001 pts)
delta H to (0.6, 0) (50001 pts)
delta H to (0.65, 0) (50001 pts)

md13.xvg: Ignoring set ‘pV (kJ/mol)’.
md13.xvg: 0.0 - 10000.0; lambda = (0.65, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.6, 0) (50001 pts)
delta H to (0.65, 0) (50001 pts)
delta H to (0.7, 0) (50001 pts)

md14.xvg: Ignoring set ‘pV (kJ/mol)’.
md14.xvg: 0.0 - 10000.0; lambda = (0.7, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.65, 0) (50001 pts)
delta H to (0.7, 0) (50001 pts)
delta H to (0.75, 0) (50001 pts)

md15.xvg: Ignoring set ‘pV (kJ/mol)’.
md15.xvg: 0.0 - 10000.0; lambda = (0.75, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.7, 0) (50001 pts)
delta H to (0.75, 0) (50001 pts)
delta H to (0.8, 0) (50001 pts)

md16.xvg: Ignoring set ‘pV (kJ/mol)’.
md16.xvg: 0.0 - 10000.0; lambda = (0.8, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.75, 0) (50001 pts)
delta H to (0.8, 0) (50001 pts)
delta H to (0.85, 0) (50001 pts)

md17.xvg: Ignoring set ‘pV (kJ/mol)’.
md17.xvg: 0.0 - 10000.0; lambda = (0.85, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.8, 0) (50001 pts)
delta H to (0.85, 0) (50001 pts)
delta H to (0.9, 0) (50001 pts)

md18.xvg: Ignoring set ‘pV (kJ/mol)’.
md18.xvg: 0.0 - 10000.0; lambda = (0.9, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.85, 0) (50001 pts)
delta H to (0.9, 0) (50001 pts)
delta H to (0.95, 0) (50001 pts)

md19.xvg: Ignoring set ‘pV (kJ/mol)’.
md19.xvg: 0.0 - 10000.0; lambda = (0.95, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.9, 0) (50001 pts)
delta H to (0.95, 0) (50001 pts)
delta H to (1, 0) (50001 pts)

md1.xvg: Ignoring set ‘pV (kJ/mol)’.
md1.xvg: 0.0 - 10000.0; lambda = (0.05, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0, 0) (50001 pts)
delta H to (0.05, 0) (50001 pts)
delta H to (0.1, 0) (50001 pts)

md20.xvg: Ignoring set ‘pV (kJ/mol)’.
md20.xvg: 0.0 - 10000.0; lambda = (1, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.95, 0) (50001 pts)
delta H to (1, 0) (50001 pts)
delta H to (1, 0.05) (50001 pts)

md21.xvg: Ignoring set ‘pV (kJ/mol)’.
md21.xvg: 0.0 - 10000.0; lambda = (1, 0.05)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0) (50001 pts)
delta H to (1, 0.05) (50001 pts)
delta H to (1, 0.1) (50001 pts)

md22.xvg: Ignoring set ‘pV (kJ/mol)’.
md22.xvg: 0.0 - 10000.0; lambda = (1, 0.1)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.05) (50001 pts)
delta H to (1, 0.1) (50001 pts)
delta H to (1, 0.15) (50001 pts)

md23.xvg: Ignoring set ‘pV (kJ/mol)’.
md23.xvg: 0.0 - 10000.0; lambda = (1, 0.15)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.1) (50001 pts)
delta H to (1, 0.15) (50001 pts)
delta H to (1, 0.2) (50001 pts)

md24.xvg: Ignoring set ‘pV (kJ/mol)’.
md24.xvg: 0.0 - 10000.0; lambda = (1, 0.2)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.15) (50001 pts)
delta H to (1, 0.2) (50001 pts)
delta H to (1, 0.25) (50001 pts)

md25.xvg: Ignoring set ‘pV (kJ/mol)’.
md25.xvg: 0.0 - 10000.0; lambda = (1, 0.25)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.2) (50001 pts)
delta H to (1, 0.25) (50001 pts)
delta H to (1, 0.3) (50001 pts)

md26.xvg: Ignoring set ‘pV (kJ/mol)’.
md26.xvg: 0.0 - 10000.0; lambda = (1, 0.3)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.25) (50001 pts)
delta H to (1, 0.3) (50001 pts)
delta H to (1, 0.35) (50001 pts)

md27.xvg: Ignoring set ‘pV (kJ/mol)’.
md27.xvg: 0.0 - 10000.0; lambda = (1, 0.35)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.3) (50001 pts)
delta H to (1, 0.35) (50001 pts)
delta H to (1, 0.4) (50001 pts)

md28.xvg: Ignoring set ‘pV (kJ/mol)’.
md28.xvg: 0.0 - 10000.0; lambda = (1, 0.4)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.35) (50001 pts)
delta H to (1, 0.4) (50001 pts)
delta H to (1, 0.45) (50001 pts)

md29.xvg: Ignoring set ‘pV (kJ/mol)’.
md29.xvg: 0.0 - 10000.0; lambda = (1, 0.45)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.4) (50001 pts)
delta H to (1, 0.45) (50001 pts)
delta H to (1, 0.5) (50001 pts)

md2.xvg: Ignoring set ‘pV (kJ/mol)’.
md2.xvg: 0.0 - 10000.0; lambda = (0.1, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.05, 0) (50001 pts)
delta H to (0.1, 0) (50001 pts)
delta H to (0.15, 0) (50001 pts)

md30.xvg: Ignoring set ‘pV (kJ/mol)’.
md30.xvg: 0.0 - 10000.0; lambda = (1, 0.5)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.45) (50001 pts)
delta H to (1, 0.5) (50001 pts)
delta H to (1, 0.55) (50001 pts)

md31.xvg: Ignoring set ‘pV (kJ/mol)’.
md31.xvg: 0.0 - 10000.0; lambda = (1, 0.55)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.5) (50001 pts)
delta H to (1, 0.55) (50001 pts)
delta H to (1, 0.6) (50001 pts)

md32.xvg: Ignoring set ‘pV (kJ/mol)’.
md32.xvg: 0.0 - 10000.0; lambda = (1, 0.6)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.55) (50001 pts)
delta H to (1, 0.6) (50001 pts)
delta H to (1, 0.65) (50001 pts)

md33.xvg: Ignoring set ‘pV (kJ/mol)’.
md33.xvg: 0.0 - 10000.0; lambda = (1, 0.65)
dH/dl & foreign lambdas:

Back Off! I just backed up bar.xvg to ./#bar.xvg.1#

Back Off! I just backed up barint.xvg to ./#barint.xvg.1#
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.6) (50001 pts)
delta H to (1, 0.65) (50001 pts)
delta H to (1, 0.7) (50001 pts)

md34.xvg: Ignoring set ‘pV (kJ/mol)’.
md34.xvg: 0.0 - 10000.0; lambda = (1, 0.7)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.65) (50001 pts)
delta H to (1, 0.7) (50001 pts)
delta H to (1, 0.75) (50001 pts)

md35.xvg: Ignoring set ‘pV (kJ/mol)’.
md35.xvg: 0.0 - 10000.0; lambda = (1, 0.75)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.7) (50001 pts)
delta H to (1, 0.75) (50001 pts)
delta H to (1, 0.8) (50001 pts)

md36.xvg: Ignoring set ‘pV (kJ/mol)’.
md36.xvg: 0.0 - 10000.0; lambda = (1, 0.8)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.75) (50001 pts)
delta H to (1, 0.8) (50001 pts)
delta H to (1, 0.85) (50001 pts)

md37.xvg: Ignoring set ‘pV (kJ/mol)’.
md37.xvg: 0.0 - 10000.0; lambda = (1, 0.85)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.8) (50001 pts)
delta H to (1, 0.85) (50001 pts)
delta H to (1, 0.9) (50001 pts)

md38.xvg: Ignoring set ‘pV (kJ/mol)’.
md38.xvg: 0.0 - 10000.0; lambda = (1, 0.9)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.85) (50001 pts)
delta H to (1, 0.9) (50001 pts)
delta H to (1, 0.95) (50001 pts)

md39.xvg: Ignoring set ‘pV (kJ/mol)’.
md39.xvg: 0.0 - 10000.0; lambda = (1, 0.95)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.9) (50001 pts)
delta H to (1, 0.95) (50001 pts)
delta H to (1, 1) (50001 pts)

md3.xvg: Ignoring set ‘pV (kJ/mol)’.
md3.xvg: 0.0 - 10000.0; lambda = (0.15, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.1, 0) (50001 pts)
delta H to (0.15, 0) (50001 pts)
delta H to (0.2, 0) (50001 pts)

md40.xvg: Ignoring set ‘pV (kJ/mol)’.
md40.xvg: 0.0 - 10000.0; lambda = (1, 1)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (1, 0.95) (50001 pts)
delta H to (1, 1) (50001 pts)

md4.xvg: Ignoring set ‘pV (kJ/mol)’.
md4.xvg: 0.0 - 10000.0; lambda = (0.2, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.15, 0) (50001 pts)
delta H to (0.2, 0) (50001 pts)
delta H to (0.25, 0) (50001 pts)

md5.xvg: Ignoring set ‘pV (kJ/mol)’.
md5.xvg: 0.0 - 10000.0; lambda = (0.25, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.2, 0) (50001 pts)
delta H to (0.25, 0) (50001 pts)
delta H to (0.3, 0) (50001 pts)

md6.xvg: Ignoring set ‘pV (kJ/mol)’.
md6.xvg: 0.0 - 10000.0; lambda = (0.3, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.25, 0) (50001 pts)
delta H to (0.3, 0) (50001 pts)
delta H to (0.35, 0) (50001 pts)

md7.xvg: Ignoring set ‘pV (kJ/mol)’.
md7.xvg: 0.0 - 10000.0; lambda = (0.35, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.3, 0) (50001 pts)
delta H to (0.35, 0) (50001 pts)
delta H to (0.4, 0) (50001 pts)

md8.xvg: Ignoring set ‘pV (kJ/mol)’.
md8.xvg: 0.0 - 10000.0; lambda = (0.4, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.35, 0) (50001 pts)
delta H to (0.4, 0) (50001 pts)
delta H to (0.45, 0) (50001 pts)

md9.xvg: Ignoring set ‘pV (kJ/mol)’.
md9.xvg: 0.0 - 10000.0; lambda = (0.45, 0)
dH/dl & foreign lambdas:
dH/dl (coul-lambda) (50001 pts)
dH/dl (vdw-lambda) (50001 pts)
delta H to (0.4, 0) (50001 pts)
delta H to (0.45, 0) (50001 pts)
delta H to (0.5, 0) (50001 pts)

Samples in time interval: 0.000 - 10000.000
Removing samples outside of: 1000.000 - 10000.000

Temperature: 298 K

Detailed results in kT (see help for explanation):

lam_A lam_B DG +/- s_A +/- s_B +/- stdev +/-
0 1 2.19 0.02 0.14 0.03 0.14 0.03 0.30 0.01
1 2 1.86 0.07 0.19 0.05 0.20 0.05 0.36 0.02
2 3 1.83 0.04 -0.17 0.08 -0.17 0.08 0.35 0.02
3 4 1.84 0.02 0.16 0.05 0.16 0.05 0.31 0.02
4 5 1.81 0.03 -0.13 0.02 -0.13 0.02 0.30 0.01
5 6 1.73 0.02 0.21 0.03 0.21 0.03 0.35 0.02
6 7 1.35 0.02 0.18 0.03 0.17 0.03 0.32 0.01
7 8 1.23 0.03 -0.06 0.01 -0.05 0.01 0.26 0.01
8 9 1.22 0.04 0.07 0.01 0.07 0.01 0.28 0.00
9 10 1.12 0.02 0.02 0.03 0.02 0.03 0.29 0.01
10 11 1.05 0.02 0.05 0.04 0.05 0.04 0.33 0.01
11 12 1.05 0.05 -0.05 0.03 -0.04 0.03 0.38 0.01
12 13 0.90 0.04 0.19 0.03 0.19 0.03 0.42 0.01
13 14 0.66 0.01 0.05 0.02 0.05 0.02 0.35 0.00
14 15 0.57 0.00 0.04 0.01 0.03 0.01 0.34 0.00
15 16 0.46 0.01 0.08 0.01 0.07 0.01 0.32 0.01
16 17 0.44 0.04 -0.05 0.02 -0.06 0.02 0.30 0.01
17 18 0.37 0.02 0.12 0.04 0.12 0.04 0.29 0.02
18 19 0.21 0.01 0.04 0.02 0.04 0.02 0.21 0.00
19 20 0.12 0.02 0.05 0.02 0.05 0.02 0.21 0.01
20 21 2.04 0.03 0.38 0.03 0.48 0.04 0.88 0.02
21 22 2.05 0.03 0.32 0.02 0.41 0.02 0.88 0.01
22 23 2.08 0.02 0.36 0.02 0.45 0.02 0.88 0.01
23 24 2.06 0.02 0.38 0.03 0.48 0.04 0.89 0.01
24 25 2.04 0.07 0.37 0.07 0.48 0.09 0.92 0.05
25 26 2.25 0.06 0.12 0.05 0.17 0.05 0.80 0.01
26 27 2.34 0.02 0.46 0.01 0.54 0.01 0.88 0.01
27 28 2.15 0.15 0.47 0.14 0.64 0.17 1.00 0.10
28 29 1.42 0.17 0.96 0.10 1.14 0.15 1.36 0.08
29 30 1.61 0.08 -0.10 0.08 -0.04 0.08 0.92 0.02
30 31 2.16 0.18 0.35 0.16 0.50 0.18 0.97 0.09
31 32 1.81 0.19 0.70 0.13 0.82 0.15 1.18 0.06
32 33 1.46 0.22 0.56 0.28 0.81 0.32 1.26 0.14
33 34 1.01 0.30 0.80 0.15 0.91 0.17 1.37 0.06
34 35 0.15 0.12 1.26 0.09 1.55 0.11 1.62 0.08
35 36 -1.54 0.14 1.90 0.11 2.10 0.13 2.20 0.11
36 37 -1.96 0.10 0.62 0.11 0.63 0.11 1.17 0.04
37 38 -1.17 0.10 0.70 0.06 0.76 0.06 1.09 0.04
38 39 -0.21 0.05 0.31 0.05 0.39 0.05 0.92 0.02
39 40 0.74 0.03 0.43 0.03 0.55 0.03 1.01 0.01

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.

Final results in kJ/mol:

point 0 - 1, DG 5.44 +/- 0.05
point 1 - 2, DG 4.60 +/- 0.17
point 2 - 3, DG 4.54 +/- 0.10
point 3 - 4, DG 4.57 +/- 0.05
point 4 - 5, DG 4.48 +/- 0.07
point 5 - 6, DG 4.29 +/- 0.05
point 6 - 7, DG 3.35 +/- 0.05
point 7 - 8, DG 3.05 +/- 0.07
point 8 - 9, DG 3.01 +/- 0.10
GROMACS reminds you: “There is no such thing as free energy. Anyone who advocates it does not know what he is talking about.” (Alireza Haghighat)

point 9 - 10, DG 2.78 +/- 0.06
point 10 - 11, DG 2.61 +/- 0.05
point 11 - 12, DG 2.61 +/- 0.13
point 12 - 13, DG 2.23 +/- 0.10
point 13 - 14, DG 1.63 +/- 0.02
point 14 - 15, DG 1.42 +/- 0.01
point 15 - 16, DG 1.14 +/- 0.03
point 16 - 17, DG 1.09 +/- 0.09
point 17 - 18, DG 0.92 +/- 0.04
point 18 - 19, DG 0.52 +/- 0.03
point 19 - 20, DG 0.31 +/- 0.05
point 20 - 21, DG 5.05 +/- 0.09
point 21 - 22, DG 5.07 +/- 0.08
point 22 - 23, DG 5.16 +/- 0.05
point 23 - 24, DG 5.10 +/- 0.05
point 24 - 25, DG 5.06 +/- 0.17
point 25 - 26, DG 5.58 +/- 0.14
point 26 - 27, DG 5.80 +/- 0.04
point 27 - 28, DG 5.33 +/- 0.38
point 28 - 29, DG 3.51 +/- 0.42
point 29 - 30, DG 3.98 +/- 0.20
point 30 - 31, DG 5.35 +/- 0.45
point 31 - 32, DG 4.50 +/- 0.46
point 32 - 33, DG 3.61 +/- 0.54
point 33 - 34, DG 2.49 +/- 0.73
point 34 - 35, DG 0.38 +/- 0.30
point 35 - 36, DG -3.83 +/- 0.34
point 36 - 37, DG -4.87 +/- 0.26
point 37 - 38, DG -2.91 +/- 0.25
point 38 - 39, DG -0.52 +/- 0.13
point 39 - 40, DG 1.84 +/- 0.07

total 0 - 40, DG 110.23 +/- 2.42

2

Hi
Sorry I do not understand which is your question.
In general different mutations may require different lambda points to get the same accuracy in the results. It can be that a ligand transformation in a protein or in water solution require different sampling to converge.
\Alessandra

3

Hi,

Thank you for the reply! My question is if there are any suggestions on how to fix the warning. After encountering the warning I have added more lambdas (~50 windows total; 10ns each) which hasn’t changed it. I attached the mdp file as the warning states that there might be something wrong with the simulation, so I wonder if there is anything odd about my mdp file. I tried it on different protein isoforms with the same ligand and still received the warning, so I am not sure what to try next. I have seen people run fep with less windows and shorter simulation per window, which makes me think this is not an issue of undersampling.

Let me know what you think!

4

Hi,
As far as I understood, you decouple a ligand bound to a protein. This results with a weakly coupled ligand wanders through our system. This is a poorly reversible situation: there are suddenly very few states that map from a weakly coupled to a more strongly coupled molecule, which will drastically reduce the accuracy of the free energy calculation.

A solution is forcing the ligand to stay at a specific position relative to the protein, for example forcing the ligand to stay at its native position even when it has been fully de-coupled. Once the free energy has been calculated, care must be taken to correct for the fact that the ligand was forcing to stay in its native position. This can easily be done analytically.

An alternative is to calculate relative binding free energy. If you want to compare the binding affinity for different ligands to same protein, you can use the thermodynamic cycle in Fig 10 (Free energy calculations — GROMACS 2021.3 documentation ). Here ligand A is mutated in ligand B in solution and in the protein. Then the relative binding free energy can be calculated.

Best regards
Alessandra

1 Like
5

Hi,

Thank you for the reply! Yes, that is correct. I decouple the ligand bound to the protein, and in order to not let it wonder around I added Boresch-style restraint to the ligand, which I account for through analytical correction. The mdp file itself doesn’t have lambda-restraints, however the restraint is applied as the ligand stays bound to the pocket, even when fully decoupled!

Thank you

6

I also get this warning for just my ligand in water system as well as my protein-ligand-water system. Why would I get this warning for just my ligand in water system. I am interested in theoretically first and then I can attach all my ligand solvation mdp files or at least one and describe nvt/npt/md differences I made for vdw-lambdas and coul-lambdas with respect to turning off vdw or electrostatics for state A → state B.