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



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:
 >
 > 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
 >
 > -= This is automatically added to each message by mailing script =-
 > CHEMISTRY -AatT- ccl.net -- To Everybody  | CHEMISTRY-REQUEST -AatT-
 ccl.net -- To Admins
 > Ftp: ftp.ccl.net  |  WWW: http://www.ccl.net/chemistry/   | Jan: jkl -AatT- ccl.net
 --
 ==========================================================================
 "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:m.brunsteiner -AatT- ucl.ac.uk
 http://www.ucl.ac.uk/~uccambr/
 --------------------------------------------------------------------------