According to AMBER manual, the potential energy function is given as
(for the TI approach and icfe = 1)
V(lambda) = (1 -lambda)^k * V_0 + [1 - (1 - lambda)^k] * V_1
and
dV/d(lambda) = k * (1- lambda)^(k-1) * (V_1 - V_0)
where V_0 and V_1 are the initial and final states, respectively, and
dV/d(lambda) is the derivative of the potential energy with respect to
lambda.
For k > 1 & icfe = 1 & lambda = 1;
V(1,k>1) = V_1
dV/d(lambda) = 0
This system is exactly the same system when u dont do any TI calculation
and just prepare the system according the final state and do the MD
simulation.
[if k = 1 and lambda = 1; then dV/d(lambda) = (V_1 - V_0)].
In runmd.f (AMBER 8), the algorithm to calculate the dv/dl is given as
if( onstep ) then
if( klambda == 1 ) then
do i=1,nren
edvdl(i) = edvdl(i) + ener(i) - ecopy(i)
edvdl_r(i) = edvdl_r(i) + ener(i) - ecopy(i)
end do
else
clfac = klambda*(1.d0 - clambda)**(klambda-1)
do i=1,nren
edvdl(i) = edvdl(i) + (ener(i) - ecopy(i))*clfac
edvdl_r(i) = edvdl_r(i) + (ener(i) - ecopy(i))*clfac
end do
end if
end if
So, when k > 1 and lambda = 1, then it is analytically zero. The algoritm
above will calculate it zero, too.
The question is, whether it converges or not (for a system) when lambda ->
1. In DAC's tutorial, when k = 4 and cut = 9, it does not converges to
zero that good. But this does not mean that when lambda = 1, u will have
a value different than zero for dv/dl. It will always be zero.
I and Harry Stern derived a new mixing function, which will make the
mixing potential, V(lambda), zero at lambda = 0 and lambda = 1. This lets
us to use dummy atoms on both states, initial and final states. I have
been doing some simulations on a 4X4 base paired system, and I am using k
= 6 and icfe = 2. The transformation is very smooth; it slowly increases
to a dv/dl value, and then decreases, and smoothly becomes zero when k =
1.
It is part of AMBER 9 but I can send u the modifed files for AMBER 8, and
u can test it if u are interested.
Best,
PS: It might be a good idea to print out the dvdl values even if lambda =
0 or 1.
On Fri, 10 Mar 2006, B. Lachele Foley wrote:
> This is my understanding. Please correct me if I'm wrong.
>
> It isn't analytically zero at lambda=0 for icfe=1 (icfe=2
> wasn't available in Amber 8).
>
> The mixing function goes to zero at lambda=1 for higher values
> of klambda, but the potential being mixed isn't nicely behaved
> -- it diverges at lambda=1. What matters is the combination
> of the two functions.
>
> If you look at David Case's example
> (http://amber.scripps.edu/tutorial/shirts/index.html), a
> favorable combination of cut values and klambda can, indeed,
> converge to zero at lambda=1, but they don't necessarily.
>
> I don't think I -need- the DV/DL to go to zero at lambda=1.
> Merely being finite is probably ok (if someone who's done more
> math on this disagrees, I'd gladly hear their arguments). I
> think I just need the value to be finite and to know what the
> value is.
>
> :-) lachele
>
>
> >I am not sure how this works with icfe=1, but when I set
> icfe=2, for
> >lambda=0 and lambda=1, the dv/dl values are not printed out.
> They are
> >analytically zero. When I run with lambda=0.05 or
> lambda=0.95, I am
> >getting dv/dl values, but they are close to zero. I am doing
> a dummy to
> >hydrogen and hydrogen to dummy transformations in the TI
> calculation. So,
> >I am guessing that the same thing happens with icfe=1 when
> lambda=1,
> >because analytically it is zero.
> >
> >On Fri, 10 Mar 2006, B. Lachele Foley wrote:
> >
> >> I've done van der Waals perturbation in Amber 8 and Amber 9.
> >> My issue (below) happens in both. I want to know if it can
> >> change for Amber 9.
> >>
> >> The context, here, is disappearing atoms in an explicit box of
> >> water by first turning off charges and then by turning off the
> >> van der Waals.
> >>
> >> When I run the charge perturbation, I get DV/DL values printed
> >> in the output file for all values of clambda.
> >>
> >> When I run a vdW perturbation, I get DV/DL values printed in
> >> the output file for all values of clambda except 1. There
> >> should be a value at 1, even if it's zero. When I say not
> >> printed, I mean that "grep DV\/DL blah_1.000.o" yields nada.
> >> Do you know why the values aren't printed? Can they be
> >> printed? I'm using klambda=6, if that's important.
> >>
> >> You'd never notice this with quadrature. For my systems,
> >> though, having endpoints seems to improve accuracy.
> >>
> >> :-) Lachele
> >> B. Lachele Foley, PhD, '02
> >> Assistant Research Scientist
> >> Complex Carbohydrate Research Center, UGA
> >> 706-542-0263
> >> lfoley.uga.edu
> >>
> >>
> >
> >--
> > Ilyas Yildirim
> > ---------------------------------------------------------------
> > - Department of Chemisty - -
> > - University of Rochester - -
> > - Hutchison Hall, # B10 - -
> > - Rochester, NY 14627-0216 - Ph.:(585) 275 67 66 (Office) -
> > - http://www.pas.rochester.edu/~yildirim/ -
> > ---------------------------------------------------------------
> >
> B. Lachele Foley, PhD, '02
> Assistant Research Scientist
> Complex Carbohydrate Research Center, UGA
> 706-542-0263
> lfoley.uga.edu
>
>
--
Ilyas Yildirim
---------------------------------------------------------------
- Department of Chemisty - -
- University of Rochester - -
- Hutchison Hall, # B10 - -
- Rochester, NY 14627-0216 - Ph.:(585) 275 67 66 (Office) -
- http://www.pas.rochester.edu/~yildirim/ -
---------------------------------------------------------------
Received on Wed Apr 05 2006 - 23:49:40 PDT
