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.gil9...@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.loeff...@stfc.ac.uk > > wrote: > > > On Tue, 16 May 2017 15:13:10 -0400 > > Dan Gil <dan.gil9...@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-requ...@gromacs.org. > > -- 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-requ...@gromacs.org.