[gmx-developers] Re: gmx-developers digest, Vol 1 #313 - 1 msg

[Berk Hess]

> 
> 
> Hi, all-
> 
> > I just derived that v(t+0.5*dt) in the constraint direction at t+0.5*dt is
> > zero. I don't know of any references, but I would say that also with
> > constraints the algorithms should be identical (unless you constrain the
> > velocities incorrectly).
> >
> > Berk Hess
> 
> OK, this kind of makes sense.  We start from the assumption that the positions
> are constrainted at all time points, then as far as I can tell from the code
> (in 3.2, update.c, line 789 - is this the wrong place?) the velocities are
> computed by taking the position differences from the constrained differences
> 
> v(t+0.5*dt) = (r(t+dt) - r(t))/dt
> 
> IF we can define r(t + 0.5*t) (which I'm not certain we can within the
> discrete algorithm and preserve symplectic-ness), we would define it within
> order dt (or perhaps even dt^2, I'm not sure) r(t + 0.5*dt) as:
> 
> r(t+0.5*dt) =~ 0.5*(r(t+dt) + r(t))

You can define anything you like.
The important point in the accuracy of the sampled coordinates r.
These are indentical for velocity verlet and leap-frog if you
choose the starting velocities appropriately.
Velocities can be determined afterwards to different accuracies using
more or less r points.
Note that although the velocity verlet integration is time reversible,
the velocity is not reversed properly: v(t) goes to v(t-1).
This means that if you change v(t) using quantities at time t, the
algorithm is no longer time reversible.

> 
> so v_12(t+0.5*dt) . r_12(t+0.5*dt) = 0.5/dt * (r_12(t+dt) - r_12(t)) . (r_12(t+dt) + r_12(t))
> 				   = 0.5/dt * |r_12(t+dt)|^2 - |r_12(t)|^2
>                                    = 0
> 
> However, this is only correct to a certain order (as far as I can tell)
> whereas with rattle, it's true that v(t).r(t) = 0 to arbitrary precision.

This is to any order, unless we want to discuss the accuracy of v, which
is neither fully accurate in leap-frog, nor in velocity verlet.

> 
> Doing the same thing at time t,
> 
> v(t-0.5*dt) = (r(t) - r(t-dt))/dt
> v(t+0.5*dt) = (r(t+dt) - r(t))/dt
> 
> so:
> 
> v(t) = 0.5 * (v(t+0.5*dt) + v(t-0.5*dt))
>      = 0.5/dt * (r(t+dt) - r(t-dt))
> 

The more appropriate expression for v(t) is:
  v(t) = v(t-0.5*dt) + 0.5*f(r(t))*dt = v(t+0.5*dt) - 0.5*f(r(t))*dt
Note that this v(t) does not obey the constraints,
but this can be corrected using for instance lincs, without affecting
the integration, at least analytically.

> and:
> 
> v_12(t) . r_12(t) = 0.5/dt * (r_12(t+dt) - r_12(t-dt)) . r_12(t)
> 
> Which I don't think necessarily equals zero . . . at least not any way I can
> see.

Correct.

> 
> So, you've convinced me that it's not that bad to use leapfrog with
> constraints but not that it's equivalent.  I still don't see any advantage
> over velocity with leapfrog, and leapfrog still has the problems with
> extendend system equations of state that affect both the pressure and velocity
> (i.e.  Parrinello-Rahman, but not Nose-Hoover).
> 
> Really, the proof is in the pudding.  The question is, does it conserve energy
> (in the case of NVE) or enthalpy (in the case of NPH) with no drift _and_ with
> variance proportional to (dt)^2.

In NVE leap frog is identical (except for the velocity definition)
to velocity verlet and both conserve energy.
In NPT (I don't know what you mean exactly with NPH) things depend very much
on the implementation. If you just use the v(t) of velocity verlet the
algorithm is not time reversible any more, but it might still conserve
kinetic energy. The current implementation with leap frog is also not
time reversible and might not have good conservation properties.
Another problem might be that the pressure from the previous step is used.

If you want to do things properly I would suggest to perform
leap frog integration with half steps for the velocity. Thus you
have the velocity at t-dt/2, t and t+dt/2. Such an approach is also
used for third-order stochastic dynamics algorithm with constraint see:
Van Gunsteren and Berendsen, Mol Sim 1, 173-185 (1988) 

I found an appendix of a thesis describing a similar approach for
Nose-Hoover with Parrinello-Rahman:
http://www.wag.caltech.edu/home-pages/gao/thesis/appendix6.pdf

I quickly read through it and it looks reasonable,
but unfortunately no order is given for the approximations
and nothing is said about the properties of the algorithm.

Berk.