[gmx-developers] Error of violate the Second Law of Thermodynamics in Free energy calculation with BAR
David van der Spoel
spoel at xray.bmc.uu.se
Fri Nov 16 16:31:30 CET 2012
On 2012-11-14 17:18, Sander Pronk wrote:
>
> On 14 Nov 2012, at 17:09 , David van der Spoel <spoel at xray.bmc.uu.se> wrote:
>
>> On 2012-11-14 15:23, Sander Pronk wrote:
>>> I just checked because I forgot, but g_bar does BAR on a linear interpolation of dH/dl if presented with only dH/dl data.
>>> There are legitimate situations where you'd want to do this, but they all obviously rely on the Hamiltonian depending linearly on lambda - hence the warning.
>>>
>> Sorry for ignorance, but is the way around this just to give it the edr files instead of the xvg files?
>
> No - the BAR method relies on Hamiltonian differences between the (lambda) state that you simulate at, and the state where you want the free energy difference to. Normally, those differences are calculated directly through the foreign_lambda mdp setting.
> Suppose, however, that you have a Hamiltonian that is H= lambda * H_1 + H_2, then you can calculate what the energy difference would be between your value of lambda and any other value of lambda. That is what this interpolation business does for you, for all values of the energy encountered in the xvg (or edr) file. It then uses these interpolated energy diferences to calculate the BAR free energy difference (which relies on the average exponent of the energy difference and therefore can't use the pre-averaged dH/dl.
>
It seems nevertheless there is a bug in g_bar. As said we're using
non-evenly spaced lambda's in order to get better sampling close to 0
and 1, to be precise we have 0.0, 0.02, 0.04, 0.07, 0.1, 0.15, 0.2 and
so on. In the histogram xvg I get
@ s0 legend "N(dH/d\xl\f{} | \xl\f{}=0)"
@ s1 legend "N(\xD\f{}H(\xl\f{}=0.05) | \xl\f{}=0)" <== non existant lambda
@ s2 legend "N(\xD\f{}H(\xl\f{}=-0.03) | \xl\f{}=0.02)" <== negative lambda
@ s3 legend "N(dH/d\xl\f{} | \xl\f{}=0.02)"
@ s4 legend "N(\xD\f{}H(\xl\f{}=0.07) | \xl\f{}=0.02)"
@ s5 legend "N(\xD\f{}H(\xl\f{}=-0.01) | \xl\f{}=0.04)" <== negativ lambda
@ s6 legend "N(dH/d\xl\f{} | \xl\f{}=0.04)"
It looks a lot like the program reads a given file at the right lambda
and then adds/subtracts 0.05. The resulting curves in the histogram do
not look good either. I could not find 0.05 in the source code however.
More clues?
>>
>>
>>> Numerical integration is already implemented in g_average if I remember correctly.
>>> You're right that it wouldn't be a bad idea to make the linearized BAR an option, and standard numerical integration the default. g_bar wouldn't be the right name for that tool though :-)
>>>
>>> Sander
>>>
>>> On 14 Nov 2012, at 13:52 , "Shirts, Michael (mrs5pt)" <mrs5pt at eservices.virginia.edu> wrote:
>>>
>>>> Sander, why the requirement that the derivative data be from linear
>>>> coupling? I should be fairly easy to put numerical integration in if you
>>>> have all of the files.
>>>>
>>>> Best,
>>>> ~~~~~~~~~~~~
>>>> Michael Shirts
>>>> Assistant Professor
>>>> Department of Chemical Engineering
>>>> University of Virginia
>>>> michael.shirts at virginia.edu
>>>> (434)-243-1821
>>>>
>>>>
>>>>> From: Sander Pronk <pronk at kth.se>
>>>>> Reply-To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
>>>>> Date: Wed, 14 Nov 2012 11:08:26 +0100
>>>>> To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
>>>>> Subject: Re: [gmx-developers] Error of violate the Second Law of
>>>>> Thermodynamics in Free energy calculation with BAR
>>>>>
>>>>> 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.
>>>>>
>>>>> --
>>>>> 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.
>>>>
>>>> --
>>>> 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.
>>>
>>
>>
>> --
>> David van der Spoel, Ph.D., Professor of Biology
>> Dept. of Cell & Molec. Biol., Uppsala University.
>> Box 596, 75124 Uppsala, Sweden. Phone: +46184714205.
>> spoel at xray.bmc.uu.se http://folding.bmc.uu.se
>> --
>> 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.
>
--
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone: +46184714205.
spoel at xray.bmc.uu.se http://folding.bmc.uu.se
More information about the gromacs.org_gmx-developers
mailing list