From chemistry-request _-at-_)server.ccl.net Thu Mar 6 12:49:42 2003 Received: from mail-d.bcc.ac.uk ([144.82.100.24]) by server.ccl.net (8.11.6/8.11.0) with ESMTP id h26Hnga22453 for ; Thu, 6 Mar 2003 12:49:42 -0500 Received: from pop-b.ucl.ac.uk by mail-d.bcc.ac.uk with SMTP (Mailer) with ESMTP; Thu, 6 Mar 2003 17:49:39 +0000 Received: from pop-d.ucl.ac.uk (pop-d.ucl.ac.uk [144.82.100.34]) by pop-b.ucl.ac.uk (8.11.6+Sun/8.9.3) with ESMTP id h26HnaM09646; Thu, 6 Mar 2003 17:49:36 GMT Sender: mb \\at// ucl.ac.uk Message-ID: <3E678762.BFE86976 #at# ucl.ac.uk> Date: Thu, 06 Mar 2003 17:37:38 +0000 From: M Brunsteiner Organization: UCL X-Mailer: Mozilla 4.76C-SGI [en] (X11; I; IRIX64 6.5 IP30) X-Accept-Language: en MIME-Version: 1.0 To: Charles Xie CC: CHEMISTRY #at# ccl.net, bob #at# concord.org Subject: Re: CCL:Energy conservation of molecular dynamics in the presence ofelectrostatic forces References: <5.1.1.6.2.20030306081753.00b185a8-!at!-secure.concord.org> Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit X-UCL-MailScanner: Found to be clean 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 %-% at %-% ccl.net -- To Everybody | CHEMISTRY-REQUEST %-% at %-% ccl.net -- To Admins > Ftp: ftp.ccl.net | WWW: http://www.ccl.net/chemistry/ | Jan: jkl#* at *#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 /at\ucl.ac.uk http://www.ucl.ac.uk/~uccambr/ --------------------------------------------------------------------------