From chemistry-request@server.ccl.net Thu Mar 6 13:11:23 2003 Received: from postoffice.mail.cornell.edu ([132.236.56.7]) by server.ccl.net (8.11.6/8.11.0) with ESMTP id h26IBNa23280 for ; Thu, 6 Mar 2003 13:11:23 -0500 Received: from cornell.edu ([128.84.182.122]) by postoffice.mail.cornell.edu (8.9.3/8.9.3) with ESMTP id NAA15418; Thu, 6 Mar 2003 13:11:16 -0500 (EST) Message-ID: <3E6790DD.4040103@cornell.edu> Date: Thu, 06 Mar 2003 13:18:05 -0500 From: Richard Gillilan User-Agent: Mozilla/5.0 (X11; U; Linux i686; en-US; rv:1.0.1) Gecko/20020830 X-Accept-Language: en-us, en MIME-Version: 1.0 To: Charles Xie CC: CHEMISTRY@ccl.net, bob@concord.org Subject: Re: CCL:Energy conservation of molecular dynamics in the presence of electrostatic forces References: <5.1.1.6.2.20030306081753.00b185a8@secure.concord.org> Content-Type: text/plain; charset=us-ascii; format=flowed Content-Transfer-Encoding: 7bit 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 > 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. > > > >