# [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