[gmx-developers] Error of violate the Second Law of Thermodynamics in Free energy calculation with BAR

[David van der Spoel]

Dear GMX developers

We're calculating some small organic molecule's desolvation free energy. 
Recently we got this error from BAR calculation. Please anyone explain 
what's wrong in here?

First:
WARNING: Using the derivative data (dH/dlambda) to extrapolate delta H 
values.
This will only work if the Hamiltonian is linear in lambda.

==> Does this mean that when using softcore, or unevenly spaced lambda 
values this approximation will break down?

==> and can this cause the second warning below?

and then Second:
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 have two MD simulation steps. For example desolvation free energy of 
chloroform and methanol:
The final result from MD1 something like this:
total   0.000 -  1.000,   DG  7.72 +/-  0.05

The final result from MD2 something like this:
total   0.000 -  1.000,   DG  3.96 +/-  0.06

Total Gibbs energy of desolvation 11.7 +/- 0.1 kJ/mol  (including with 
these two warnings)

Something wrong in MD?
MD1.mdp
;Run control
integrator               = sd       ; Langevin dynamics
tinit                    = 0
dt                       = 0.002
nsteps                   = 1000000  ; 2 ns
nstcomm                  = 100
; Output control
nstxout                  = 500
nstvout                  = 500
nstfout                  = 0
nstlog                   = 500
nstenergy                = 500
nstxtcout                = 0
xtc-precision            = 1000
; Neighborsearching and short-range nonbonded interactions
nstlist                  = 10
ns_type                  = grid
pbc                      = xyz
rlist                    = 1.0
; Electrostatics
coulombtype              = PME
rcoulomb                 = 1.0
; van der Waals
vdw-type                 = switch
rvdw-switch              = 0.8
rvdw                     = 0.9
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                  = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.12
; EWALD/PME/PPPM parameters
pme_order                = 6
ewald_rtol               = 1e-06
epsilon_surface          = 0
optimize_fft             = no
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps                  = system
tau_t                    = 0.2
ref_t                    = 300
; Pressure coupling is on for NPT
Pcoupl                   = Parrinello-Rahman
tau_p                    = 5
compressibility          = 4.5e-05
ref_p                    = 1.0
; Free energy control stuff
free_energy              = yes
init_lambda              = 0.0
delta_lambda             = 0
foreign_lambda           = 0.05
sc-alpha                 = 0.5
sc-power                 = 1.0
sc-sigma                 = 0.3
couple-lambda0           = vdw-q      ; only van der Waals interactions
couple-lambda1           = vdw     ; turn off everything, in this case 
only vdW
couple-intramol          = no
nstdhdl                  = 10
; Do not generate velocities
gen_vel                  = no
; options for bonds
constraints              = all-bonds  ;
; Type of constraint algorithm
constraint-algorithm     = lincs
; Constrain the starting configuration
; since we are continuing from NPT
continuation             = yes
; Highest order in the expansion of the constraint coupling matrix
lincs-order              = 12

MD2.mdp
;Run control
integrator               = sd       ; Langevin dynamics
tinit                    = 0
dt                       = 0.002
nsteps                   = 1000000  ; 2 ns
nstcomm                  = 100
; Output control
nstxout                  = 500
nstvout                  = 500
nstfout                  = 0
nstlog                   = 500
nstenergy                = 500
nstxtcout                = 0
xtc-precision            = 1000
; Neighborsearching and short-range nonbonded interactions
nstlist                  = 10
ns_type                  = grid
pbc                      = xyz
rlist                    = 1.0
; Electrostatics
coulombtype              = PME
rcoulomb                 = 1.0
; van der Waals
vdw-type                 = switch
rvdw-switch              = 0.8
rvdw                     = 0.9
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                  = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.12
; EWALD/PME/PPPM parameters
pme_order                = 6
ewald_rtol               = 1e-06
epsilon_surface          = 0
optimize_fft             = no
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps                  = system
tau_t                    = 0.2
ref_t                    = 300
; Pressure coupling is on for NPT
Pcoupl                   = Parrinello-Rahman
tau_p                    = 5
compressibility          = 4.5e-05
ref_p                    = 1.0
; Free energy control stuff
free_energy              = yes
init_lambda              = 0.0
delta_lambda             = 0
foreign_lambda           = 0.05
sc-alpha                 = 0.5
sc-power                 = 1.0
sc-sigma                 = 0.3
couple-lambda0           = vdw      ; only van der Waals interactions
couple-lambda1           = none     ; turn off everything, in this case 
only vdW
couple-intramol          = no
nstdhdl                  = 10
; Do not generate velocities
gen_vel                  = no
; options for bonds
constraints              = all-bonds  ;
; Type of constraint algorithm
constraint-algorithm     = lincs
; Constrain the starting configuration
; since we are continuing from NPT
continuation             = yes
; Highest order in the expansion of the constraint coupling matrix
lincs-order              = 12

Thank you so much.

Khatnaa
-- 
gmx-users mailing list    gmx-users at gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists