[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
Wed Nov 14 17:09:23 CET 2012
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?
> 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
More information about the gromacs.org_gmx-developers
mailing list