Free Energy Calculation general questions and LINCS Warning

User discussions
1

GROMACS version: 2018
GROMACS modification: No
Hello Gromacs-Users,

I’m still doing Free Energy Calculation and it is getting better. Anyway I have some general questions. I absolutley understand, when this questions should answered by paper-research. But I find it hard to get clear answers, so I hope it’s okay, when I try to find answers here.

Background: I try do find the octanol-water partition coefficient for some aroma compounds. I follow the Gromacs Free Energy Tutorial for this in general, but I use GROMOS96 54a7 an I plan to do an absolute Free Energy Calculation (vdw and coulomb), as soon as the current probems are solved.

My questions are:

  1. I have received documents from a paper author. They did a normal simulation until production run and only then they involved the Free-Energy mdp-options. In the tutorial the Free Energy options are involved from the energy minimization. Is both possible?
  2. The tutorial told me to set all charges to 0. Is this also allowed for the solvent (Octanol)?
  3. My simulation is a validation for finding the partition coefficient in a lipid and sugar system. I read, sometimes Gromacs won’t work for some Free Energy calculations. Is this the case for me?
  4. Did I understand the thermodynamic cycle for finding the partition coeffiecient: I do 2 simulations, each simulation is changing from state A to state B because of the different lambda steps?
2

There are many ways to prepare the initial state of a system. There is no “one size fits all” approach to any MD simulation, though I still think it wise to equilibrate under the new Hamiltonian for each value of λ

NO! The tutorial is not for solvation free energy, only for the vdW component. In reality, it isn’t even necessary to set charges to zero given the flexibility of the .mdp keywords. I did it to keep things simple and reproduce a very specific quantity. You should not set any charges to zero in the topology if you are computing a solvation free energy.

That seems like a weird and vague assertion, and without more detail on the specific criticism, no one can possibly comment on this.

You’ll need to do this in both water and octanol solvents. Note that octanol in a biphasic system also has some small amount of water (0.25 mol fraction or so, IIRC) so you may need to address this in your own simulations, though it may not matter within error of the computed ΔG.

1 Like
3

Wow, thank you! That was very helpful! There’s another question out of interest. In the Advanced Applications part of the Tutorial is said, for an absolute Free Energy Calculation I should set
couple-lambda0 = none
couple-lambda1 = vdw-q
I don’t understand the difference between couple-lambda0 and couple-lambda1. Why can’t I set couple-lambda0 = vdw-q?

4

In this process, a “dummy” entity becomes the fully interacting ligand as a function of λ, thus the resulting quantity is ΔG of solvation.

Here, a fully interacting species becomes a non-interacting dummy and the resulting quantity is ΔG of desolvation.

5

I understand! Thank you again!

6

Never ending questions: My octanol-water-partition-coefficient is for validation of another system. The aim is to find the partition coefficient of small molecules in a lipid/crystal sugar system. So I can do the Free Energy Calculation easily in Lipids. But is it working in a system, which is ‘solvated’ with little sugar crystals (recent finding, that solved sugar won’t work)? Another idea would have been, to put the small molecules on a sugar surface and try to do the Free Energy Calculation. But actually this is against the physical rules - because there’s actually no solvent. So I’m not sure if I can do a Free Energy Solvation with crystals. Actually I should use umbrella sampling here, shouldn’t I? (Hope my question is understandable)

7

Alchemical free energy calculations don’t require solvent. They require that the nonbonded interactions of one molecule be physically decoupled from whatever else is in the system. This is typically solvent, but does not need to be.

8

Your support is so great! Thanks a lot, that was a very important information!

9

My Free Energy Simulation worked now. But I’m struggling again.
I got an output with gmx bar for each simulation, one with water and one with octanol. I will post both outputs below. I have done everything the same, expect: I used ‘solvate’ for the water molecules and ‘insert-molecules’ for the octanol-molecules at the beginning. But for the octanol-simulation I got a 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.
I read about it and there is said, that I can’t have negative entropy. I understand that. But I have no clue, what I can do against it or what is the reason. Again I would appreciate help very much!

