Value too high for SFE (charmm-gui files)

User discussions
,
1

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 ?

2

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.

3

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 ?

4

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.

5

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

6

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 !

1 Like