[gmx-users] (dH/dl) calculation

[David Mobley]

Lei,

> I am also confused about the details of dgdl (or dHdl) calculation in
> GROMACS. For example, if the lambda is fixed at certain value, how did the
> program determine dHdl?

See below.

> My initial understanding was kind of similiar to that of Mauricio that dHdl
> equals the difference in total potential energy between two time frame
> (snapshots).

That's never even approximately true unless you are changing lambda as
a function of time, and even then, it is not really true, as dH/dl
isn't evaluated using that approach, and it also contains an extra
factor of 1/dl relative to the difference in total potential energy,
plus there are dynamics going on...

 However, this seems does not make sense when lambda is fixed.
> Could you give more details about the calculation of the derivative dHdl at
> each snapshot? Thank you.

Try reading the manual, for example, equation 4.116 and 4.120 in the
Gromacs 3.3 manual. If you read that, and my previous e-mail, where I
talked about taking the derivative, and you're still confused, write
back and try to be more specific about the nature of your confusion.
You're talking about evaluating dH/dl using a finite-difference
scheme, which is not how it's done -- it's done analytically, since
H(x,p,l) is known.

Take this example: Suppose f(l)=sin(l), and you want to know df/dl.
You ask, how do we evaluate this while keeping l fixed? Simple: Take
the derivative of sin(l) with respect to l. You could *also* do it
with a finite difference scheme (f(l+dl)-f(l))/dl or some such), but
there's no reason to do that, since the derivative approach is
simpler.

Isn't this clear from the manual? As much as I think the manual is
unclear sometimes, this isn't one of those times...

David

> Best regards,
>
> Lei Zhou
>
>
>
> On 11/6/06, David Mobley <dmobley at gmail.com> wrote:
> > Mauricio,
> >
> > I'm somewhat confused by your question and notation. However, I think
> > the basic answer is something like this: In molecular dynamics, you
> > know the Hamiltonian from which you are sampling; call it H(x,p, l),
> > where x denotes all of the positions, p the momentums, and l lambda.
> > This, of course, is closely linked to the potential energy. Anyway, at
> > any snapshot, you can simply take the derivative dH(x,p,l)/dl, and you
> > have dH/dl at that snapshot. This is usually straightforward since you
> > know the dependence of all of the terms in your Hamiltonian on lambda,
> > so you actually have the functional form for dH/dl as well -- so it
> > just involves taking the appropriate combination of positions,
> > momentums, etc. This is of course all handled internally by the code.
> > <dH/dl>, then, is just the time-average of dH/dl, which can be
> > evaluated every step by the code.
> >
> > I am not sure if that's helpful at all, as I'm not entirely sure what
> > problem you're having. After all, whenever you do TI calculations in
> > GROMACS, the code gives you back dG/dl (or dH/dl, or dA/dl) for every
> > snapshot in an xvg output file. Are you just confused about how the
> > code gets this (I think I just answered that above), or are you trying
> > to figure out how to use it? If you're confused about how to use it,
> > try to ask a question that relates to the specific issue you're
> > confused about.
> >
> > Best wishes,
> > David Mobley
> > UCSF
> >
> >
> > On 11/5/06, Mauricio Sica <msica at unq.edu.ar > wrote:
> > > Dear experts
> > >
> > > I am doing FEP (thermodynamic integration method) simulations.
> > > I have a questions about <dH/dl> calculation in GROMACS.
> > > Take in mind equation 3.77 from the GROMACS 3.3 manual.
> > > There, dA/dl is calculated as
> > >
> > > dA/dl = SS{ (dH/dl) exp()dp dq } / SS{ exp()dp dq = <dA/dl>NVT;l }
> > >
> > > where SS are doble integrals (sorry for the notation).
> > >
> > > My question is: how is (dH/dl) (in the middel-term of the equation)
> > > calculated?
> > > My idea is that the difference V(L=1)-V(L=0) is calculated for every
> time
> > > step (irrespective of the lambda value of the simulation) and <dG/dl> is
> > > the time average of that difference.
> > >
> > > <dG/dl> = < V(L=1)(i)-V(L=0)(i)/1 >
> > >
> > > Is this correct?
> > >
> > >
> > > Thanks
> > >
> > >
> > >
> > > _______________________________________________
> > > gmx-users mailing list    gmx-users at gromacs.org
> > > http://www.gromacs.org/mailman/listinfo/gmx-users
> > > 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/mailing_lists/users.php
> > >
> > _______________________________________________
> > gmx-users mailing list    gmx-users at gromacs.org
> > http://www.gromacs.org/mailman/listinfo/gmx-users
> > 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/mailing_lists/users.php
> >
>
>
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> 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/mailing_lists/users.php
>
>