# Value too high for SFE (charmm-gui files)
**URL:** https://gromacs.bioexcel.eu/t/value-too-high-for-sfe-charmm-gui-files/12463
**Category:** User discussions
**Tags:** mdp-parameters, simulation-setup
**Created:** [July 30, 2025, 11:17am UTC](https://gromacs.bioexcel.eu/t/value-too-high-for-sfe-charmm-gui-files/12463 "2025-07-30T11:17:54Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@Usual\_Protein\_User](https://gromacs.bioexcel.eu/u/Usual_Protein_User)
#### Post date: [July 30, 2025, 11:17am UTC](https://gromacs.bioexcel.eu/t/value-too-high-for-sfe-charmm-gui-files/12463/1 "2025-07-30T11:17:54Z")
GROMACS version: 2025.2
GROMACS modification: No
Here post your question
i did a calculation for the Desolvation Energy of Tryptophan and got quite helpful tips for my Lambda-Settings and the result (Thanks @Prof.Hess and Imullender). However still somethings seems to be wrong in my setup:
I used the charmm-gui setup with the charmm forcefield
My result says:
Final results in kJ/mol:
point 0 - 1, DG 110.08 +/- 0.27
point 1 - 2, DG 77.99 +/- 0.60
point 2 - 3, DG 48.99 +/- 0.29
point 3 - 4, DG 15.16 +/- 0.10
point 4 - 5, DG 10.09 +/- 0.04
point 5 - 6, DG 5.75 +/- 0.05
point 6 - 7, DG 1.90 +/- 0.02
point 7 - 8, DG 6.14 +/- 0.02
point 8 - 9, DG 5.80 +/- 0.06
point 9 - 10, DG 5.44 +/- 0.07
point 10 - 11, DG 4.57 +/- 0.09
point 11 - 12, DG 3.47 +/- 0.12
point 12 - 13, DG 1.70 +/- 0.15
point 13 - 14, DG -0.82 +/- 0.56
point 14 - 15, DG -3.49 +/- 0.52
point 15 - 16, DG -6.13 +/- 0.42
point 16 - 17, DG -7.18 +/- 0.23
point 17 - 18, DG -5.90 +/- 0.08
point 18 - 19, DG -3.50 +/- 0.05
point 19 - 20, DG -1.08 +/- 0.06
total 0 - 20, DG 268.98 +/- 1.16
## Different Papers say the value should be arround 20/25 KJ/mol
My min\_mdp is:
; Run control
integrator = steep
nsteps = 5000
; EM criteria and other stuff
emtol = 100
emstep = 0.01
niter = 20
nbfgscorr = 10
; Output control
nstlog = 1
nstenergy = 1
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 1
ns\_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme\_order = 6
ewald\_rtol = 1e-06
epsilon\_surface = 0
; Temperature and pressure coupling are off during EM
tcoupl = no
pcoupl = no
; Free energy control stuff
free\_energy = yes
init\_lambda\_state = 0
delta\_lambda = 0
calc\_lambda\_neighbors = 1 ; only immediate neighboring windows
couple-moltype = PROA ; name of moleculetype to decouple
couple-lambda0 = vdw-q ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
; 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.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.70 0.75 0.80 0.85 0.90 0.95 1.0
coul\_lambdas = 0.0 0.2 0.4 0.6 0.7 0.8 0.9 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
; We are not transforming any bonded or restrained interactions
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 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 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1)
mass\_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
; Not doing simulated temperting here
temperature\_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
; No velocities during EM
gen\_vel = no
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Do not constrain the starting configuration
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
My equilibration\_nvt is:
; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 50000 ; 100 ps
nstcomm = 100
; Output control
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
nstxout-compressed = 0
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 20
ns\_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme\_order = 6
ewald\_rtol = 1e-06
epsilon\_surface = 0
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc\_grps = system
tau\_t = 1.0
ref\_t = 298
; Pressure coupling is off for NVT
Pcoupl = No
tau\_p = 0.5
compressibility = 4.5e-05
ref\_p = 1.0
; Free energy control stuff
free\_energy = yes
init\_lambda\_state = 0
delta\_lambda = 0
calc\_lambda\_neighbors = 1 ; only immediate neighboring windows
couple-moltype = PROA ; name of moleculetype to decouple
couple-lambda0 = vdw-q ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
; 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.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.70 0.75 0.80 0.85 0.90 0.95 1.0
coul\_lambdas = 0.0 0.2 0.4 0.6 0.7 0.8 0.9 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
; We are not transforming any bonded or restrained interactions
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 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 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1)
mass\_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
; Not doing simulated temperting here
temperature\_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
; Generate velocities to start
gen\_vel = yes
gen\_temp = 300
gen\_seed = -1
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Do not constrain the starting configuration
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
My equilibration\_npt is:
; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 50000 ; 100 ps
nstcomm = 100
; Output control
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
nstxout-compressed = 0
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 20
ns\_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme\_order = 6
ewald\_rtol = 1e-06
epsilon\_surface = 0
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc\_grps = system
tau\_t = 1.0
ref\_t = 298
; Pressure coupling is on for NPT
Pcoupl = Parrinello-Rahman
tau\_p = 1.0
compressibility = 4.5e-05
ref\_p = 1.0
; Free energy control stuff
free\_energy = yes
init\_lambda\_state = 0
delta\_lambda = 0
calc\_lambda\_neighbors = 1 ; only immediate neighboring windows
couple-moltype = PROA ; name of moleculetype to decouple
couple-lambda0 = vdw-q ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
; 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.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.70 0.75 0.80 0.85 0.90 0.95 1.0
coul\_lambdas = 0.0 0.2 0.4 0.6 0.7 0.8 0.9 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
; We are not transforming any bonded or restrained interactions
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 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 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1)
mass\_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
; Not doing simulated temperting here
temperature\_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
; Do not generate velocities
gen\_vel = no
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Constrain the starting configuration
; since we are continuing from NVT
continuation = yes
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
My production\_mdp is:
; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 50000 ; 100 ps
nstcomm = 100
; Output control
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
nstxout-compressed = 0
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 20
ns\_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme\_order = 6
ewald\_rtol = 1e-06
epsilon\_surface = 0
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc\_grps = system
tau\_t = 1.0
ref\_t = 298
; Pressure coupling is on for NPT
Pcoupl = Parrinello-Rahman
tau\_p = 1.0
compressibility = 4.5e-05
ref\_p = 1.0
; Free energy control stuff
free\_energy = yes
init\_lambda\_state = 0
delta\_lambda = 0
calc\_lambda\_neighbors = 1 ; only immediate neighboring windows
couple-moltype = PROA ; name of moleculetype to decouple
couple-lambda0 = vdw-q ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
; 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.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.70 0.75 0.80 0.85 0.90 0.95 1.0
coul\_lambdas = 0.0 0.2 0.4 0.6 0.7 0.8 0.9 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
; We are not transforming any bonded or restrained interactions
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 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 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1)
mass\_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
; Not doing simulated temperting here
temperature\_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
; Do not generate velocities
gen\_vel = no
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Constrain the starting configuration
; since we are continuing from NVT
continuation = yes
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
Has anyone an idea what might be missing here ?
---
### Author:  [@lmullender](https://gromacs.bioexcel.eu/u/lmullender)
#### Post date: [July 31, 2025, 1:48pm UTC](https://gromacs.bioexcel.eu/t/value-too-high-for-sfe-charmm-gui-files/12463/2 "2025-07-31T13:48:28Z")
Hi!
From what I can see, your mdp parameters look sensible to me, but your result is indeed far off from what it should be. The first thing that I would try is to extend your production simulations, because 100 ps is a bit short for a more complicated molecule like tryptophan. I would go at least 1 or 2 ns per lambda, to be safe.
If that doesn’t help, the next possible cause is the spacing of the lambda windows. A good thing to look out for is the error estimate for each lambda window, given in the `gmx bar` output. Where that estimate is very high, you might consider adding an additional intermediate lambda value. As to what “too high” means, I’m not super experienced myself, but anything above 0.2 - 0.3 or so would make me suspicious. Also, very high values of the free energy difference itself like you have for your first couple of lambda windows are kind of suspicious.
---
### Author:  [@Usual\_Protein\_User](https://gromacs.bioexcel.eu/u/Usual_Protein_User)
#### Post date: [August 4, 2025, 12:52pm UTC](https://gromacs.bioexcel.eu/t/value-too-high-for-sfe-charmm-gui-files/12463/3 "2025-08-04T12:52:25Z")
Hi !
Thanks again for your answer and your time to help ! I replied just now, after completing my MD-run with the new settings. I increased my running time up to 1 ns to get a feeling for the new setup and it took over 20 hours to finish (i believe i didn’t install my parallelization correctly, but that’s something for another topic).
However these are my new results, they haven’t really changed - so i would suggest, that running longer doesn’t do more for my result ?
Final results in kJ/mol:
point 0 - 1, DG 110.20 +/- 0.25
point 1 - 2, DG 77.27 +/- 0.17
point 2 - 3, DG 48.45 +/- 0.20
point 3 - 4, DG 15.29 +/- 0.05
point 4 - 5, DG 10.07 +/- 0.02
point 5 - 6, DG 5.78 +/- 0.03
point 6 - 7, DG 2.06 +/- 0.03
point 7 - 8, DG 6.06 +/- 0.03
point 8 - 9, DG 5.73 +/- 0.04
point 9 - 10, DG 5.29 +/- 0.02
point 10 - 11, DG 4.65 +/- 0.03
point 11 - 12, DG 3.73 +/- 0.04
point 12 - 13, DG 2.08 +/- 0.06
point 13 - 14, DG -1.28 +/- 0.10
point 14 - 15, DG -2.97 +/- 0.08
point 15 - 16, DG -5.95 +/- 0.06
point 16 - 17, DG -7.51 +/- 0.08
point 17 - 18, DG -6.03 +/- 0.03
point 18 - 19, DG -3.60 +/- 0.01
point 19 - 20, DG -1.14 +/- 0.02
total 0 - 20, DG 268.19 +/- 0.68
Do you think spacing now might work here ? It feels like the error is somewhere else, but i am not sure about that. Or is there a better way to calculate the SFE, maybe the regular way via the lambda states is too much effort for this usecase - i just read something about gmx\_MMPBSA. Do you know something about that ?
---
### Author:  [@lmullender](https://gromacs.bioexcel.eu/u/lmullender)
#### Post date: [August 4, 2025, 1:43pm UTC](https://gromacs.bioexcel.eu/t/value-too-high-for-sfe-charmm-gui-files/12463/4 "2025-08-04T13:43:17Z")
Hi! Ok, in that case it doesn’t seem like it was the simulation time. Like I said, my next hunch would be to look at the lambda spacings. In particular, I would add another one at vdw\_l=0.65, and at coul\_l=0.1,0.3,0.5 . Probably you can lower the simulation time again then to save on computational effort.
I haven’t used MM/PBSA so I can’t comment on that one. From what I understand, it might be a bit cheaper, but also less accurate than BAR.
---
### Author:  [@Usual\_Protein\_User](https://gromacs.bioexcel.eu/u/Usual_Protein_User)
#### Post date: [August 4, 2025, 2:38pm UTC](https://gromacs.bioexcel.eu/t/value-too-high-for-sfe-charmm-gui-files/12463/5 "2025-08-04T14:38:59Z")
Right, i will try it now with the new spacing method :) Was also thinking about reducing the simulation time as well. Okay thanks, read similar things about the accuracy.
I would like to ask another question: Due to the fact that i am new to gromacs, i don’t understand the whole Calculation process with the lambdas - so in theory, when i add another lambda settings (intramol=yes) it should need more computational power right ? But i only experienced a longer simulation time by increasing my **nsteps** in my prod.mdp. So when i increase my steps in my equilibration and minimization files it should take additionally more time, right ? Are the settings with the steps the only way to set my computational effort ? When my steps with this little molecule don’t play a role apparently, does it mean my steps are already fine enough to calculate a protein later ?
I hope i made my confusion clear
---
### Author:  [@Usual\_Protein\_User](https://gromacs.bioexcel.eu/u/Usual_Protein_User)
#### Post date: [August 6, 2025, 9:02pm UTC](https://gromacs.bioexcel.eu/t/value-too-high-for-sfe-charmm-gui-files/12463/6 "2025-08-06T21:02:12Z")
I got help from a research group in this field and finally understood my mistake. I rushed the papers and the detailed description of the molecules. In the papers i was looking at the results for the Tryptophan **Sidechain** not the zwitterionic form of tryptophan .. So the results between these two differ a lot. But it seems that my results for the zwitterionic did quite well (arround 260 KJ/mol in the paper). So thanks again for your time and effort and sorry for the confusion !