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

[Michael Shirts]

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))

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.

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))

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.

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.

Cheers,
Michael