Water output:
lam_A lam_B DG +/- s_A +/- s_B +/- stdev +/-
0 1 0.29 0.01 0.04 0.01 0.04 0.01 0.27 0.00
1 2 0.90 0.01 0.07 0.01 0.08 0.01 0.40 0.00
2 3 1.54 0.01 0.15 0.02 0.17 0.02 0.58 0.00
3 4 2.01 0.06 0.38 0.04 0.41 0.04 0.90 0.02
4 5 1.88 0.06 0.77 0.07 0.70 0.07 1.25 0.03
5 6 1.30 0.03 0.78 0.05 0.61 0.04 1.21 0.03
6 7 0.77 0.03 0.65 0.04 0.50 0.04 1.09 0.02
7 8 0.29 0.01 0.57 0.02 0.44 0.02 0.97 0.01
8 9 -0.08 0.01 0.38 0.03 0.29 0.02 0.83 0.00
9 10 -0.34 0.03 0.33 0.02 0.25 0.01 0.73 0.02
10 11 -0.52 0.02 0.21 0.02 0.16 0.01 0.62 0.01
11 12 -0.61 0.01 0.18 0.01 0.14 0.01 0.57 0.01
12 13 -0.69 0.00 0.16 0.01 0.13 0.00 0.53 0.00
13 14 -0.76 0.01 0.14 0.00 0.12 0.00 0.49 0.01
14 15 -0.82 0.00 0.12 0.00 0.10 0.00 0.46 0.00
15 16 -0.86 0.00 0.12 0.00 0.10 0.00 0.44 0.00
16 17 -0.90 0.00 0.09 0.00 0.08 0.00 0.42 0.00
17 18 -0.92 0.01 0.09 0.01 0.08 0.01 0.41 0.00
18 19 -0.95 0.01 0.09 0.01 0.07 0.01 0.39 0.00
19 20 -0.97 0.00 0.08 0.00 0.07 0.00 0.38 0.00
21 22 -0.01 0.00 0.01 0.00 0.01 0.00 0.14 0.00
22 23 -0.03 0.00 0.01 0.00 0.01 0.00 0.14 0.00
23 24 -0.05 0.00 0.01 0.00 0.01 0.00 0.14 0.00
24 25 -0.07 0.00 0.01 0.00 0.01 0.00 0.15 0.00
25 26 -0.09 0.00 0.01 0.00 0.01 0.00 0.15 0.00
26 27 -0.11 0.00 0.01 0.00 0.01 0.00 0.15 0.00
27 28 -0.13 0.00 0.01 0.00 0.01 0.00 0.16 0.00
28 29 -0.16 0.00 0.01 0.00 0.01 0.00 0.16 0.00
29 30 -0.19 0.00 0.01 0.00 0.01 0.00 0.17 0.00
30 31 -0.21 0.00 0.01 0.00 0.01 0.00 0.17 0.00
31 32 -0.24 0.00 0.02 0.00 0.02 0.00 0.18 0.00
32 33 -0.28 0.00 0.01 0.00 0.01 0.00 0.18 0.00
33 34 -0.30 0.00 0.01 0.00 0.01 0.00 0.18 0.00
34 35 -0.34 0.00 0.02 0.00 0.02 0.00 0.19 0.00
35 36 -0.38 0.01 0.02 0.01 0.02 0.01 0.20 0.00
36 37 -0.42 0.01 0.02 0.01 0.02 0.01 0.20 0.00
37 38 -0.46 0.00 0.02 0.01 0.02 0.00 0.21 0.00
38 39 -0.51 0.00 0.02 0.00 0.02 0.00 0.21 0.00
39 40 -0.56 0.00 0.03 0.00 0.03 0.00 0.22 0.00
40 41 -0.61 0.00 0.02 0.00 0.02 0.00 0.22 0.00

Final results in kJ/mol:

