Subject: SUMMARY:predictor-corrector/quaternion for SPC

From: "Tahmid I. Mizan" <[email-masked]>
Date: Tue, 23 Aug 1994 11:39:46 -0400

Thanks to everybody who responded to my question regarding the predictor- 
corrector/quaternion method as appled to SPC water. The quaternion method, 
of course, unlike SHAKE, solves the rigid body equations of motion for the 
molecule as a whole. Analytic stability analysis of MD equations is not
possible. Some sort of approximate analysis may be possible by linearizing
the differential equations or by Lyapunov methods. I did not attempt this, 
but rather took a trial-and-error approach. It turns out that for 2 fs or
larger time step, iterating on the evaluation-correction step solves the 
problem. Of course, this means an extra force evaluation per time-step which
is undesirable. However, for 3 fs timestep a P(EC)**2 (two force evaluations
per time step) appears to be stable so that there is some computational
advantage over a 1 fs time-step ( 0.67 force evaluations per fs as opposed
to 1 force evaluation per fs). The 3 fs runs give the correct values of 
all static properties ( I have not calculated autocorrelation functions).
I should mention that for a 4 fs time step, the situation deteriorates and
a P(EC)**3 is required for stability.

T. I. Mizan
University of Michigan

*******************************************************************************
ORIGINAL QUESTION

I am using the predictor-corrector/quaternion subroutines from Allen & Tildesley
to simulate SPC water. When I use a 1 femtosecond timestep the program runs    
properly and gives the correct results for energy, pressure g(r) etc. However, 
when I use a 2 fs timestep the program 'blows up'. Apparently the OH bond lengths
become larger and larger and it just blows up. So it appears to be some kind
of numerical instability. 

Has anybody else had a similar experience ? Or is it just some error in my  
implementation.

Is there a fix ?

  -- Tahmid Mizan
*******************************************************************************
>From [email-masked] Tue Aug 16 21:21:04 1994

Perhaps I can answer your questions about your water dynamics.  I also
do simulations of SPC water, although not with predictor-corrector,
or quaternions.  So I can't say whether there's a bug in Allen & Tildesley's
code.  But timesteps of 1-2 femtoseconds are about as large as you can
usually use on rigid water.  Some people use as high as 2 fs, I can't go
any higher than 1 fs, but it depends a little on your algorithm.  So this
makes it sound like it's not a bug but just an instability at large timesteps.
But the fact that your OH bond(s) keep growing makes it sound like there
might be some other problem.  That isn't the sort of blowup I'd expect
>from timesteps, although I haven't used quaternions, so I could be mistaken.

Hope this helped.

-Steve Stuart
[email-masked]
*******************************************************************************
>From [email-masked] Tue Aug 16 21:21:14 1994

Are you using SHAKE on your water bonds? It seems to be necessary
with models that use a larger O vdw to enclose the H's, because
the margin of coverage for the H charges is too small for bond 
stretching to be acceptable w/ longer time steps.

Bill Ross
*******************************************************************************
>From [email-masked] Wed Aug 17 00:01:02 1994

Tahmid,

Are you constraining the geometry of the water molecule (bond and angle)
with SHAKE?  With SPC, your hydrogens have no lennard-jones terms and
but they do have charge.  If you don't constrain these hydrogens they
can easily cause your calc. to blow up.

Mark

Mark A. Thompson                    
Sr. Research Scientist              email:  [email-masked]
Molecular Science Research Center   FAX  :  509-375-6631
Pacific Northwest Laboratory        voice:  509-375-6734
PO Box 999, Mail Stop K1-90
Richland, WA.  99352
*******************************************************************************
>From [email-masked] Wed Aug 17 01:50:59 1994

Hello

I am responding to your question about stability of A&T algorithm for
water. Although I have not worked with this, my officemate has, but
is in taiwan until the end of the month. If you are still interested
in comunicating with him, his e-mail is [email-masked]

(his name is Felix Li)

Hope this little bit helps...

Keith Ball
U of C physics grad
*******************************************************************************
>From [email-masked] Wed Aug 17 03:41:07 1994

Hello Tahmid,

