[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