Subject: Re: CCL:Help on SHAKE constraints

From: [email-masked] (Hinsen Konrad)
Date: Mon, 17 Oct 94 15:21:09 -0400

I wish to ask about the use of constraints like SHAKE in molecular
   dynamics.  What are the reasons that one would use SHAKE?  What are
   the deficiencies of SHAKE?  When SHAKE (or the like) is used and you
   experience problems: 1)What are the symptoms of "failure"?, 2)Under
   what conditions does this failure occure? and 3)What is the workaround
   for this failure?

When doing MD simulations, you want to make your time step as big as
possible. However, you must make them small enough to sample correctly
the fastest motions in your system, which are typically bond vibrations.
The idea behind constraints is to eliminate the fastest motions by
simply fixing some internal coordinates (like bond lengths). Once
they are fixed, there are no more fast motions, and you can increase
you time steps.

Of course this means that you are simulating a different system. Some
physical quantities will change as a consequence of introducing
constraints; an obvious example are bond length fluctuations, which
become zero. But if the eliminated degrees of freedom are sufficiently
well uncoupled from the remaining ones, you still get accurate
results for the quantities you are really interested in.

Some people dislike the use of constraints, arguing that since
real molecules do not have these constraints, the simulation shouldn't
have them either. To me, this is not at all obvious. Any MD
simulation makes some rather drastic simplifying assumptions;
the most obvious one is the use of classical mechanics instead
of quantum mechanics. Whether or not to use constraints is simply
part of a definite model, just as the choice of this or that
force field. Whether a model is good or bad needs to be tried
out. Unfortunately, there have been next to no systematic studies
on the effect of constraints on the interesting quantities in
real systems. It should be kept in mind that using a larger
time step not only reduces cost, but also roundoff errors, so
the results might iin some respect get better by using constraints.

If you decide to use constraints, you have to get suitable
equations of motion for the constrained system. There are
two well-established methods from classical mechanics:
1) Add the constraints to the unconstrained equations of
   motion using Lagrange multipliers.
2) Introduce new coordinates that automatically satisfy the
   constraints.
In theory both methods are equivalent, but of course not
in numerical applications. If you follow the first approach,
you have to solve the problem that the constraints are
satisfied only for infinite time steps. To get them exactly
satisfied, you have to solve the constraint equations explicitly
for a finite time step. You then get a set of quadratic equations,
which has to be solved iteratively. This is exactly what SHAKE
does. And that also explains how it might fail: There is no
way to prove that the iteration converges, and sometimes it
doesn't. Most implementations simply stop after some number
of iterations. I have been told that the "solution" is simply
to make a few steps without SHAKE and then go on as if nothing
had happened.

The second approach has rarely been tried, probably it looks
more complicated at first glance. However, it is possible
to derive general equations of motion with just a few types
of motion (translation and some kinds of rotation), and it
is possible to let the computer find these coordinates
automatically. The advantage is that you needn't worry
about the exactness of your constraints. Another advantage
is that you can treat a much wider range of constraints
(rigid subgroups, planar molecules, etc.) Besides, the
coordinates themselves may prove useful for analysis of
the simulation results.
On the other hand, there is a price to pay: you get 
velocity-dependent forces (Coriolis and centrifugal
forces), and the implementation is a bit more
complicated.
If you are interested in this approach, have a look at
the paper "Generalized Euler equations for linked rigid
bodies", by G. Kneller and K. Hinsen, in the September
issue of Phys. Rev. E. It contains the derivation of
the general equations of motion. We will publish something
on the implementation and some results soon.

-------------------------------------------------------------------------------
Konrad Hinsen                     | E-Mail: [email-masked]
Departement de Chimie             | Tel.: +1-514-343-6111 ext. 3953
Universite de Montreal            | Fax:  +1-514-343-7586
C.P. 6128, succ. A                | Deutsch/Esperanto/English/Nederlands/
Montreal (QC) H3C 3J7             | Francais (phase experimentale)
-------------------------------------------------------------------------------