From chemistry-request -x- at -x- server.ccl.net Thu Mar 6 10:02:05 2003 Received: from web.concord.org ([4.19.234.32]) by server.ccl.net (8.11.6/8.11.0) with ESMTP id h26F24a16690 for ; Thu, 6 Mar 2003 10:02:05 -0500 Received: from fc.concord.org (fc.concord.org [4.19.234.60]) by web.concord.org (8.12.5/8.12.5) with ESMTP id h26F20Vw010942 for ; Thu, 6 Mar 2003 10:02:04 -0500 Received: from dhcp221.concord.org (4.19.234.221) by fc.concord.org (FirstClass Mail Server v7.1) with ESMTP (Sender: qxie -x- at -x- concord.org) transient id 20; Thu, 06 Mar 2003 10:02:00 -0500 Message-Id: <5.1.1.6.2.20030306081753.00b185a8(+ at +)secure.concord.org> X-Sender: qxie %-% at %-% secure.concord.org X-Mailer: QUALCOMM Windows Eudora Version 5.1.1 Date: Thu, 06 Mar 2003 09:20:27 -0500 To: CHEMISTRY "-at-" ccl.net From: Charles Xie Subject: Energy conservation of molecular dynamics in the presence of electrostatic forces Cc: bob #at# concord.org Mime-Version: 1.0 Content-Type: text/plain; charset="us-ascii"; format=flowed 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