[gmx-users] Simulations showing that mass contributes to the solvation free energy in a thermodynamic cycle
Dan Gil
dan.gil9973 at gmail.com
Fri Aug 11 17:31:25 CEST 2017
Hi,
I have been attempting TI calculations for a while, and here I am again. I
have asked questions about this here before, but I wasn't able to find a
satisfactory solution :(. I am determined to figure this problem out though.
Here is the thermodynamic cycle:
Step A. [Heptane to Perfluoroheptane] in an ideal gas state
Step B. [Heptane to Perfluoroheptane] in an aqueous solution
Step C. [Heptane] in an ideal gas state to [Heptane] in an aqueous solution
Step D. [Perfluoroheptane] in an ideal gas state to [Perfluoroheptane] in
an aqueous solution.
I am doing step A and B in simulation to find the solvation free energy
difference of heptane and perfluoroheptane in water.
I prepared the input configuration by inserting 7833 SPC/E water molecules
around one OPLS-AA heptane molecule. I equilibrated with the PR barostat (1
bar) and obtained the average volume and box dimensions. The topology file
for heptane/perfluoroheptane used are at the bottom of this message.
All thermodynamic integrations were done with an NVT ensemble. The mdp file
is also listed at the bottom of this message.
I used gmx bar to analyze the results. The result is (somehow) that the
mass of the molecule largely contributes to the solvation energy.
I've also tried playing around with other options. I used all-bond
constraints (against warnings), without DispCorr, with tau-t = 0.1, with
NPT ensemble, using sd integrator versus md, varying more than one lambda
at a time, using an older version of Gromacs (4.6 vs 5.1), and other things
as well. Every time, the mass had large contributions to the solvation
energy.
I plotted the output of gmx bar that I show below. Both the derivative and
the actual free energy follow a rather smooth curve. The plot of Step A is
consistently higher than Step B, which makes me think that this unlikely to
be a sampling error.
I also did simulations with bonded-lambdas changing only. Here, the free
energy change canceled out almost exactly, which is very surprising because
it goes against published literature results. Bonded interactions are
supposed to have the largest effect on the solvation free energy.
I've detailed out everything I can, and I really appreciate that you've
taken the time to read this far! Please let me know if I am missing
something.
Dan
Step A:
lam_A lam_B DG +/- s_A +/- s_B +/- stdev +/-
0 1 26.08 0.52 18.97 0.17 10.91 0.49 21.50 2.04
1 2 11.58 0.11 3.59 0.12 2.57 0.05 2.93 0.06
2 3 7.60 0.06 1.41 0.08 1.14 0.08 1.65 0.04
3 4 5.68 0.04 0.79 0.04 0.62 0.05 1.22 0.02
4 5 4.52 0.05 0.54 0.02 0.43 0.04 0.97 0.01
5 6 3.74 0.03 0.35 0.05 0.27 0.03 0.77 0.02
6 7 3.22 0.03 0.26 0.04 0.20 0.02 0.67 0.01
7 8 2.81 0.02 0.21 0.02 0.15 0.04 0.60 0.01
8 9 2.49 0.02 0.17 0.03 0.14 0.02 0.53 0.01
9 10 2.23 0.01 0.12 0.01 0.09 0.03 0.47 0.01
Final results in kJ/mol:
point 0 - 1, DG 65.06 +/- 1.30
point 1 - 2, DG 28.88 +/- 0.27
point 2 - 3, DG 18.97 +/- 0.16
point 3 - 4, DG 14.16 +/- 0.11
point 4 - 5, DG 11.27 +/- 0.13
point 5 - 6, DG 9.34 +/- 0.07
point 6 - 7, DG 8.02 +/- 0.08
point 7 - 8, DG 7.01 +/- 0.06
point 8 - 9, DG 6.21 +/- 0.04
point 9 - 10, DG 5.57 +/- 0.02
total 0 - 10, DG 174.46 +/- 1.51
Step B:
lam_A lam_B DG +/- s_A +/- s_B +/- stdev +/-
0 1 15.71 1.35 11.89 2.28 5.88 0.87 5.71 2.90
1 2 8.87 0.67 0.95 0.54 0.40 0.50 1.97 0.14
2 3 6.10 0.21 2.37 0.36 1.70 0.26 2.16 0.25
3 4 4.72 0.14 -0.32 0.09 -0.38 0.08 1.08 0.04
4 5 4.52 0.09 0.58 0.11 0.52 0.11 0.98 0.05
5 6 2.97 0.12 1.03 0.10 0.95 0.09 1.17 0.05
6 7 2.29 0.22 -0.28 0.16 -0.33 0.17 0.76 0.09
7 8 2.59 0.20 0.03 0.15 0.02 0.15 0.71 0.03
8 9 2.28 0.05 0.29 0.10 0.27 0.09 0.60 0.05
9 10 1.95 0.04 0.06 0.10 0.06 0.11 0.53 0.03
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 39.18 +/- 3.37
point 1 - 2, DG 22.13 +/- 1.67
point 2 - 3, DG 15.23 +/- 0.51
point 3 - 4, DG 11.77 +/- 0.35
point 4 - 5, DG 11.27 +/- 0.21
point 5 - 6, DG 7.40 +/- 0.31
point 6 - 7, DG 5.72 +/- 0.55
point 7 - 8, DG 6.46 +/- 0.49
point 8 - 9, DG 5.69 +/- 0.12
point 9 - 10, DG 4.85 +/- 0.10
total 0 - 10, DG 129.70 +/- 5.91
TOPOLOGY:
[ moleculetype ]
HEPT 3
[ atoms ]
1 opls_135 1 HEPT CH3 1 -0.180 12.011
opls_961 0.360 12.011
2 opls_136 1 HEPT CH2 2 -0.120 12.011
opls_962 0.240 12.011
3 opls_136 1 HEPT CH2 3 -0.120 12.011
opls_962 0.240 12.011
4 opls_136 1 HEPT CH2 4 -0.120 12.011
opls_962 0.240 12.011
5 opls_136 1 HEPT CH2 5 -0.120 12.011
opls_962 0.240 12.011
6 opls_136 1 HEPT CH2 6 -0.120 12.011
opls_962 0.240 12.011
7 opls_135 1 HEPT CH3 7 -0.180 12.011
opls_961 0.360 12.011
8 opls_140 1 HEPT H 1 0.060 1.008
opls_965 -0.120 18.998
9 opls_140 1 HEPT H 1 0.060 1.008
opls_965 -0.120 18.998
10 opls_140 1 HEPT H 1 0.060 1.008
opls_965 -0.120 18.998
11 opls_140 1 HEPT H 2 0.060 1.008
opls_965 -0.120 18.998
12 opls_140 1 HEPT H 2 0.060 1.008
opls_965 -0.120 18.998
13 opls_140 1 HEPT H 3 0.060 1.008
opls_965 -0.120 18.998
14 opls_140 1 HEPT H 3 0.060 1.008
opls_965 -0.120 18.998
15 opls_140 1 HEPT H 4 0.060 1.008
opls_965 -0.120 18.998
16 opls_140 1 HEPT H 4 0.060 1.008
opls_965 -0.120 18.998
17 opls_140 1 HEPT H 5 0.060 1.008
opls_965 -0.120 18.998
18 opls_140 1 HEPT H 5 0.060 1.008
opls_965 -0.120 18.998
19 opls_140 1 HEPT H 6 0.060 1.008
opls_965 -0.120 18.998
20 opls_140 1 HEPT H 6 0.060 1.008
opls_965 -0.120 18.998
21 opls_140 1 HEPT H 7 0.060 1.008
opls_965 -0.120 18.998
22 opls_140 1 HEPT H 7 0.060 1.008
opls_965 -0.120 18.998
23 opls_140 1 HEPT H 7 0.060 1.008
opls_965 -0.120 18.998
[ bonds ]
1 2 1
1 8 1
1 9 1
1 10 1
2 3 1
2 11 1
2 12 1
3 4 1
3 13 1
3 14 1
4 5 1
4 15 1
4 16 1
5 6 1
5 17 1
5 18 1
6 7 1
6 19 1
6 20 1
7 21 1
7 22 1
7 23 1
[ pairs ]
1 4 1
1 13 1
1 14 1
2 5 1
2 15 1
2 16 1
3 8 1
3 9 1
3 10 1
3 6 1
3 17 1
3 18 1
4 11 1
4 12 1
4 7 1
4 19 1
4 20 1
5 13 1
5 14 1
5 21 1
5 22 1
5 23 1
6 15 1
6 16 1
7 17 1
7 18 1
8 11 1
8 12 1
9 11 1
9 12 1
10 11 1
10 12 1
11 13 1
11 14 1
12 13 1
12 14 1
13 15 1
13 16 1
14 15 1
14 16 1
15 17 1
15 18 1
16 17 1
16 18 1
17 19 1
17 20 1
18 19 1
18 20 1
19 21 1
19 22 1
19 23 1
20 21 1
20 22 1
20 23 1
[ angles ]
1 2 3 1
1 2 11 1
1 2 12 1
2 1 8 1
2 1 9 1
2 1 10 1
2 3 4 1
2 3 13 1
2 3 14 1
3 2 11 1
3 2 12 1
3 4 5 1
3 4 15 1
3 4 16 1
4 3 13 1
4 3 14 1
4 5 6 1
4 5 17 1
4 5 18 1
5 4 15 1
5 4 16 1
5 6 7 1
5 6 19 1
5 6 20 1
6 5 17 1
6 5 18 1
6 7 21 1
6 7 22 1
6 7 23 1
7 6 19 1
7 6 20 1
8 1 9 1
8 1 10 1
9 1 10 1
11 2 12 1
13 3 14 1
15 4 16 1
17 5 18 1
19 6 20 1
21 7 22 1
21 7 23 1
22 7 23 1
[ dihedrals ]
1 2 3 4 3
1 2 3 13 3
1 2 3 14 3
2 3 4 5 3
2 3 4 15 3
2 3 4 16 3
3 2 1 8 3
3 2 1 9 3
3 2 1 10 3
3 4 5 6 3
3 4 5 17 3
3 4 5 18 3
4 3 2 11 3
4 3 2 12 3
4 5 6 7 3
4 5 6 19 3
4 5 6 20 3
5 4 3 13 3
5 4 3 14 3
5 6 7 21 3
5 6 7 22 3
5 6 7 23 3
6 5 4 15 3
6 5 4 16 3
7 6 5 17 3
7 6 5 18 3
8 1 2 11 3
8 1 2 12 3
9 1 2 11 3
9 1 2 12 3
10 1 2 11 3
10 1 2 12 3
11 2 3 13 3
11 2 3 14 3
12 2 3 13 3
12 2 3 14 3
13 3 4 15 3
13 3 4 16 3
14 3 4 15 3
14 3 4 16 3
15 4 5 17 3
15 4 5 18 3
16 4 5 17 3
16 4 5 18 3
17 5 6 19 3
17 5 6 20 3
18 5 6 19 3
18 5 6 20 3
19 6 7 21 3
19 6 7 22 3
19 6 7 23 3
20 6 7 21 3
20 6 7 22 3
20 6 7 23 3
MDP:
;Integration Method and Parameters
integrator = md
nsteps = 100000
dt = 0.002
nstenergy = 100
nstcalcenergy = 100
nstlog = 5000
;Output Control
nstxout = 100
nstvout = 0
;Cutoff Schemes
cutoff-scheme = group
rlist = 1.0
vdw-type = cut-off
rvdw = 1.0
DispCorr = EnerPres
;Coulomb interactions
coulombtype = pme
rcoulomb = 1.0
fourierspacing = 0.4
;Constraints
constraints = none
;Temperature coupling
tcoupl = v-rescale
tc-grps = System
tau-t = 0
ref-t = 300
;Free energy calculation
free-energy = yes
init-lambda-state = INITLAMBDA
delta-lambda = 0
calc-lambda-neighbors = 1
coul-lambdas = 0 0 0 0 0 0 0 0 0 0 0
vdw-lambdas = 0 0 0 0 0 0 0 0 0 0 0
bonded-lambdas = 0 0 0 0 0 0 0 0 0 0 0
restraint-lambdas = 0 0 0 0 0 0 0 0 0 0 0
mass-lambdas = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
nstdhdl = 100
More information about the gromacs.org_gmx-users
mailing list