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

Sander Pronk pronk at kth.se
Wed Nov 14 11:08:26 CET 2012


On Nov 14, 2012, at 08:42 , David van der Spoel <spoel at xray.bmc.uu.se> wrote:

> 
> 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.

This simply means that g_bar is integrating dH/dlambda values, not using BAR (the preferred method). 
g_bar then has to assume that the dH/dlambda values it gets from the dhdl.xvg file are the result of a standard linear decoupling - if you just set the lambda values like you did the coupling is linear.

Most likely, not all lambda points had 'foreign_lambda' values for all neighboring lambda points. g_bar outputs all lambda values it finds in all dhdl.xvg files. 

It might also be that you're not giving g_bar all of the dhdl.xvg files.

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

no - it just means that you shouldn't set anything manually at the individual points that couples to lambda in a non-linear way.


> ==> 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.

This is triggered when the shared entropy in one or the other direction is negative: that is an unphysical result that can be the result of very few sampling points leading to large fluctuations.


> 
> 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
> 
> 
> -- 
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-developers-request at gromacs.org.




More information about the gromacs.org_gmx-developers mailing list