[gmx-users] gmx 4.6.1, Expanded ensemble: weird balancing factors

[Joakim Jämbeck]

Dear Gromacs users,

I am trying to compute the free energy of hydration of ethane in order to get more familiar with the implementation of expanded ensemble in Gromacs v. 4.6.1.

The simulations runs fine and the balancing/weight factors are equilibrated after around 1.6 ns (a short equilibration time but this should be ok for a rough estimate of dG).

After this the simulation is continued for another 14 ns. At the end I look at the balancing factors to get a fast estimate of dG and see that they are tiny:

             MC-lambda information
  N  CoulL   VdwL    Count   G(in kT)  dG(in kT)
  1  0.000  0.000    16290    0.00000   -0.09268
  2  0.200  0.000    14940   -0.09268   -0.01647
  3  0.500  0.000    14627   -0.10914   -0.03867 <<
  4  1.000  0.000    14080   -0.14781    0.02986
  5  1.000  0.200    14309   -0.11795    0.07408
  6  1.000  0.400    15170   -0.04386    0.01711
  7  1.000  0.600    15276   -0.02675   -0.02304
  8  1.000  0.800    15061   -0.04979   -0.06110
  9  1.000  1.000    14328   -0.11090    0.00000

When adding these to get dG they value is way off what I get if I use the more traditional approach of separate simulations at each lambda value and using BAR at the end.

Does anyone have any ideas why I get these values?

Here is my mdp-file:

integrator         = sd
tinit                    = 0
dt                       = 0.001
nsteps              = 15000000
comm-mode   = Linear
nstcomm         = 10

nstlog               = 1000
nstenergy        = 100
nstlist                = 10
ns_type            = grid
pbc                    = xyz
rlist                     = 1.0
cutoff-scheme  = verlet

coulombtype    = PME
rcoulomb          = 1.0
vdw-type           = cutoff
rvdw                   = 1.0
fourierspacing  = 0.12
pme_order        = 4
ewald_rtol         = 1e-04
DispCorr            = EnerPres

tc-grps                 = System
tau_t                    = 5.0
ref_t                     = 300
Pcoupl                = berendsen
tau_p                   = 10.0
compressibility  = 4.5e-5
ref_p                    = 1.013
gen_vel               = no

constraints          = all-bonds
constraint-algorithm = lincs

; Free energy/expanded ensemble
free-energy            = expanded
sc-alpha                 = 0.5
sc-r-power              = 6
init-lambda-state   = 0
coul-lambdas         = 0.0 0.2 0.5 1.0 1.0 1.0 1.0 1.0 1.0
vdw-lambdas         = 0.0 0.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0
calc-lambda-neighbors = 9
nstdhdl                     = 200
dhdl-print-energy   = yes
couple-moltype       = C1X
couple-lambda0      = vdw-q
couple-lambda1      = none
couple-intramol        = no
nstexpanded             = 100
lmc-stats                    = wang-landau
lmc-move                   = metropolis
lmc-weights-equil     = wl-delta
weight-equil-wl-delta = 0.001
wl-scale                     = 0.7
wl-ratio                      = 0.8
init-wl-delta               = 1.0
wl-oneovert               = yes

Thanks in advance!

/ Joakim