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

From: Charles Xie <[email-masked]>
Date: Thu, 06 Mar 2003 09:20:27 -0500

Dear CCLers,

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

---

From: "David A. Case" <[email-masked]>
Date: Thu, 6 Mar 2003 09:20:21 -0800

On Thu, Mar 06, 2003, Charles Xie wrote:
> der Waals forces.

I think something is wrong here: the electrostatic forces for reasonable
chemical models are not larger than vdW forces, and all of the significant
contributions to the force come from nearby atoms.  Is it possible you are
confusing energy with force?

> seem to do energy conservation, even though the NVE protocol is literally 

This program should be conserving energy with a 1 fs time step.  Certainly
the nearly parallel CHARMM program does, as do other popular simulation codes.
Of course, this could depend critically on parameters you would specify,
particularly the cutoffs.

..hope this helps...dac

-- 

==================================================================
David A. Case                     |  e-mail:      [email-masked]
Dept. of Molecular Biology, TPC15 |  fax:          +1-858-784-8896
The Scripps Research Institute    |  phone:        +1-858-784-9768
10550 N. Torrey Pines Rd.         |  home page:                   
La Jolla CA 92037  USA            |    http://www.scripps.edu/case
==================================================================

---

From: Charles Xie <[email-masked]>
Date: Thu, 06 Mar 2003 12:47:05 -0500

Thank you for your answer.

>confusing energy with force?

If the electrostatic energy is greater, shouldn't the forces be greater? I
understand that there are cancellations because of positive and negative
charges scattered around. But overall, the forces should (??) be larger.
For example, suppose we play with a Lennard-Jones system. If epsilon
is increased, we will definitely observe an increasing of melting point.
The implication of this is that the interatomic forces are greater, so
the atoms cooperatively resist the raising of temperature.

I am currently dealing with gaseous phase reactions. What I did is very
simple. I started from a Lennard-Jones gas, added charges to them.
And I found it hard to get total energy conservation.

>particularly the cutoffs.

Is there any rule in determining the charges of atoms in CHARMm's force
fields? They look to a novice like me randomly chosen.

>..hope this helps...dac

Definitely. Thanks a lot.

---

From: M Brunsteiner <[email-masked]>
Date: Thu, 06 Mar 2003 17:37:38 +0000

dear Charles,

I do not agree with what you say.
Energy is usually conserved well, provided you

a) account for the long range interactions
   properly with Ewald or any equivalent method
   or at least use a force-shifted or switched 
   coulomb potential (It is well knwon that applying a 
   cutoff to the (particularly Coulomb) forces
   will heat your system up.
b) use a symplectic integration scheme
   such as any version of the Verlet algorithm
   (and NOT Runge Kutta or Gear, ...)

If you do so and still end up with a divergent energy
then there's something wrong with the program
or your input, if you ask me.

If you allow for chemical reactions, however, then
the energy will of course not be conserved (and if so than this 
would be a coincidence) unless you treat the whole system on an 
ab-initio basis, because, for example, the chemical binding energy 
of a covalent bond is of course not properly included in a 
classical MD Hamiltonian. And you seem to talk about classical MD 
in your mail, don't you ? Chemical reactions are not supposed to 
happen in such a system.

best wishes,
mic

Charles Xie wrote:
> -= 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]

-- 
==========================================================================
"Emotions are alien to me.  I'm a scientist."
(Spock, "This Side of Paradise", stardate 3417.3)

--------------------------------------------------------------------------
Michael Brunsteiner
Centre for Theoretical and Computational Chemistry
University College London
mailto:[email-masked]
http://www.ucl.ac.uk/~uccambr/
--------------------------------------------------------------------------

---

From: Konrad Hinsen <[email-masked]>
Date: Thu, 6 Mar 2003 18:11:12 +0100

> energy to conserve. This means a time step of 0.001-0.01 femtosecond

That doesn't seem quite right. A long range does not imply a strong force, it 
just means that the force on any given atom depends on the positions of many 
others. The contributions can (and often do, in practice) cancel out. 
Moreover, large forces do not automatically imply fast motions, other 
parameters (e.g. temperature) enter as well.

