[gmx-developers] Re: gmx-developers digest, Vol 1 #313 - 1 msg
Michael Shirts
mrshirts at stanford.edu
Tue Apr 6 01:32:11 CEST 2004
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
More information about the gromacs.org_gmx-developers
mailing list