by: Yuekai Sun and Professor Jun Lou
Abstract
Molecular dynamics (MD) is a computer simulation technique that studies particle models by stimulating the time evolution of a system of interacting particles. MD is commonly employed in materials science, nanotechnology and biochemistry to study various processes at the atomic scale. The techniques of MD have also been adapted for use in fluid dynamics and astrophysics. In this article, we introduce the basic science and computational methods behind MD and explore some of its current applications.
Given a system consisting of n particles subjected to interactions described by a given force field, can the trajectory of every particle be predicted? 18th century French astronomer and mathematician Pierre-Simon Laplace said of solving the above problem, called the n-body problem, “Given for one instant an intelligence which could comprehend all the force by which nature is animated and the respective situation of the beings who compose it—an intelligence sufficiently vast to submit these data to analysis—it would embrace in the same formula the movements of the greatest bodies of the universe and those of the lightest atom; for it, nothing would be uncertain and the future, as the past, would be present to its eyes.”
Molecular dynamics (MD) solves the n-body problem by numerically integrating the equations of motion for every particle. There are two common forms of MD, ab-initio MD and classical MD. Ab-initio MD solves the many-body problem for a system governed by the Schrödinger equation with the Born-Oppenheimer approximation (without this approximation, the Schrödinger equation cannot be solved even numerically for many systems). At every timestep, the Schrödinger equation for electrons is solved to compute the potential. The position and velocity of every nucleus are then propagated by numerically integrating the classical equations of motion. To solve the Schrödinger equation for electrons, various approximations derived with methods such as the tight-binding model or density functional theory are used. Ab-initio MD can simulate quantum mechanical phenomena, such as bonding between atoms and heat conduction in metals, but the computational cost limits the size of systems that can be studied with this model to a few thousand atoms.
Classical MD simplifies the ab-initio approach by replacing the quantum mechanical potential computations with a parameterized empirical potential function that depends only on the position of nuclei. This potential function is either fitted to experimental data or derived from the underlying quantum mechanical principles. Classical MD can be used to study systems with up to millions of atoms but cannot simulate phenomena that depend on electron behavior. We will restrict our discussion to classical MD.
The name molecular dynamics is misleading because given the almost universal nature of the many-body problem, MD can model any system that involves interaction between particles. MD is employed in materials science to study a variety of phenomena, including the behavior of materials under stress, crack propagation, and how various defects affect the strength of materials. MD is also commonly used in biochemistry to study the behavior of macromolecules. The techniques of MD have also been generalized to thermal science and astrophysics, where particle models are used to study various hydrodynamic instabilities and phase transitions in thermal science and the structure of the universe in astrophysics.2
Interaction computations
The computation of interactions between atoms is typically the most computationally expensive segment of a MD simulation. For a pair-wise potential (a potential described by a potential function that does not account for the environment of atoms), the interaction computations scale with order O(N2). Various methods have been developed that reduce the order to O(N) for short-range potentials (potentials that decay faster with respect to distance than the dimension of the system are considered short-range) and O(NlogN) for long-range potentials. The most common methods are the linked-cell and neighbor list methods for short range potentials and Ewald summation method for long-range potentials.
Both the linked-cell and neighbor-list method assume that the potential and its spatial derivative at the location of an atom can be approximated by the superposition of the potential contributions from atoms within a certain interaction cutoff distance. Mathematically, the potential function U(rij) can be approximated by the truncated potential function Û(rij).
But to compute the distance between every atom pair costs O(N2) calculations. The linked-list method decomposes the system spatially into cells of side length greater than or equal to the cutoff distance. Every atom only interacts with atoms in the same cell or an adjacent cell so only distances between atom pairs in the same and adjacent cells must be computed. Atoms are resorted into cells every 10 to 20 timesteps. This reduces the complexity of the computation to O(N) .2
The linked-list method is commonly used in parallel molecular dynamics codes because it simplifies communications between threads. Every thread is assigned a cell and is passed a list of the atoms in its cell and every adjacent cell every time atoms are resorted into cells.
Figure 1 Linked-cell method: To compute the interactions of the filled in atom, only atoms in the same (dark-shaded) and adjacent (light-shaded) cells are considered. Only atoms within the shaded disk interact with the darkened atom.
The neighbor list method creates lists of all atoms within a certain distance of every atom so every atom only interacts with atoms in its neighbor list. The neighbor list cutoff distance must be greater than the interaction cutoff distance and the difference determines the frequency at which the neighbor lists must be recreated. The neighbor list method also reduces the complexity of the computation to O(N) (2).
Figure 2 Neighbor list method: The neighbor list for the filled-in atom contains all the atoms in the light-shaded disk. Only atoms within the dark-shaded disk interact with the filled in atom.
For long-range potentials, the interaction cutoff distance must be very large for the potential at the location of an atom to be approximated by the superposition of the potential contributions from atoms within the cutoff distance. The Ewald summation method is the most common method and reduces the computational complexity to O(Nlog(N)) 3 by decomposing the system spatially into cells and assuming the interactions of an atom with every other atom in the system can be approximated by its interactions with every other atom in the same cell and their images in periodic cells tessellated along all dimensions.
Mathematically, this can be expressed as:
where N is the number of atoms in the same cell, nx, ny, nz are half the number of periodic cells in the x, y, z directions respectively, Ø(rij) is the part of the function that does not decay with respect to distance, L is the side length of the cell and d is the dimension of the system. The Ewald summation method then decomposes the above sum into a short-range term that can be summed using techniques for short-range potentials and a long-range term that can be summed with Fourier methods. A detailed mathematical introduction to the Ewald summation method is beyond the scope of this article and we refer the reader to 3 for a detailed introduction to the Ewald summation method.
Integration methods
The solution to an n-body problem can be chaotic and it is fruitless to seek an accurate trajectory beyond several “collision steps”. The main criteria for choosing an integration method for the equations of motion in MD are therefore energy conservation, the ability to accurately reproduce thermodynamic ensemble properties and computational requirement. Interaction computations are the most computationally expensive segment of an MD simulation and an integration method requiring more than 1 interaction computation per timestep is inefficient unless it can maintain the same accuracy with a proportionally larger timestep. This criterion disqualifies the common Runge-Kutta methods.
Two commonly used families of integration methods are the Verlet algorithm and its variations and predictor-corrector methods. The Verlet algorithm can be derived from the Taylor expansion of the position variable (r):
Adding the expressions for r(t+h) and r(t-h) and solving for r(t+h):
One drawback of the Verlet algorithm is the lack of a velocity term in the expressions, making it difficult to compute the atomic velocities. Common variations of the Verlet algorithm designed to correct this drawback are the leap-frog algorithm and velocity-Verlet algorithm.
The predictor-corrector family of methods (PC methods) all consists of 3 steps. First, atomic positions and velocities are propagated according to the Taylor expansions of the position and velocity variables. Then, the interactions are computed based on the new positions and the accelerations are compared to the accelerations predicted with the Taylor expansion of the acceleration variable. Finally, the positions and velocities are “corrected” according to the difference between the computed and predicted accelerations. Although PC methods require 2 interaction computations per timestep, they can maintain the same accuracy with a timestep that is more than twice as long.3
PC methods have mostly been discarded in favor of the simpler Verlet methods because Verlet methods give better energy conservation and are easier to implement. We refer the reader interested in a derivation of the basic PC algorithm to. 1
Constraint dynamics
Geometric constraints can be implemented in the context of the above integration methods to simulate rigid molecules. The rigid molecule approach is not appropriate for large molecules. Methods to simulate flexible molecules have been developed although we will not introduce them in this article. We refer the reader to1 for a detailed introduction to the simulation of flexible molecules. The most common method is the SHAKE algorithm and its variations. The SHAKE algorithm first advances the atomic positions and velocities ignoring the constraints. The new positions and velocities are then corrected by introducing a constraint force, the magnitude of which is determined by the displacement of the new positions from their constrained positions.
The kth constraint constraining the distance between two atoms can be expressed mathematically as:
where ri and rj are the positions of atoms i and j and dk is the bond length between them. These constraint equations represent forces that are described by their action on the dynamics of the system and are not derived from the potential. These forces can be expressed within the equations of motion as constraint potentials:
where n is the number of constraints and λ is a constant that must be computed in order for the constraint to be satisfied. We numerically integrate the equations of motion with the Verlet algorithm.
where 
is the unconstrained position of atom i. ri and rj must obey the constraint equations σk:
where rk is the constrained bond length. The value of λk can be solved for with standard root-finding techniques. In the SHAKE algorithm, the system of k equations is solved with iterative methods [1].
Thermostats and barostats
MD with the classical equations of motion simulates a system at constant volume and energy. Every timestep, a new state of the system, called a microstate, is generated. Because every microstate generated has the same volume and energy, MD with the classical equations of motion can be thought of as sampling from a set of microstates with a given number of atoms, volume and energy, called a microcanonical/NVE ensemble (so called because the Number of particles, Volume and Energy are constant).
It is often desirable to sample from a canonical/NVT (Number of particles, Volume and Temperature are constant) or isothermal-isobaric/NPT (Number of particles, Pressure and Temperature are constant) ensemble to maintain compatibility with experimental studies. Various methods have been developed to constrain temperature and pressure and are called thermostats and barostats respectively. Most modify the system dynamics by modifying the equations of motion. Common thermostat algorithms include the Berendsen thermostat and Nosé-Hoover thermostat and common barostats include the Parrinello-Rahman method.
The equations of motion for most thermostat algorithms take the general form:
where γ can be thought of as a coefficient that adjusts the system kinetic energy to match the desired temperature. The effect of a thermostat algorithm on system dynamics can also be described as a scaling of internal velocities (velocities relative to center of mass) at every timestep by a factor λ. Mathematically, this can be expressed as:
where 
is the unscaled internal velocity and h is the timestep.
Because the temperature cannot change instantaneously, λ(t,0) = 1.
If we let 
in (*), then the above expression agrees with (*).
The Berendsen thermostat is based on Newton’s law of cooling:
where T0 is the heat bath temperature and τ is an empirical coupling parameter that quantifies the strength of the coupling between the system and heat bath. The finite difference quotient is:
The velocity scaling factor 
at every timestep is determined by the ratio between the system temperature and heat bath temperature and the coupling parameter.
Note that the coupling parameter affects the system dynamics like a damping parameter. The Berendsen thermostat, therefore, does not generate a true canonical velocity distribution because the coupling parameter damps out system temperature fluctuations.
Both the Nosé-Hoover themostat and the Parrinello-Rahman method are based on the extended system approach. In this approach, additional virtual degrees of freedom are introduced into the system that can be used to regulate system properties.
A detailed introduction to the Nosé-Hoover themostat and the Parrinello-Rahman method is beyond the scope of this article. We refer the reader 4 for a detailed introduction to the Nosé-Hoover thermostat and2
for a detailed introduction to the Parrinello-Rahman method.
Conclusion
According to American physicist Richard Feynman, “all things are made of atoms and everything that living things do can be understood in terms of the jiggling and wiggling of atoms”. MD is a computer simulation technique that simulates the “jiggling and wiggling of atoms” and has been hailed as a way of predicting the future by animating nature’s forces. In this article, we introduced the basic method behind MD to provide the reader with the necessary background knowledge to understand an MD simulation.
References
[1]. Rapaport, D. C.The Art of Molecular Dynamics Simulations. Cambridge : Cambridge University Press, 2004.
[2]. Griebel, Michael, Knapek, Stephan and Zumbusch, Gerhard. Numerical Simulation in Molecular Dynamics: Numerics, Algorithms, Parallelization, Applications. Berlin : Springer, 2007.
[3]. Leach, A. R.Molecular Modeling: Principles and Applications. Upper Saddle River : Prentice Hall, 2001.
[4]. Frenkel, Daan and Smit, Berend.Understanding Molecular Simulation: From Algorithms to Applications. San Diego : Academic Press, 2002.
[5]. Hünenberger, Philippe. Thermostat Algorithms for Molecular Dynamics Simulations. Advanced Computer Simulation. Berlin : Springer, 2005, pp. 105-149.