[gmx-users] Perturbation Thermodynamic Integration
Hannes Loeffler
Hannes.Loeffler at stfc.ac.uk
Fri Jun 2 12:16:34 CEST 2017
What I notice is:
You don't use bond-constraints but you have set your time step to 2 fs.
What happens if you (significantly) increase simulation time per
lambda? Your variances are really large.
What is the motivation of using such a long cutoff radius for vdW?
On Fri, 26 May 2017 11:13:08 -0400
Dan Gil <dan.gil9973 at gmail.com> wrote:
> 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