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.