point 0 - 1, DG 0.72 +/- 0.01
point 1 - 2, DG 2.24 +/- 0.01
point 2 - 3, DG 3.82 +/- 0.03
point 3 - 4, DG 4.98 +/- 0.14
point 4 - 5, DG 4.66 +/- 0.15
point 5 - 6, DG 3.23 +/- 0.06
point 6 - 7, DG 1.90 +/- 0.08
point 7 - 8, DG 0.71 +/- 0.02
point 8 - 9, DG -0.21 +/- 0.03
point 9 - 10, DG -0.84 +/- 0.08
point 10 - 11, DG -1.28 +/- 0.05
point 11 - 12, DG -1.50 +/- 0.02
point 12 - 13, DG -1.70 +/- 0.01
point 13 - 14, DG -1.89 +/- 0.01
point 14 - 15, DG -2.03 +/- 0.00
point 15 - 16, DG -2.14 +/- 0.01
point 16 - 17, DG -2.24 +/- 0.01
point 17 - 18, DG -2.29 +/- 0.02
point 18 - 19, DG -2.35 +/- 0.02
point 19 - 20, DG -2.41 +/- 0.00
point 21 - 22, DG -0.01 +/- 0.00
point 22 - 23, DG -0.07 +/- 0.00
point 23 - 24, DG -0.11 +/- 0.00
point 24 - 25, DG -0.17 +/- 0.00
point 25 - 26, DG -0.22 +/- 0.00
point 26 - 27, DG -0.27 +/- 0.01
point 27 - 28, DG -0.33 +/- 0.00
point 28 - 29, DG -0.39 +/- 0.01
point 29 - 30, DG -0.46 +/- 0.00
point 30 - 31, DG -0.52 +/- 0.00
point 31 - 32, DG -0.60 +/- 0.00
point 32 - 33, DG -0.68 +/- 0.01
point 33 - 34, DG -0.75 +/- 0.00
point 34 - 35, DG -0.83 +/- 0.01
point 35 - 36, DG -0.94 +/- 0.01
point 36 - 37, DG -1.04 +/- 0.02
point 37 - 38, DG -1.14 +/- 0.01
point 38 - 39, DG -1.25 +/- 0.01
point 39 - 40, DG -1.39 +/- 0.01
point 40 - 41, DG -1.51 +/- 0.01

total 0 - 41, DG -11.34 +/- 0.30

