[gmx-users] Perturbation Thermodynamic Integration
Dan Gil
dan.gil9973 at gmail.com
Fri May 26 17:13:13 CEST 2017
Hello,
I've run Gromacs following the suggestions, but I still get nonzero
contributions from the mass changing. I will outline in-depth what I did.
Thank you very much for your help!
(1) Heptane molecule was inserted into a box of size 10x10x10 nm3 with 7833
water molecules using Gromacs insert-molecules. Bad contacts were resolved
using the steep algorithm. The pressure was increased to a 10 bar using the
Berendsen pressure couple until a condensed liquid phase was formed. The
system was equilibrated at standard conditions (1 bar, 300 K) using the
Parinello-Rahman couple and velocity-rescaling thermostat. Step size of 2
femtoseconds with the stochastic integrator was used.
(2) New directories were created with copies of the configuration generated
in step (1). Each directory also contained grompp.mdp which defined the
lambda value (i.e. 1 through 10) for each directory. The configurations
were further equilibrated at this stage for 2 ns. The topology is shown
below. The topology is defined so that "A" is a hydrocarbon and "B" is a
fluorocarbon.
(3) Data collection simulations began after step (2), using grompp,mdp
shown below.
(4) Steps (2) and (3) were repeated for heptane in the gas phase with no
water molecules.
(5) Results: I've ran two alchemical transformation steps. A(aq)->B(aq) and
A(gas)->B(gas). The free energy change difference of the two steps is the
difference in solvation free energy. That is, deltaF(A(aq)->B(aq)) -
deltaF(A(gas)->B(gas)). Mass contribution is large. I will also leave a
matrix of the mean and variance of each contribution at the bottom.
--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
...
-- MDP Options --
;Integration Method and Parameters
integrator = sd
nsteps = 1000000
dt = 0.002
nstenergy = 1000
nstlog = 5000
;Output Control
nstxout = 0
nstvout = 0
;Cutoff Schemes
cutoff-scheme = group
rlist = 1.0
vdw-type = cut-off
rvdw = 2.0
;Coulomb interactions
coulombtype = pme
rcoulomb = 1.0
fourierspacing = 0.4
;Constraints
constraints = none
;Temperature coupling
tcoupl = v-rescale
tc-grps = system
tau-t = 0.1
ref-t = 300
;Pressure coupling
pcoupl = parrinello-rahman
ref-p = 1.01325
compressibility = 4.5e-5
tau-p = 5
;Free energy calculation
free-energy = yes
init-lambda-state = 0
delta-lambda = 0
fep-lambdas =
calc-lambda-neighbors = 1
vdw_lambdas = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
coul_lambdas = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
bonded_lambdas = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
restraint_lambdas = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
mass_lambdas = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
nstdhdl = 10
--Results--
A(aq)->B(aq)
Means
; mass-lambda | coul-lambda | vdw-lambda | bonded-lambda | restraint-lambda
| dH/dlambda to before | dH/dlambda to this | dH/dlambda to next | pV |
lambda
585.258 29.695 138.389 -106.111 0.000 0.000 0.000 78.180 14.434 0.0
194.905 27.926 109.373 -113.567 0.000 -8.514 0.000 35.428 14.433 0.1
201.785 23.966 90.522 -123.135 0.000 -5.856 0.000 32.987 14.437 0.2
101.067 21.244 76.286 -130.204 0.000 6.727 0.000 20.621 14.435 0.3
67.782 19.301 64.038 -136.090 0.000 12.170 0.000 15.391 14.436 0.4
54.102 17.403 54.579 -140.360 0.000 15.207 0.000 12.567 14.438 0.5
45.855 14.376 46.120 -146.966 0.000 17.949 0.000 10.041 14.437 0.6
39.769 11.953 39.144 -151.300 0.000 20.040 0.000 8.168 14.438 0.7
35.384 11.090 33.911 -151.862 0.000 21.252 0.000 7.172 14.439 0.8
31.411 9.077 28.113 -154.575 0.000 22.809 0.000 5.829 14.440 0.9
28.398 9.974 24.842 -151.347 0.000 23.126 0.000 0.000 14.439 1.0
Variances
; mass-lambda | coul-lambda | vdw-lambda | bonded-lambda | restraint-lambda
| dH/dlambda to before | dH/dlambda to this | dH/dlambda to next | pV |
lambda
17764.592 28.101 889.212 7624.759 0.000 0.000 0.000 265.656 0.002 0.0
1580.485 25.381 492.489 7085.318 0.000 93.157 0.000 95.406 0.002 0.1
4269.905 26.694 347.400 7025.386 0.000 118.370 0.000 120.552 0.002 0.2
1329.740 26.648 259.096 7083.151 0.000 85.371 0.000 87.603 0.002 0.3
249.921 24.406 195.233 7149.494 0.000 77.239 0.000 79.441 0.002 0.4
121.269 17.111 156.460 7039.120 0.000 74.030 0.000 76.206 0.002 0.5
87.000 15.901 128.267 7054.185 0.000 73.556 0.000 75.723 0.002 0.6
65.818 14.577 104.742 7096.497 0.000 73.225 0.000 75.402 0.002 0.7
55.143 14.713 91.124 7081.810 0.000 72.768 0.000 74.919 0.002 0.8
41.118 13.472 76.214 7198.105 0.000 73.496 0.000 75.675 0.002 0.9
33.836 15.572 65.528 7162.912 0.000 72.932 0.000 0.000 0.002 1.0
A(gas)->B(gas)
Means
; mass-lambda | coul-lambda | vdw-lambda | bonded-lambda | restraint-lambda
| dH/dlambda to before | dH/dlambda to this | dH/dlambda to next | pV |
lambda
1110.166 32.509 114.451 -115.748 0.000 0.000 0.000 127.586 4.377 0.0
398.592 26.708 101.772 -133.086 0.000 -26.055 0.000 52.957 3.912 0.1
237.629 27.141 84.628 -137.244 0.000 -7.766 0.000 34.879 5.057 0.2 84.897
20.358 76.144 -133.012 0.000 8.728 0.000 18.621 14.436 0.3 131.475 20.205
64.964 -158.091 0.000 7.811 0.000 19.736 3.927 0.4 107.720 20.655 57.763
-156.171 0.000 10.776 0.000 16.984 4.815 0.5 91.054 25.769 51.617 -146.044
0.000 11.632 0.000 16.326 4.020 0.6 78.003 21.697 45.907 -154.841 0.000
14.905 0.000 13.273 4.114 0.7 69.346 15.392 41.194 -167.694 0.000 18.277
0.000 10.139 6.244 0.8 62.263 17.129 37.528 -164.429 0.000 18.954 0.000
9.667 6.133 0.9 56.582 18.052 34.094 -162.977 0.000 19.727 0.000 0.000
4.456 1.0
Variances
; mass-lambda | coul-lambda | vdw-lambda | bonded-lambda | restraint-lambda
| dH/dlambda to before | dH/dlambda to this | dH/dlambda to next | pV |
lambda
40348.430 3.943 514.703 7085.063 0.000 0.000 0.000 481.126 0.030 0.0
6978.993 1.390 352.654 7972.572 0.000 151.629 0.000 154.128 0.000 0.1
2160.691 4.261 219.874 7020.561 0.000 94.026 0.000 96.219 0.263 0.2
301.074 32.826 249.228 7137.750 0.000 78.224 0.000 80.446 0.002 0.3
694.265 1.765 123.119 7225.511 0.000 80.863 0.000 83.117 0.000 0.4
452.408 7.903 97.154 7151.103 0.000 78.148 0.000 80.364 0.099 0.5
347.707 1.124 69.253 7024.693 0.000 74.678 0.000 76.844 0.004 0.6
230.115 4.723 56.425 7208.842 0.000 75.247 0.000 77.443 0.009 0.7
194.552 3.213 56.178 7169.195 0.000 74.497 0.000 76.669 2.824 0.8
154.875 4.498 47.516 7205.079 0.000 74.482 0.000 76.649 0.711 0.9
128.781 1.370 36.130 7231.535 0.000 74.006 0.000 0.000 0.155 1.0
On Tue, May 16, 2017 at 3:50 PM, Hannes Loeffler <Hannes.Loeffler at stfc.ac.uk
> wrote:
> On Tue, 16 May 2017 15:13:10 -0400
> Dan Gil <dan.gil9973 at gmail.com> wrote:
>
>
> > If you do this via decoupling ("absolute" transformation) you do that
> > > once for molecule A and once for molecule B.
> >
> >
> > I believe you are referring to the BAR method? I am trying to see if
> > I get the same results as another group. They used thermodynamic
> > integration from A to B, so I decided to spend some more time with
> > this. I will try later if I get consistent results for both methods.
>
> No, I am not referring to a specific analysis methods. I commented on
> an alternative pathway.
>
>
> > transforming the masses may interact badly
> > > with bond constraints that are applied to alchemically transformed
> > > bonds
> >
> >
> > Thank you for this information! Do you know if there are there
> > publications that refer to this issue?
>
> Not that I know of.
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>
More information about the gromacs.org_gmx-users
mailing list