[gmx-developers] Re: gmx-developers Digest, Vol 16, Issue 12

[Michael Shirts]

Hi, all-

> Although by itself correct, there is another issue here.

Oh, there are always other issues.  Hence, I bring it up on the
discussion list before I just change the CVS code :)
 
> Although Velocity Verlet avoid kinetic energy and pressure complications,
> it introduces complications with the constraint algorithms as the velocities
> need to be constrained as well, whereas there are automatically correct
> with leap-frog.

SHAKE already has a version for velocity verlet -- RATTLE.  I've coded
that up already in my local version.  Including the bookkeeping for
both leapfrog and velocity verlet is one of the things that messes up
the bookkeeping.

> I guess there is no conceptual problem and one can always just constrain
> the full-step velocities. But this does involve coding a velocity
> version for all constraint algorithms. (I have already implemented LINCS code for this).

As far as I can tell in the current code, the velocities are simply
constrained by subtracting the positions.  This is essentially
equivalent to position verlet -- the coordinates are just the same as
any other verlet, but the velocities lose some precision,   Can't this
just be done in velocity verlet as well?  That's also what I've done
in my version.  I'm not seeing what the problem is here, though I
admit that these things are -very- subtle, and I could be missing
something.

By the way, there are other, more complicated issues with constraints
when it comes to pressure control and scaling the system, (which would
occur with any integration algorithm) as I understand (Michel, you
could probably explain more!), we should be using the molecular
virial, as opposed to the atomic virial, and scaling the distance
between molecular centers, not simply atomic centers.  I'll read up on
this a little more, but if anybody has looked into this issue more
carefully, let me know :)

Best,
Michael