[no subject]
Here is my promised summary of responses to my question about SHAKE.
Another question was raised for implementors of molecular dynamics:
does your program have multiple time steps? If not, do you plan to
implement them? I'll summarize any answers here.
Bill Ross
UCSF
From: ross(-(at)-)zeno.mmwb.ucsf.edu (Bill Ross)
To: chemistry(-(at)-)ccl.net
Subject: SHAKE failure
I'm collecting interesting stories of SHAKE failure in molecular
dynamics runs - cases that were never figured out as well as ones
that were. References to anything written on this subject would
be welcome too. I'll summarize to the reflector.
Bill Ross
UCSF
[The inspiration for this question was occasional SHAKE failure
in Amber. Dave Pearlman diagnosed this as stemming from periodic
boundary conditions (constant pressure) where ions are treated
as part of the solute, all solute-solute interactions are included
(no cutoff applied) and so the solute is not imaged with itself:
when an ion crosses the edge of the box it is translated to the
other side, and if another ion is close by, they suddenly "see"
each other and the resulting spike in the virial causes an
expansion of the box to equalize the "pressure," which puts some
bond lengths beyond the ability of SHAKE to recover. In my
opinion, the best way to get around this is to either run with
constant volume or treat the ions as part of the solvent so they
are imaged. Neither way is completely satisfactory, since in
constant volume ions will still suddenly appear in proximity
and when ions are treated as part of the solvent cutoffs are
applied and long-range electrostatics are lost. I have always
run with enough water so that it hasn't happened to me, but
this is partially a matter of luck.]
>From STOUTEN(-(at)-)embl-heidelberg.de Wed Nov 13 01:44:46 1991
Subject: Re: SHAKE failure
To: ross(-(at)-)zeno.mmwb.ucsf.edu
Hi Bill,
I don't have any interesting stories. I wonder, however, why you do your
investigation since SHAKE is being phased out and will be replaced by
multiple time step algorithms. Even Wilfred takes this stand.
Cheers, Pieter Stouten
>From STOUTEN(-(at)-)embl-heidelberg.de Thu Nov 14 16:26:56 1991
Subject: Re: SHAKE failure
To: ross(-(at)-)zeno.mmwb.ucsf.edu
Hi Bill,
On Wed, 13 Nov 91 09:11:32 PST you wrote:
>why even ask? - sociology/folklore of science.
>
That sounds like a good reason. Still, I don't have exciting stories. When
being far from equilibrium I had often problems or when applying heavy
torsion constraints (this seems unrelated but is not). Then I just did not
use shake.
>how soon do you expect multiple time steps to take over?
>
Hard to say. I am a bit away from the field now. I know that already 3 years
ago people were talking about implementing it. As for me, I basically use
GROMOS and considering how busy Wilfred c.s. are I don't know when the first
official release after GROMOS 87 will see the light.
Cheers, Pieter Stouten.
#### # # ### # European Molecular Biology Laboratory
# ## ## # # # Biocomputing Programme
### # # # ### # Meyerhofstrasse 1, D-6900 Heidelberg, Germany
# # # # # # e-mail: stouten(-(at)-)embl-heidelberg.de
#### # # ### #### phone: +49-6221-387 472, fax: 387 517
>From balbes(-(at)-)osiris.rti.org Wed Nov 13 05:40:27 1991
To: ross(-(at)-)zeno.mmwb.ucsf.edu (Bill Ross)
Subject: Re: SHAKE failure
Okay, this won't be much help, but...
I had shake fail after about 14 ps. I had Tom Darden look at it, and he
said that because I was saving the steps in binary form, any restart
would fail exactly the same way. He said (I think) that he saves
things in ascii, so that roundoff on restarting will get around any
failures of this type. I have since been using tom's fast amber 3a,
and haven't had any more problems of this type. Of course the files
have been long since purged.
This was quite awhile ago, so details are fuzzy.
Lisa
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% standard disclaimer %%%%
Lisa M. Balbes, Ph.D. phone: 919-541-6563
Research Triangle Institute, PO Box 12194 vmail: 919-541-6767, xt 6563
Research Triangle Park, NC 27709-2194 email:
balbes(-(at)-)osiris.rti.org
- This came directly from a computer and should not be doubted or disbelieved.-
>From harris(-(at)-)athena.mit.edu Wed Nov 13 05:47:02 1991
To: ross(-(at)-)zeno.mmwb.ucsf.edu (Bill Ross)
Subject: Re: SHAKE failure
This example concerns SHAKE's cousin, RATTLE applied to simulation of 50
hexane molecules at liquid density. It is more an example of the effects
of the timestep than a complete failure. Using a timestep of .007 ps, if one
compares the pressure and temperature computed with the molecular virial
with the similar quantities computed using the atomic virial discrepencies
are obversed. In a 1.120 nanosecond (160 000 time steps), the average
molecular temperature is 303.5 compared to an atomic temperature of 300.5
(there is a thermostat attached to maintain the average temperature at 300K by
a weak coupling to an external bath (i.e Berendsen et al. JCP 81, p3684 1984).
The pressure from the molecular virial is also about 20 atm higher than
that computed from the atomic virial. Cutting the timestep in half reduces
the difference between the two methods of computing the temperature and
pressure to 1 degree and 5 atm. It appears that the atomic versions of the
temperature and pressure are more accurate, but the statistics are too poor
to bury the question.
Jonathan G. Harris, H. P. Meissner Assistant Professor,
Department of Chemical Engineering, MIT Rm 66-450
25 Ames Street, Cambridge, MA 02139
harris(-(at)-)athena.mit.edu (617)253-5273 Fax 253-9695
-------
From: nobody(-(at)-)kodak.com
To: "amber(-(at)-)cgl.ucsf.edu"(-(at)-)kodak.com
Subject: RE: SHAKE nightmares.
>From: NAME: Adi M. Treasurywala
FUNC: Biophys. & Compu. Chem.
TEL: (518)445-7042
<TREASURYWALA(-(at)-)A1(-(at)-)DSRGVJ>
To: NAME: Edward P. Jaeger <JAEGEREP(-(at)-)A1(-(at)-)DSRGVJ>,
"amber(-(at)-)cgl.ucsf.edu"(-(at)-)KODAKR(-(at)-)MRGATE(-(at)-)WPC
Bill,
We have completed a fairly detailed study in AMBER3.0A of this
problem. We ran into it because we were trying to do REAL
classical dynamics (ie not constant temperature!!!). What happened
was that in the initial runs on the molecule that we were
interested in we found that the total energy simply did not
stabilize but decreased steadily over the run at a fairly
precipitous rate. There seemed to be no way to stabilize it.
We therefore took a small cyclic peptide WFGLMQ and essentially
banged the heck out of it by trying many different conditions to
narrow down the reason for this "bug". One of us at least
was
convinced that it must be something that WE were doing wrong
rather than that AMBER could have a bug in it that was so glaring!
Anyway, we discovered a LOT of other interesting things along
the way such as objective criteria (vs semi religeously based
recipies) for knowing when your thermalization was done, and some
real reasons to start and finish the collection of MD data
(incidentaly it is rather system size dependant...). After much
trial and error we found (thanks are due here to George Seibel who
REALLY put us on to this) that it was indeed the SHAKE option that
had caused us all this grief.
If we turned SHAKE OFF on an otherwise identical run that had
failed with SHAKE ON, it went just fine!! Everything stabilized
and we were able to run real classical dynamics (this is sort of
an emotional issue for us here!!).
We are in the process of preparing a manuscript covering this
work. Let me just ask the gurus on this net... Would that be a
good thing to do? Would it help anyone? Would it insult or offend
anyone? Incidentaly I learned that the real problem with SHAKE was
the integrator (leapfrog just can't handle classical MD with SHAKE
turned on I was told by someone.).
Hope to hear the other stories soon and to hear any discussions
about classical vs constant temp runs!
Adi T & Ed Jaeger.