You don't say what kind of systems you study. In liquids and biomolecules at 
room termperature, a time step of 1 fs works perfectly well with 
electrostatic interactions. One does see improvement of energy conservation 
when going to 0.5 fs, but there is no need to go to time steps as small as 
those you quote. So either your system is very different, or there is a bug 
in your force calculation code.

> 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).

That depends on the thermostat that is used. Personally I prefer the extended 
systems method, precisely because there is still a conserved energy that can 
be used to monitor the behaviour of the simulation.
-- 
-------------------------------------------------------------------------------
Konrad Hinsen                            | E-Mail: [email-masked]
Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24
Rue Charles Sadron                       | Fax:  +33-2.38.63.15.17
45071 Orleans Cedex 2                    | Deutsch/Esperanto/English/
France                                   | Nederlands/Francais
-------------------------------------------------------------------------------

---

From: Richard Gillilan <[email-masked]>
Date: Thu, 06 Mar 2003 13:18:05 -0500

Charles:  I am a bit confused by your results. The factor that most
influences integration step size is how rapidly the potential changes
with distance, not its magnitude. Electrostatic forces are much smoother
and change more slowly than van der Waal's terms (1/r vs 1/r^10 or 1/r^14).
A step size of 1fs should give you good energy conservation. You did not
say what code you were using. The usual cause of energy conservation 
problems
is the "stiff" nature of C-H and other similar bonds.

Here are a few possible problems:

1. Numerically check the correctness of your derivatives in the force 
calculation. Symmetric
    finite difference is simple and usually works very well for this and 
can help track down
    problems. If you are using someone else's code, then you can skip this.

2. Are you using periodic boundary conditions? If one atom of a neutral 
dipole pair crosses the
    periodic boundary, you get a sudden discontinuity in the energy 
since the charges are now
    separated by a large distance. You need to apply periodic boundaries 
to neutral charge groups
    as a whole. A traditional way of dealing with this is by using 
cutoffs with a smooth ramp
    to zero at some large distance based on distance between centers of 
charge groups.
    Nowdays, folks use particle mesh Ewald or some other summation 
method to handle electrostatics.

3. Does the energy drift or does it just fluctuate a lot around an average?
    Methods like Gear and Runge-Cutta  are not used in MD because,
    though they give very good short-time accuracy, they lack long term 
stability you get with
    a symplectic integration scheme like Verlet. If you come from a 
Physics background
    you may be accustomed to running much more accurate integrations 
with tighter error
    control.  This is usually not appropriate for large-scale MD since 
the exponential divergence
    of trajectories from their initial conditions renders coordinate 
accuracy meaningless in less than
    a picosecond. MD is used mainly for sampling coordinates. 
Nonetheless, your energy should be
    stable and not drift with time. Running equivalent time periods with 
different time steps should
    show improvement in conservation as steps get small.

Good luck.

   Richard Gillilan
   MacCHESS, Cornell

> dt is the time step. If this term is too large, the numerical error 
> will rapidly
> temperature
> to a desired value, thus remove the numerical errors resulted from 
> large a*dt^2).
> 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.
>

---

From: Konrad Hinsen <[email-masked]>
Date: Thu, 6 Mar 2003 20:21:11 +0100

> charges scattered around. But overall, the forces should (??) be larger.

But cancellation is the key point. In condensed matter, the charge density 
averaged over a region containing a few atoms is always very close to zero, 
as otherwise the electrostatic energy would be so large as to make the system 
unstable.

> I am currently dealing with gaseous phase reactions. What I did is very

That may be the cause of your problem, in gas phase, I would expect energies 
and forces to vary much more strongly over time, and charge cancellation to 
happen only on very large length scales.

I wonder if standard MD is the best choice for that kind of system. Some 
adaptive step-size integrator that slows down only in case of near-collisions 
might be better. Or perhaps even a scheme similar to those used for 
hard-sphere systems, where the time step is obtained by estimating the time 
until the next collision.
-- 
-------------------------------------------------------------------------------
Konrad Hinsen                            | E-Mail: [email-masked]
Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24
Rue Charles Sadron                       | Fax:  +33-2.38.63.15.17
45071 Orleans Cedex 2                    | Deutsch/Esperanto/English/
France                                   | Nederlands/Francais
-------------------------------------------------------------------------------