[gmx-developers] constraints and BAR free energy

[Shirts, Michael (mrs5pt)]

> Well, the integration error is the whole thing.
> BAR is significantly better than free energy integration
> and it is unbiased.

It depends.  If you are introducing atomic sites in, then BAR is much better
than TI; you can get the same result with significantly (1/2 to 1/3) as many
states.  If you are changing a harmonic bond or changing a charge linearly,
the advantage is much smaller.

And BAR is not technically unbiased -- it is asymptotically unbiased, which
means that is only unbiased as the number of samples goes to infinity.
However, in most molecular cases, the bias is swamped by the uncertainty --
it goes to zero pretty fast.

I would suggest performing the transformation in two steps; first, do
everything but the constrained lengths with BAR, then change the constraints
with TI.

Please keep in touch with me on how this goes; I'm finishing up a major
expansion of the free energy code, and I'd like to make things like this as
easy to do as possible,

If there is some uncertainty about the current constraint code, I'd perform
a test case of solvation of a diatomic dipole with different constraint
lengths to see if you can get the thermodynamic cycle to converge.

Best,
~~~~~~~~~~~~
Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434)-243-1821


> From: "hess at sbc.su.se" <hess at sbc.su.se>
> Reply-To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
> Date: Wed, 20 Oct 2010 08:33:39 -0400
> To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
> Subject: Re: [gmx-developers] constraints and BAR free energy
> 
> Well, the integration error is the whole thing.
> BAR is significantly better than free energy integration
> and it is unbiased.
> 
> We are thinking about doing BAR properly with constraints,
> but that's not trivial.
> 
> Berk
> 
>> What's the approximation for case 1), apart of course from the integration
>> error? Isn't it just an analytical derivative dV/dl?
>> Agreed that in most normal cases the error for case 2) should be small but
>> my
>> superstition against perturbing constraints does seem somewhat justified,
>> especially since in most cases with 2fs time steps it can be avoided by
>> not
>> mutating hydrogen to non-hydrogen...
>> thanks!
>> 
>> Floris
>> 
>> 
>> 
>> 
>> ----- Original Message ----
>> From: "hess at sbc.su.se" <hess at sbc.su.se>
>> To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
>> Sent: Wed, 20 October, 2010 11:03:51
>> Subject: Re: [gmx-developers] constraints and BAR free energy
>> 
>> Hi,
>> 
>> Strictly speaking BAR is impossible with changing constraint lengths.
>> Since the constrained distances are different for different lambda's,
>> the ensembles never overlap.
>> 
>> So there are two options:
>> 1) I forget about BAR and do free energy integration
>> 2) I use Delta_lambda*dV/dlambda for the constraint contribution
>> 
>> 1) Has the same "approximation" for the contraints as 2) and additionally
>> has larger and systematic errors. So I would go for 2).
>> 
>> In most situations constraint lengths change very little between
>> adjacent lambda points, so 2) should be a very good approximation.
>> 
>> Berk
>> 
>>> Hello,
>>> 
>>> I have a question about the perturbation of constraints. My
>>> understanding
>>> is
>>> that dv/dlambda is easily calculated for constraints, so that constraint
>>> perturbation is fine for thermodynamic integration - but that
>>> calculating
>>> the
>>> deltaV to another lambda value as required for e.g. BAR doesn't make
>>> sense
>>> or
>>> isn't practical. Here's the relevant comment in force.c:
>>> 
>>> /* Notes on the foreign lambda free energy difference evaluation:
>>>  * Adding the potential and ekin terms that depend linearly on lambda
>>>  * as delta lam * dvdl to the energy differences is exact.
>>>  * For the constraint dvdl this is not exact, but we have no other
>>> option.
>>>  */
>>> 
>>> So far I've always avoided perturbing constraints but I now need to
>>> recalculate
>>> some results based on topologies with a ring-bound hydrogen mutating to
>>> a
>>> methyl
>>> group. Is there an exact way to do this in gromacs using BAR?
>>> 
>>> More generally I'd be concerned about the accuracy of the linear
>>> extrapolation
>>> of constraint dv/dl - depending on how you set out your transformation
>>> pathway I
>>> would imagine you could end up skipping over a lot of detail in exactly
>>> this
>>> kind of situation where a methyl group is 'pushing' out from a hydrogen.
>>> 
>>> Best wishes,
>>> 
>>> Floris
>>> 
>>> 
>>> 
>>> --
>>> 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.
>> 
> 
> -- 
> 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.