Date: Tue, 24 Sep 2002 07:14:01 -0700
From: "David A. Case"
Subject: Re: expression of V(lambda) in GIBBS

On Tue, Sep 24, 2002, hagop demirdjian wrote:
>
>
> What is the expression of V(lambda) in GIBBS (I could not find it in the
> manual, so I assumed it is linear : V(l) = l*V1 + (1-l)*V0

The mixing depends on the value of IOLEPS (see p. 173 of the manual); I
agree that the documentation could be more specific here. Also, if you
are doing thermodynamic integration in gibbs, the mixing is always linear,
as you write above.

Note also that, in sander, lambda=0 is the "regular" state in LEaP, and
lambda=1 is the perturbed state. The reverse convention is used in gibbs
(see. p. 158, section 7.1).

>
> I simulated a mutation where atoms disapear. I computed DeltaGs in GIBBS and
> SANDER : in GIBBS, dV/dl values are all negative and decrease from -16 to
> -30 kCal/mol/l in SANDER, I used klambda = 4, and got values decreasing
> quickly from 350 to -20 and increasing to 0 at lambda =1.
>
> Both integrations give the same DeltaGs. I assume that differences are
> due to a difference in the expressions of the V(l). Am I right ?

This should be the case, although (see below) one would really need to know
more about how you did the gibbs calculations to say more. I'd be very
interested in seeing your results (and input files) if you feel comfortable
with that. We'd like to generate a tutorial that compares making atoms
disappear in sander vs. gibbs, to understand better what the potential
advantages and pitfalls are in each approach.

>
> And finally, what in your opinion is best suited for TI with vanishing
> atoms ? : SANDER with klambda = 4 or GIBBS with an appropriate value ISDX0 ?
>

This is quite a complex question. Certainly, gibbs has been the "default" way
in Amber for some time, but you need to be careful about setting a number of
parameters, in addition to ISDX0. Care needs to be taken in gibbs about
handling perturbations that change bond lengths. In April, Dave Pearlman
posted the following:

In general, for ANY Gibbs free energy simulation you should turn on SHAKE for
all bonds, set NCORC=1 and set INTPRT = 5. These should be the defaults for
a free energy simulation.

Sander is much simpler, which has both good a bad points. With klambda=4, the
variation in <dV/dl> from one lambda value to another can be pretty steep,
esp. for small lambda. This may make the integral difficult to evaluate, and
we are still working on deciding the best way to handle this. Furthermore,
since the internal contributions are not removed, you have to be sure to
subtract two comparable sides of a thermodynamic cycle (e.g. make the atoms
disappear in both gas and solution, and take the difference.) On the plus
side, TI in sander runs in parallel, uses the latest advances in PME
methodology, and can be run with the latest force fields. But if the ISDX0
"trick" used in gibbs (or something similar) is really the best, then it
should be moved into sander as well; what the best course to take here is
still a matter of discussion among the developers.

Disclaimer: most of the free energy calculations I do are simpler, "charging"
calculations, where one is just changing from a deprotonated to a protonated
form of a side chain (or similar). In a sense, this is what sander was
designed to be "good" at, but it should also give correct results for a much
wider range of perturbations. Comments from people with more experience here
would be very welcome.

..good luck....dac