[gmx-developers] Error of violate the Second Law of Thermodynamics in Free energy calculation with BAR
Shirts, Michael (mrs5pt)
mrs5pt at eservices.virginia.edu
Wed Nov 14 13:52:54 CET 2012
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.
More information about the gromacs.org_gmx-developers
mailing list