[no subject]

[nicholas #*at*# mulliken.mcs.anl.gov (John Nicholas)]

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 #*at*# tcg.anl.gov