[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 15:23:08 CET 2012
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.
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.
More information about the gromacs.org_gmx-developers
mailing list