Octanol output:
lam_A lam_B DG +/- s_A +/- s_B +/- stdev +/-
0 1 0.27 0.01 0.14 0.02 0.15 0.02 0.51 0.01
1 2 0.66 0.02 0.29 0.04 0.31 0.04 0.76 0.02
2 3 0.75 0.04 0.46 0.05 0.44 0.05 1.01 0.02
3 4 0.52 0.07 0.60 0.08 0.49 0.07 1.11 0.01
4 5 0.12 0.08 0.62 0.12 0.48 0.10 1.05 0.06
5 6 -0.27 0.08 0.47 0.04 0.36 0.04 0.92 0.03
6 7 -0.59 0.05 0.43 0.05 0.34 0.04 0.82 0.03
7 8 -0.83 0.02 0.24 0.03 0.18 0.02 0.68 0.01
8 9 -0.92 0.01 0.20 0.03 0.15 0.02 0.64 0.01
9 10 -0.98 0.02 0.21 0.02 0.16 0.02 0.62 0.01
10 11 -1.08 0.01 0.19 0.01 0.15 0.01 0.58 0.01
11 12 -1.16 0.01 0.18 0.02 0.14 0.01 0.55 0.01
12 13 -1.22 0.01 0.12 0.02 0.09 0.01 0.50 0.01
13 14 -1.25 0.01 0.14 0.01 0.12 0.01 0.50 0.01
14 15 -1.28 0.01 0.10 0.01 0.08 0.01 0.47 0.01
15 16 -1.29 0.01 0.11 0.01 0.09 0.01 0.45 0.01
16 17 -1.32 0.01 0.11 0.01 0.09 0.01 0.43 0.01
17 18 -1.35 0.01 0.10 0.00 0.08 0.00 0.42 0.01
18 19 -1.38 0.01 0.09 0.01 0.08 0.01 0.41 0.01
19 20 -1.40 0.01 0.09 0.00 0.08 0.00 0.40 0.01
21 22 -0.01 0.00 0.00 0.00 0.00 0.00 0.05 0.01
22 23 -0.01 0.00 -0.00 0.00 -0.00 0.00 0.04 0.01
23 24 -0.01 0.00 0.00 0.00 0.00 0.00 0.04 0.00
24 25 -0.01 0.00 -0.00 0.00 -0.00 0.00 0.05 0.00
25 26 -0.01 0.00 0.00 0.00 0.00 0.00 0.05 0.01
26 27 -0.02 0.00 0.00 0.00 0.00 0.00 0.05 0.01
27 28 -0.02 0.01 -0.00 0.00 -0.00 0.00 0.05 0.01
28 29 -0.02 0.00 0.00 0.00 0.00 0.00 0.05 0.00
29 30 -0.03 0.01 0.00 0.01 0.00 0.01 0.05 0.01
30 31 -0.03 0.01 -0.00 0.01 -0.00 0.01 0.05 0.00
31 32 -0.03 0.01 0.00 0.00 0.00 0.00 0.06 0.00
32 33 -0.03 0.01 -0.01 0.00 -0.01 0.00 0.06 0.00
33 34 -0.03 0.00 0.00 0.00 0.00 0.00 0.04 0.01
34 35 -0.04 0.01 0.01 0.00 0.01 0.00 0.05 0.00
35 36 -0.05 0.01 -0.00 0.01 -0.00 0.01 0.05 0.00
36 37 -0.04 0.01 -0.01 0.00 -0.01 0.00 0.05 0.01
37 38 -0.06 0.01 0.03 0.00 0.03 0.00 0.05 0.00
38 39 -0.07 0.01 -0.01 0.01 -0.01 0.01 0.06 0.00
39 40 -0.06 0.01 -0.00 0.00 -0.00 0.00 0.07 0.01
40 41 -0.04 0.01 -0.01 0.01 -0.01 0.01 0.07 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 0.66 +/- 0.01
point 1 - 2, DG 1.64 +/- 0.04
point 2 - 3, DG 1.86 +/- 0.09
point 3 - 4, DG 1.29 +/- 0.18
point 4 - 5, DG 0.30 +/- 0.20
point 5 - 6, DG -0.67 +/- 0.19
point 6 - 7, DG -1.47 +/- 0.12
point 7 - 8, DG -2.06 +/- 0.05
point 8 - 9, DG -2.28 +/- 0.04
point 9 - 10, DG -2.42 +/- 0.05
point 10 - 11, DG -2.67 +/- 0.03
point 11 - 12, DG -2.88 +/- 0.03
point 12 - 13, DG -3.02 +/- 0.03
point 13 - 14, DG -3.10 +/- 0.03
point 14 - 15, DG -3.17 +/- 0.01
point 15 - 16, DG -3.20 +/- 0.02
point 16 - 17, DG -3.27 +/- 0.01
point 17 - 18, DG -3.35 +/- 0.02
point 18 - 19, DG -3.41 +/- 0.02
point 19 - 20, DG -3.48 +/- 0.03
point 21 - 22, DG -0.01 +/- 0.01
point 22 - 23, DG -0.02 +/- 0.00
point 23 - 24, DG -0.03 +/- 0.00
point 24 - 25, DG -0.03 +/- 0.00
point 25 - 26, DG -0.03 +/- 0.01
point 26 - 27, DG -0.05 +/- 0.01
point 27 - 28, DG -0.05 +/- 0.01
point 28 - 29, DG -0.05 +/- 0.01
point 29 - 30, DG -0.07 +/- 0.01
point 30 - 31, DG -0.08 +/- 0.02
point 31 - 32, DG -0.08 +/- 0.02
point 32 - 33, DG -0.08 +/- 0.02
point 33 - 34, DG -0.06 +/- 0.01
point 34 - 35, DG -0.09 +/- 0.01
point 35 - 36, DG -0.12 +/- 0.02
point 36 - 37, DG -0.09 +/- 0.02
point 37 - 38, DG -0.14 +/- 0.02
point 38 - 39, DG -0.19 +/- 0.02
point 39 - 40, DG -0.14 +/- 0.03
point 40 - 41, DG -0.11 +/- 0.02

total 0 - 41, DG -36.23 +/- 0.45

If needed, Production mdp:
; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 500000 ; 1 ns
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.4
; Electrostatics
coulombtype = PME
rcoulomb = 1.4
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.3
rvdw = 1.4
; 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 = tetra ; name of moleculetype to decouple
couple-lambda0 = none
couple-lambda1 = vdw-q
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 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
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 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
coul_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.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 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 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 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 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 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 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 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 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 = all-bonds
; Type of constraint algorithm
constraint-algorithm = lincs
; Constrain the starting configuration
; since we are continuing from NPT
continuation = yes
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12

10

I have similar question: should we have couple-intramol = yes, since octanol is not a small molecules, for my case I am coupling vdw and coulomb as follows:
couple-lambda0 = none
couple-lambda1 = vdw-q
couple-intramol = yes

; 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
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 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
coul_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.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Am I doing it correctly? first increasing LJ and then coulomb coupling, is it safe to have 0.05 step for vdw and 0.1 for coulomb, I read in shirts et al. that one should also need uneven steps, but to start with is it ok to use this step sizes?