Error: Positive value for solvation free energy

User discussions
1

GROMACS version: 2025.2
GROMACS modification: No

I did a free energy calculation for tryptophan. For testing i just viewed my lambda16-18 run.
Somehow i got positive values at the End, shouldnt it be otherwise ?

Temperature: 298 K

Detailed results in kT (see help for explanation):

lam_A lam_B DG +/- s_A +/- s_B +/- stdev +/-
16 17 2.62 0.09 0.47 0.12 0.45 0.12 0.84 0.05
17 18 3.69 0.07 0.04 0.05 0.03 0.06 0.61 0.02

Final results in kJ/mol:

point 16 - 17, DG 6.49 +/- 0.22
point 17 - 18, DG 9.15 +/- 0.16

total 16 - 18, DG 15.65 +/- 0.20

These are my mdp. settings (for example my minimization.mdp)

; 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.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
; 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 = yes ; 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

2

Hi! Someone more experienced may correct me, but a couple of points:

  1. To get a meaningful result, you have to perform the BAR analysis over all lambda windows, so you should include them before you start reasoning about your results.
  2. Perhaps the main point: Why do you expect a negative result? We are just talking about free energy differences here, which can very well be positive. In the case of SFE calculations, a negative result indicates that the molecule prefers being in solution, while a positive result means that it prefers the gas phase.
  3. The way you have set it up, lambda=0 is the fully coupled state, while lambda=1 is the decoupled one.This will also affect the sign of the result.
  4. Instead of adding a soft-core potential for the Coulomb potential (sc-coul), you might want to consider first turning off Coulomb interactions before decoupling the Lennard Jones terms, to avoid problems with charges getting too close to each other. You can do this by first letting the coul-lambda go from 0 to 1 while keeping vdw-lambda at 0, and then ramping up the vdw-lambda while keeping coul-lambda at 1.
  5. Nit: To reduce clutter, there is no need to explicitly state the lambdas for bondeds, restraints etc. with all zeros, as that’s the default anyway.

Hope this helps!

3

Your solvation free energy is negative.

4

Thanks a lot for the reply !

I am really new to the MD-World. I did my first simulation in gromacs (SFE) with tryptophan by reproducing the workflow of SFE for methanol by the Tutorial from J.Lemkul. So i tried to compare my last lambda results here from the ones i got from the tutorial run (where it was negative at the end - including the total SFE).

From a paper where they also did SFE for aminoacids i found a negative total value of something arround -20;-25 KJ/mol for tryptophan, so iam trying to reproduce this right now.

Question for 4.) Definitely will try this. Do you think my steps are good as they are now ?

And Question for 5.) So you mean i can basically delete this rows then in this case ?

Thanks again - glad i found this great platform. For a beginner like me this gold here and LLMs get to their limit at some point. I will report how it went !

5

To be a bit less cryptic: your setup computes the DEsolvation energy, as the lambda=1 state is the decoupled state. To get the solvation free energy you need to multiply your result with -1.

1 Like
6

Thanks ! Now i got it.

7

I see, I thought you were saying that you would expect any SFE to be negative.
Your lambda steps are ok, you might not even need to go that fine. A step of 0.1 should be enough, except maybe between lambda=0.6-0.9. For the coulomb lambdas even 0.2 might be enough. You should always check that your lambda windows have sufficient overlap, which can vary quite a bit between systems. A good indicator is the standard deviation given for the individual lambda windows.
And yes, you can just delete the other rows.

8

No, but i understand the confusion, should have specified that. Okay, thanks ! Will try this run with these settings:

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

Is that just about what you meant ? Will report if it worked out !