In the case of the Verlet algorithm and nonbonded atoms the blowing up
is a typical behaviour when the time step is increased. See for
example the User Guide, Part 1 for Discover 2.9/3.1, January 1993,
Biosym, pp.  1/2-58 - 1/2-61: As the two nonbonded atoms approach each
other and the time step is long they may "deeply interpenetrate each
other's steeply repulsive wall between steps. The resulting force is
now so large that the atoms fly off at a speed about 1 A/fs (or 10^5
m/s). An equivalent temperature would be in the hundreds of thousands
of degrees".

If a Morse potential is used for the bond-streching term in the
forcefield also the bonded atoms probably behave by this way, because
this potential allows the bonded atoms to drift unrealistically far
apart. If you are using harmonic potentials then I don't know what
happens.

There is some discussion about possible errors in a MD program in the
Section 3.5 "Check on accuracy" (pp. 98-101) in the Allen &
Tildeslay's book "Computer Simulations of Liquids"

Regards,

   Raimo

------------------------------------------------------------------------
Raimo Uusvuori                       |  Phone:  +358-0-457 3210
Scientific Software Specialist       |  FAX:    +358-0-457 2302
Chemistry                            |
                                     |
Center for Scientific Computing      |  Street address:
P.O. Box 405                         |       Tietotie 6, FIN-02150 ESPOO
FIN-02101 ESPOO                      |  E-mail:    [email-masked]
FINLAND                              |             [email-masked]
------------------------------------------------------------------------
*******************************************************************************
>From [email-masked] Wed Aug 17 09:16:40 1994

This is in reply to your question about 1 vs 2 fs timesteps. I assume you are
talking about the integration timestep.
I think 2 fs is getting too long for a timestep in water since the vibrational
period of the highest-frequency normal mode is about 9 fs (3756 cm-1) If you
increase the timestep further, does the program still blow up? does it work if
you decrease the timestep below 1 fs?

Hope this helps,
Alex Matro
*******************************************************************************
>From [email-masked] Wed Aug 17 09:24:21 1994

With regards to your recent posting to the CCL:

I think that the predictor-corrector/quaternion subroutines you are
using are guaranteed to exhibit the kind of numerical instability you
have observed when the timestep is "too large." How large "too large" is
will depend entirely on the system you are modeling (a timestep for water
might be overkill for handling, say, liquid argon). 

This is based on limited experience (i.e., I haven't implemented the
subroutines in Allen/Tildesley you are using). I will be interested to
hear if someone tells you differently.

Best wishes!

robert topper

 Prof. Robert Q. Topper               internet:   [email-masked]
 Department of Chemistry              phone:      (212) 353-4378
 The Cooper Union                     FAX:        (212) 353-4341 
 Cooper Square                        subway:     take the N/R to 8th/NYU 
 New York, NY 10003 USA                           or the 6 to Astor Place

 The Cooper Union for the Advancement of Science and Art, established 
 by Peter Cooper in 1859, is a private institution of higher learning 
 where all students receive full-tuition scholarships.                 
*******************************************************************************
>From [email-masked] Wed Aug 17 12:11:45 1994

Tahmid,
I would like to have some info on the SPC water model.  Will you please
provide some References?
Thanks.

By the way, I do not know if my answer is too simplistic, but it appears
to me that P/C quaternion is unstable at 2 fs , so then your total energies
and other vbles go crazy (no Energy conservation, I'll suspect).  Play with it
>from 0.5 fs to 1.5 fs and see if you get conservation of E.
(for example at if you run 150 steps at 0.5 fs it should match the energy
at the 50 th step of the 1.5 fs, if it does NOT, then 1.5 fs is NO GOOD, andso
forth).  Let me know how it goes

regards,
khalid
*******************************************************************************
>From [email-masked] Mon Aug 22 18:19:16 1994

Dear Tahmid,

Do you know what is the maximum allowable time step
for the integrator you are using ?

Note that most probably the method you are using is
only "conditionally stable" i.e. stable for small
time step and unstable for time step grater then
certain treshold.

For more detail analysis you may refer any
handbook devoted to MD and/or numerical
solutions for differential equations
e.g. J.M. Haile "Molecular Dynamics Simulation: Elementary Methods"
J. Wiley & Sons, Inc.
and most probably Allen & Tildesley as well.

Of course this does not exclude possibility
of the bug in your implementation - unfortunately.

Sincerely,

Ryszard Czerminski, MSI
*******************************************************************************

-------   Regard
          Tahmid Mizan