21 Dec 06 Molecular Dynamics Simulation
 Alder and Wainwright first introduced molecular Dynamics (MD) in the
late 1950’s to study the interaction of hard spheres. Many insights
were made about the interactions of simple liquids. In 1964, a real
potential was used to simulate liquid argon, and in 1974 the first
realistic simulation was performed to study water. Since then
molecular dynamics simulations have been used to study proteins, DNA
complexes, and other biological research. Recently, mixed quantum
mechanical and classical simulations are being used to study
molecular dynamics including experimental procedures, such as X-ray
crystallography and NMR (nuclear magnetic resonance) structure
determination [1].Â
MD is a method for studying classical statistical mechanics of
well-defined systems through numerical solutions of Newton’s
equations of motion. Ideally we would want to use quantum mechanics
and solve the Schroedinger Equation (SE) to determine the mechanics
of atoms and thus molecules. However, the solution to this equation
for complicated potentials and for molecules with many electrons is
very difficult and for many systems it has not been solved as of
yet. Thus approximations must be made, and MD was created. The
first step, the Born-Oppenheimer approximation, is to assume a
solution to the SE that separates the solutions into one part for
the electrons, and one part for the nuclei. Thus for each set of
nuclear positions, we can solve the electron part of the SE, and
find the electronic energy. This is known as the inter-atomic
potential. Combining this with the nuclear-nuclear repulsions from
the other part of the approximation we can find the total potential
energy of the system. We now have an energy that is a function of
atomic positions, and we call this the potential energy surface.
Next we treat the system of nuclei as classical particles moving
along this potential energy surface.  Newton’s equations of motion
can now be used. The last approximation is to assume an analytic
function that describes the potential energy surface [1].
To perform an MD simulation several components of Statistical
Mechanics need to be taken into account. First, a thermostat is
needed to control the temperature of the system. There are several
types of thermostats that can be used, but the one that has been
proven to work the best, is a combination of the Nose and Hoover
thermostats [1]. This couples the system to a heat bath using a
fictional dynamical variable in the equations of motion. Moldy
gives the user three options of thermostats. One is a simple
scaling parameter, which rescales the system every few timesteps, so
that the temperature remains constant. The two others are
thermostats that couple the system to a heat bath as described
above. These are the Gaussian thermostat and the Nose-Hoover
thermostat. Secondly, an integrator is needed to numerically solve
Newton’s Equations. The typical method is to Taylor-expand the
equations of motion with some alterations. Moldy uses a modified
version of the Beeman algorithm to integrate [3]. Thirdly, a
potential energy function is needed to calculate the forces acting
on the system. This is the function that would come from the
Born-Oppenheimer approximation. For large systems of atoms or
molecules, this is usually divided into an electrostatic and an
inter-molecular potential. For the inter-molecular one, there are
several to choose from depending on what type of system one is
simulating. A typical function is a combination of the nuclear
repulsion and valence electron attractions. The electrostatic part
is handled using the Ewald Sum. This sum modifies the basic concept
that the system is of charged point ions mutually interacting via
the Coulomb potential. The first modification is that each ion is
neutralized at a long range by the superposition of a spherical
gaussian cloud of opposite charge centered on each ion. This
combination of point ions and gaussian charges make up the real-part
term of the Ewald Sum. The second modification is to superimpose a
second gaussian cloud with the same charge as the original point
ions, thus nullifying the effects of the first cloud (i.e. we are
left with the original system). Poisson’s equation is used to
calculate this effect in Fourier or reciprocal space. Additional
terms such as a self-energy correction due to the fictitious
presence of the gaussian clouds are added as well [2]. With these
potential functions, one can compute the forces acting on the
system, and then numerically solve Newton’s Equations of Motion.
[1]Â http://www.ch.embnet.org/MD_tutorial/
[2]Â http://www.dl.ac.uk/TCS/Software/DL_POLY/USRMAN/node64.html
[3] Refson, Keith Moldy User’s Manual, Revision 2.26.2.6, May 24,
2001
Leave a Comment
