From chemistry-request@ccl.net Tue Nov 19 02:21:15 1991
Date: Tue, 19 Nov 91 01:21:14 CST
From: nicholas@mulliken.mcs.anl.gov (John Nicholas)
To: chemistry@ccl.net
Status: R

Dear Dr. Bill Ross and other SHAKE users: 

I have some experience
with the SHAKE algorithm in MD that pertains to the recent discussion.
I have used SHAKE in my own MD code for some time and regularly do 
runs of hundreds of thousands of time steps with excellent energy 
conservation, etc. These are my opinions about some of the problems 
mentioned:

1. Dr. Adi M. Treasurywala writes that SHAKE cannot be done with the
leapfrog algorithm. I have used SHAKE with van Gunsteren's
formulation of the leapfrog algorithm (J. Comp. Chem.) with no 
problems. The implimentation of SHAKE within AMBER should be the same
as the dynamics code is essentially the same as van Gunsteren's (GROMOS).
I'm not really an AMBER user so this might vary with different
versions of the program.
The problem with energy conservation is related to the tolerance
used in SHAKE. The value suggested by van Gunsteren (1.0e-4) is
too small to prevent persistent energy drift. Changing to a value
of 1.0e-6 should cure these problems. I believe CHARMM uses 1.0e-9
which does not result in an improvement in energy conservation in my
work over 1.0e-6. The price you pay for a smaller tolerance is an
increased number of iterations for SHAKE to converge. Even with
very small tolerances SHAKE will usually converge within about
20 iterations. 

2. Dr. Lisa Balbes writes about problems restarting with binary files. 
If a true
restart is desired with a continuous trajectory one certainly would
not want to disturb the run by applying an arbitrary round-off of the
positions and/or velocities. I have had no trouble with
restarts using double precision binary files and there is nothing
inherant in the use of SHAKE that should prevent this.

3. Several researchers mention mutiple time step methods. I have used
the method of Otto Teleman (J. Comp. Chem.) and Allen and Tildlesly
describe another method in their book "Computer Simulations of
Liquids". Both methods work fine. However, SHAKE is much easier
to code if all you want to do is constrain bond lengths, which 
gives the most increase in time step (best bang for the buck).

4. The mention of problems with constant pressure simulations can
be alleviated by coding the rescaling of coordinates to maintian
bond lengths. For example, moving the centers of mass of molecules
rather than actually stretching or contracting them to respond
to pressure. I haven't done much with this so I am not clear if
a particular method gives better or worse results. Some of the other
problems with periodic boundary conditions that Dr. Dave Pearlman
describes are just code bugs that could be easily fixed.


John Nicholas
nicholas@tcg.anl.gov