[gmx-users] g_bar ... the right answer for the wrong reason?

[Michael Brunsteiner]

Hi,
i just performed a free energy TI calculation, to get the
free energy of solvation of water in water (the chemical potential of water)
i stuck closely to the templates given in the tutoral by justin lemkul.
The final result i get with g_bar looks good, and the number is within the error-bars
of both experiment and literature data - still i keep getting a warning ...
i first turn off the charges and then LJ ... but for the LJ part (only) g_bar keeps saying:

WARNING: Some of these results violate the Second Law of Thermodynamics: 
         This is can be the result of severe undersampling, or (more likely)
         there is something wrong with the simulations.

i can't think of anything that's wrong with my simulations, nor do i think
this is undersampling (21 windows each 2 nano seconds)
also the difference between my mdp and justin's are minor (as i believe, see below)
why do i get the correct answer but still this warning keeps coming?

at lambda=0 my LJ interactions are fully turned on, and for the DeltaG(0->0.05)
i,e, the first row in the table below i get from g_bar:

 lam_A  lam_B      DG   +/-     s_A   +/-     s_B   +/-   stdev   +/- 
 0.000  0.050   -2.90  0.03    2.15  0.03 9647088805553410.00 9570465417587940.00    3.30  0.05
 0.050  0.100   -3.36  0.03    1.48  0.02    1.47  0.03    1.95  0.02
 0.100  0.150   -0.49  0.00    0.03  0.00    0.03  0.00    0.22  0.00
 0.150  0.200    0.12  0.00    0.00  0.00    0.00  0.00    0.05  0.00
etc...

there appears to be a problem with the phase space overlap, but then
all other values, and the end-result, look ok ... can it be that there is an issue
with the gmx implementations of soft core potentials and g_bar?


cheers
mic


mdp:
;
integrator               = sd
nsteps                   = 1000000
;dt                       = 0.002

pbc                      = xyz
nstcomm                  = 100
comm_grps                = System
;
nstxtcout                = 0
nstenergy                = 100
nst_log                  = 100
nstxout                  = 0
nstvout                  = 0
;
rlist                    = 1.0
;
coulombtype              =
PME
rcoulomb                 = 1.0
vdw-type                 = switch
rvdw-switch              = 0.8
rvdw                     = 0.9
;
DispCorr                 = EnerPres
;
fourierspacing           = 0.12
pme_order                = 6
ewald_rtol               = 1e-06
epsilon_surface        = 0
optimize_fft             = no
;
tc_grps                  = system
tau_t                    = 1.0
ref_t                    = 300
;
Pcoupl                   = Berendsen
tau_p                    = 0.5
compressibility          = 4.5e-5
ref_p                    = 1.0
pcoupltype               = isotropic
;
constraints              = hbonds
lincs_iter               = 8
;
free_energy              = yes
init_lambda              = 0.00
delta_lambda             = 0.0
foreign_lambda           = 0.05 0.10
dhdl_derivatives         = yes
sc_alpha                 = 0.5
sc_power                 = 1
sc_sigma                 = 0.5   ; i also tried 0.3 same problem
nstdhdl                  = 100
couple-moltype           = w1
couple-lambda0           = vdw
couple-lambda1           = none
couple-intramol          = no




===============================
Why be happy when you could be normal?