Subject: Re: CCL:Energy conservation of molecular dynamics in the presence of electrostatic forces

From: Rick Venable <[email-masked]>
Date: Thu, 6 Mar 2003 17:00:38 -0500

I feel obliged to comment that CHARMM most certainly *does* conserve
energy with NVE, and that is put the test frequently as new methods are
developed and added.  Normally, a 1 fs timestep is fine, and one can
even use a 2 fs timestep if bonds to hydrogen are constrained via SHAKE.
Most of the testing and usage is done with condensed phase systems,
mostly fluid.

I should note that systems with a small number of particles (e.g. less
than a few hundred water molecules) behave poorly and may appear not to
conserve energy because of large fluctuations in energy and forces.

There was a minor bug that led to non-conservation of energy for
extended system dynamics (NPT, NVT, etc.) in the initial CHARMM
implementation, but that has been fixed.  I am currently running some
NPT simulations with 2 fs timestep and SHAKE which conserve the total
system Hamiltonian quite nicely.  The system includes numerous
zwitterions with full formal +1 and -1 charges (phosphocholine lipids).

Allen & Tildesley discuss energy conservation in section 3.5, "Checks on
Accuracy"; it's the first topic mentioned.

On Thu, 6 Mar 2003, Charles Xie wrote:
> I have a general question about totoal energy conservation of classical
> molecular dynamics (MD) in the presence of electrostatic forces.
>
> Total energy is perfectly conserved when we do a typical MD with
> only the wan der Waals (Lennard-Jones) potentials, using a time
> step of one femtosecond (e.g. for the Argon gas as a classical example).
>
> Problems, however, arise when molecules are charged. Because
> electrostatic forces are much more long-ranged (1/r dependence) and
> typically stronger, the magnitude of electrostatic forces on atoms can
> be several hundred (or even thousand) times greater than that of van
> der Waals forces. As a result, the time step has to be accordingly
> reduced for two or three orders of magnitude, in order for the total
> energy to conserve. This means a time step of 0.001-0.01 femtosecond
> has to be applied.
>
> (EXPLAIN: Most standard algorithms, such as the Verlet method, the Gear
> predictor-corrector method and the Runge-Cutta method, depend on using
> the term a*dt^2 to compute the solution stepwisely, where a is the
> acceleration,
> dt is the time step. If this term is too large, the numerical error will
> rapidly
> propogate, and the solution will diverge. Compared with pure Lennard-Jones
> simulations, adding the electrostatic forces increase accelerations.
> Therefore,
> dt must be decreased in order to keep a*dt^2 down.)
>
> I can see only two ways to solve this problem, both of which turn out to be
> impractical, as you will see later.
>
> One is to accept the reality, use a smaller time step. But, reducing the time
> step to 0.001 femtosecond will slow down the simulation 1000 times (in
> comparison to pure Lennard-Jones simulations). You will probably never see
> any emerging behavior of the model because losing your temper.
>
> The other way around is to use smaller charges, as most molecular mechanics
> force fields do. For example, instead of using 1 for a free radical's charge,
> one can use 0.01 (and self-explain that this is an effective-field
> approximation).
> While this would maintain energy conservation, you might not see the
> emerging behaviors due to charges either --- because the electrostatic forces
> are so weaker than they ought to be, they might not be able to produce
> results that only strong interactions can exhibit, such as self-organization.
>
> Energy conservation is generally not an concern for most classic molecular
> dynamics simulations that assume a heat bath (for controlling the temperature
> to a desired value, thus remove the numerical errors resulted from large
> a*dt^2).
> However, getting the energetics right is of paramount importance to modeling
> exoergicity/endoergicity of chemical reactions. Unfortunately, I am not aware
> of anyone mentioning this before. All the molecular simulation books I have
> (Allen and Tildesley, Leach, Rapaport and so on) do not even mention it at all.
> My prior experience with CHARMm v27, if I was doing it right, is that it
> doesn't
> seem to do energy conservation, even though the NVE protocol is literally
> specified.
>
> I would greatly appreciate any advice and opinion.
>
> Thank you,
>
> Charles Xie
>
> International Center
> The Concord Consortium
> www.concord.org
>
>
>
> -= This is automatically added to each message by mailing script =-
> [email-masked] -- To Everybody  | [email-masked] -- To Admins
> Ftp: ftp.ccl.net  |  WWW: http://www.ccl.net/chemistry/   | Jan: [email-masked]
>
>
>
>
>
>

=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=
"Don't blame me, I voted for Kang."
                         Homer
-------------------------------------
Rick Venable           29/500
FDA/CBER/OVRR Biophysics Lab
1401 Rockville Pike    HFM-419
Rockville, MD  20852-1448  U.S.A.
(301) 496-1905   [email-masked]
ALT email:  [email-masked]
=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=