Molecular Dynamics: An Introduction, Lloyd Fosdick Sept. 5, 1995

http://www.sissa.it/furio/md/md

http://cmm.info.nih.gov/modeling/guide_documents/molecular_dynamics_document.html

Analysis of Lennard-Jones Clusters: Methods and Results, me 1998

In the broadest sense, Molecular Dynamics (MD) is concerned with simulating the motion of molecules to gain a deeper understanding of chemical reactions, fluid flow, phase transitions, droplet formation, and other physical phenomena that derive from molecular interactions. These studies include not only the motion of many molecules as in a fluid, but also the motion of a single large molecule consisting of hundreds or thousands of atoms, as in a protean. Motion is inherent to all chemical processes. Simple vibrations, like bond stretching and angle bending, give rise to IR spectra. Chemical reactions, hormone-receptor binding, and other complex processes are associated with many kinds of intra- and intermolecular motions. Much of this work uses simple classical Newtonian mechanics, most notably Newton's law:

For each atom i in a system constituted by N atoms. Here, m_i is the atom mass, a_i = d²r_i/dt² its acceleration, and F_i the force acting upon it, due to the interactions with other atoms. Therefore, in contrast with the other methods, molecular dynamics is a deterministic technique: given an initial set of positions and velocities, the subsequent time evolution is in principle completely determined. In a pictorial sense, virtual atoms will "move" in the computer. They will bump into each other, wander around (if the system is fluid), oscillate in waves in concert with their neighbors, perhaps even evaporating away from the system if there is a free surface, in a way very similar to what atoms in a real substance would do.

Statistical physics is the link between the microscopic behavior and thermodynamics. In the limit of very long simulation times, one could expect the phase space to be fully sampled, and in that limit this averaging process would yield the thermodynamic properties. In practice, the runs are always of finite length, and one should exert caution to estimate when the sampling may be good ("system at equilibrium") or not. In this way, MD simulations can be used to measure thermodynamic properties and therefore evaluate, say, the phase diagram of a specific material.

The forces that the atoms of the system exert on each other as described by Newton’s equation govern the rate and direction of motion (velocity). In practice, the atoms are assigned initial velocities that conform to the total kinetic energy of the system, which in turn, is dictated by the desired simulation temperature. Slowly "heating" the system (initially at absolute zero) and then allowing the energy to equilibrate among the constituent atoms carries this out. The basic ingredients of molecular dynamics are the calculation of the force on each atom, and from that information, the position of each atom throughout a specified period of time (typically on the order of picoseconds = 10^-12 seconds).

The force on an atom can be calculated from the change in energy between its current position and its position a small distance away. This can be recognized as the derivative of the energy with respect to the change in the atom's position:

Energies can be calculated using either molecular mechanics or quantum mechanics methods. Molecular mechanics energies are limited to applications that do not involve drastic changes in electronic structure such as bond making/breaking. Quantum mechanical energies can be used to study dynamic processes involving chemical changes.

Knowledge of the atomic forces and masses can then be used to solve for the positions of each atom along a series of extremely small time steps (on the order of femtoseconds = 10^-15 seconds). The resulting series of snapshots of structural changes over time is called a trajectory. The use of this method to compute trajectories can be more easily seen when Newton's equation is expressed in the following form:

In practice, trajectories are not directly obtained from Newton's equation due to lack of an analytical solution. First, the atomic accelerations are computed from the forces and masses. The velocities are next calculated from the accelerations based on the following relationship:

Lastly, the positions are calculated from the velocities:

Two states are connected by some reaction path and a trajectory between them (i.e. r_i -> r_j) can be subdivided into a series of sub-states (r_k) separated by a small time step, "delta t" (e.g. 1 femtosecond).

The initial atomic positions at time "t" are used to predict the atomic positions at time "t + delta t". The positions at "t + delta t" are used to predict the positions at "t + 2*delta t", and so on.

The "leapfrog" method is a common numerical approach to calculating trajectories based on Newton's equation. The steps can be summarized as follows:

The method derives its name from the fact that the velocity and position information successively alternate at 1/2 time step intervals.

MD has no defined point of termination other than the amount of time that can be practically covered. Unfortunately, the current picosecond order of magnitude limit is often not long enough to follow many kinds of state to state transformations, such as large conformational transitions in proteins.

Beyond this "traditional" use, MD is nowadays also used for other purposes, such as studies of non-equilibrium processes, and as an efficient tool for optimization of structures overcoming local energy minima (simulated annealing).

Molecular dynamics calculations can be performed using AMBER, CHARMM, CHARMM/GAMESS, Discover, QUANTA/CHARMm, and SYBYL.

Computers are a critically important tool in the study of MD because there simply is no other way to trace the motion of a large number of interacting particles. The earliest of these computations were done in the 1950s by Berni Alder and Tom Wainwright at Lawrence Livermore National Laboratory. They studied the distribution of molecules in a liquid, using a model in which the molecules are represented as hard spheres, which interact like billiard balls. Using the fastest computer at that time, and IBM 704, they were able to simulate the motions of 32 and 108 molecules in computations requiring 10 to 30 hours. Now it is possible to perform hard sphere computations on systems of over a million particles. A hard sphere model, with millions of molecules, has been used at NASA Ames Research Center by Leonardo Dagum to simulate hypersonic flow conditions encountered by flight vehicles at high altitudes. Computations using a more realistic molecular model known as Lennard-Jones have been performed at IBM Kingston to study the flow of fluids. In these computations approximately approximately 10,000 interacting molecules represented the fluid. Even though this is miniscule compared with the number of molecules in a gram of water, the behavior of the flow was like that in a real fluid.

Reviewing even a tiny subset of the applications of molecular dynamics is far beyond the purpose of these notes. I will only briefly mention a few areas of current interest where MD has brought and/or could bring important contributions. The list should not be considered at all exhaustive:

Liquids.

Defects.

Fracture.

Surfaces.

Friction.

Clusters.

Biomolecules.

Electronic properties and dynamics.

Models of particle systems are characterized by the nature of the interactions between the particles. Generally it is assumed that the forces between the particles are conservative, two-body forces; that is, energy is conserved and the total force acting on a particle due to the other particles is the sum of the forces between pairs of particles. Thus the force acting on particle i is given by an expression of the form:

Where f_i is the total force on particle i due to the other particles, f_i,j is the force on particle i due to particle j, and n is the number of particles in the system. Force is a vector quantity, so the sum in equation (6) is a vector sum. The order of the indices is important: the first index identifies the particle action on, the second index identifies the particle causing the action. Newton’s third law states that:

There is an important relation between potential energy and force in a conservative system. If r is the position of a particle, f(r) the force acting on it, and f (r) its potential energy, then

Thus we can describe a model in terms of the force or the potential energy. Since potential energy is a scalar quantity it is often more convenient to describe the model in terms of its potential energy function, f .

At a point of minimum potential energy the partial derivatives of the potential energy function are zero, and it is a point at which all of the forces are zero. Accordingly we call this point an equilibrium point.

The Lennard-Jones pairwise potential is a simple yet reasonably accurate mathematical model of physical low temperature collections of heavy rare gas atoms such as Argon or Xenon. Such atoms are chemically inert and do not usually enter into chemical bonds, but rather interact primarily through a van der Waals force. Thus computational results from the mathematical model can be compared directly with laboratory experiments on low temperature rare gas atoms. Consequently a considerable literature has arisen regarding LJ clusters, particularly with respect to the global optimization problem of determining the lowest energy configuration as a function of N (the number of atoms in the cluster).

Collections of atoms interacting pairwise via the LJ pair potential are called Lennard-Jones clusters

The Lennard-Jones 12-6 potential energy function between two particles i and j is:

As depicted in the following graph.

This potential has an attractive tail at large r, it reaches a minimum around 1.122s , and it is strongly repulsive at shorter distance, passing through 0 at r = s and increasing steeply as r is decreased further.

The term ~1/r¹², dominating at short distance, models the repulsion between atoms when they are brought very close to each other. Its physical origin is related to the Pauli principle: when the electronic clouds surrounding the atoms starts to overlap, the energy of the system increases abruptly. The exponent 12 was chosen exclusively on a practical basis. In fact, on physical grounds an exponential behavior would be more appropriate.

The term ~1/r⁶, dominating at large distance, constitute the attractive part. This term gives cohesion to the system. The 1/r⁶ attraction is originated by van der Waals dispersion forces, originated by dipole-dipole interactions in turn due to fluctuating dipoles. These are rather weak interactions, which however dominate the bonding character of closed-shell systems, that is, rare gases such as Ar or Kr. Therefore, these are the materials that a LJ potential could mimic fairly well. The parameters e and s are chosen to fit the physical properties of the material.

On the other hand, a LJ potential is not at all adequate to model situations with open shells, where strong localized bonds may form (as in covalent systems), or where there is a delocalized "electron sea" where the ions sit (as in metals). In these systems the two-body interactions scheme itself fails very badly

However, regardless of how well it is able to model actual materials, the LJ 12-6 potential constitutes nowadays an extremely important model system. There is a vast body of papers who investigated the behavior of atoms interacting via LJ on a variety of different geometries (solids, liquids, surfaces, clusters, two-dimensional systems, etc). One could say that LJ is the standard potential to use for all the investigations where the focus is on fundamental issues, rather than studying the properties of a specific material. The simulation work done on LJ systems helped us (and still does) to understand basic points in many areas of condensed matter physics, and for this reason the importance of LJ cannot be underestimated.

The total potential of a system containing many atoms is then given by:

When using the LJ potentials in simulation, it is customary to work in a system of units where s = 1 and e = 1. Thus the following "simplified" Lennard-Jones potential energy function is often used:

This "simplified" model works as well as the more complicated model. The potential has a single well which reaches a minimum at r = 1. For r < 1 the potential is repulsive (negative derivative), and the potential energy increases rapidly toward positive infinity as r -> 0. This repulsive force keeps the two atoms from approaching each other at distance significantly less than r = 1. For r > 1, the potential is attractive. The equilibrium separation between two atoms occurs at r_eq = 1, at which point the attractive and repulsive forces are in balance and the potential energy is minimized. Note that coefficients of the two terms in the LJ potential have been selected so that the minimum occurs at r_eq = 1 with u(r_eq) = -1. Other choices are possible; however, it is easily shown that as long as the coefficient of the r^-12 term remains positive and that of the r^-6 term negative, such variations simply rescale the u and r coordinate axes, but do not otherwise fundamentally change the potential.

A problem of major interest in computational science is the determination of the stable configuration states, and in particular, the lowest energy configuration, of systems of interacting particles. In chemistry, this is sometimes referred to as the molecular conformation problem, where the particles correspond to atoms and the interactions are the various forces involved in the formation of chemical bonds as well as non-bonded interactions such as electrostatic and van der Walls forces. Determining the energetically optimal geometric configuration of a molecule is often key to understanding its chemical interactions.

Here we consider one of the simplest molecular conformation problems that of N identical atoms interacting via the Lennard-Jones model. In the model, let atom i have coordinates r_i=(x_i, y_i, z_i) in Euclidean 3-dimensional space. The interaction between atom i and j is assumed to be governed by a pair potential energy function of the form

u = u(r_i,r_j) = u(|r_i - r_j|) = u(r_ij)

where and the pair potential u(r_ij) between atoms i and j is a function only of the Euclidean distance between the atoms, and does not depend on the orientation or absolute position of the atoms (i.e. the potential function is invariant with respect to rotations and translations that act on both atoms simultaneously in 3-dimensional space).

The configuration of all N atoms can be represented by the 3N-dimensional coordinate vector r = (r₁, r₂, …, r_N). By the additivity assumption, the total potential energy function U(r) is the sum of all N(N-1)/2 pairwise interactions, i.e.

The function U(r) is often referred to as the potential energy surface (PES).

The 3N-dimensional vector of forces f(r) acting on all N atoms is given by

f(r) = -U(r)

Of particular importance are the stationary, or equilibrium, points of U(r) where the gradient, and hence the net force on each atom, vanishes. These include local minima and maxima as well as saddle points. Of these, only local minima are mechanically stable with respect to small displacements, and the location of the local minima is perhaps the most important feature of the PES.

In the previous section, we introduced the cluster potential energy function U(r) and the pair potential energy function u(r_i,r_j). In the discussions below, we shall be dealing with a mixture of scalars (real numbers), finite dimensional vectors with real components, and matrices with real elements. The following notational conventions will be used:

Scalars.

Vectors and Matrices.

Norm.

Where g = g(r₀), H = H(r₀), and O(|r-r₀|³) denotes cubic and higher order terms.

Let f : R -> R be a differentiable function of a real variable x, and let q be a constant. Then by the chain rule for differentiation,

By differentiating we obtain the formula for the second derivative in the univariate case:

The analogous formula for the gradient of a function of a quadratic form x^TQx, where x is an N-dimensional vector and Q is an N-by-N symmetric square matrix:

f(x^TQx) = 2f'(x^TQx)Qx

The analogous multivariate formula for the Hessian is

H(f(x^TQx)) = 2f'(x^TQx)Q + 4f''(x^TQx)Qxx^TQ

These last two equations provide a very elegant and compact matrix-vector representation of the gradient and Hessian of such functions without explicitly resorting to partial derivatives. Moreover, they are easily implemented in high-level computer languages. They are directly relevant to Lennard-Jones clusters because, as shown below, the Lennard-Jones pair potential function u(r_i,r_j)can be written as a function of a quadratic form in the vector(r₁,r₂), namely the square of the distance between atoms i and j.

Let be a 6-by-6 matrix composed of 3-by-3 identity matrices I. Then define

Thus the quadratic form q is the square of the distance between atoms i and j, and the pair potential can be written as:

with first and second derivatives f'(q) = -6(q^-7 - q^-4) and f''(q) = -42q^-8 - 24q^-5.

Thus by applying the quadratic form of the chain rule, we have that the gradient and Hessian of the pair potential are given by

By the linearity of partial differentiation, we can formally write

where u(r_i,r_j) is regarded as a 3N-dimensional vector equal to zero everywhere except in the 6 components corresponding to atoms i and j. Similarly, we have

where is a 3N-by-3N matrix obtained by augmenting the 6-by-6 matrix computed with zeros in the elements not corresponding to interactions between atoms i and j.

Generally the Lennard-Jones potential function U(r) is quite well behaved, and almost any standard multivariate optimization technique can be used to find local minima. Since we have developed analytic formulas for the gradient and Hessian, both first order (steepest decent and conjugate gradient) as well as second order (Newton-Raphson) methods can easily be applied. We sketch the three main local-minimization scenarios below:

1. Starting from a transition point, using a steepest descent algorithm to walk down an approximate 'reaction path' connecting the two adjacent local minima;
2. Starting close to a local minimum, using second-order information in the form of a Newton-Raphson step to converge quadratically to the local minimum, and
3. Starting from a random point in configuration space, generally not close to a local minimum or transition point, using a conjugate gradient technique to converge to a local minimum.

The following algorithm is used to 'walk down' from a transition point r₀ to its two connected adjacent local minima. Let v be the eigenvector of H(r₀) corresponding to the single negative eigenvalue, where v is assumed to have norm equal to 1. For sufficiently small values of |b|,

f(x + bv) < U(r₀)

since U(r₀) is a local maximum with respect to the normal coordinate along v. Starting at r₁ = x + bv for a small b > 0, a steepest (gradient) descent step of length b along the unit negative gradient vector is taken in the iteration

For sufficiently small positive values of b, this will be a descent step as long as r_old is not a stationary point.

The path generated by iterating approximates the 'reaction path' defined by the differential equation

In general, this fixed stop length algorithm will not converge exactly to a local minimum since the iterates are always a fixed distance apart. Rather, the algorithm will follow a descending path to the vicinity of a local minimum and is terminated when a step is attempted which no longer results in a lower value of U(r). At this point a Newton-Raphson algorithm (see below) can be applied to converge to the local minimum within a tight tolerance.

We have found, by experimentation that b = 0.01, generally works well for small Lennard-Jones clusters in the size range N £ 147. The steepest descent phase usually terminates at a point from which two or three Newton-Raphson iterations will converge to a local minimum which has converged to at least 10 significant digits in the function value U(r) and for which |g(r)| < 10^-8.

The procedure is then repeated using a starting point of r = r₀ - lv to reach the other adjacent local minimum connected to r₀.

From the multivariate Taylor series; if expansion around a point r₀ is terminated after the quadratic term and differentiated, we obtain

g(r) = g(r₀) + H(r₀)(r-r₀)

Since we are looking for a minimum (i.e. where the gradient vanishes), setting g(r) = 0 immediately lead to the Newton-Raphson iteration implicitly defined by the 3N-by-3N set of linear equations

H(r_old)(r_new-r_old) = -g(r_old)

It can be shown that under mild assumptions (H is invertible) this iteration converges quadratically to a local minimum if started sufficiently close to that point (i.e. in a region where the Taylor's series truncation is quite accurate). Similarly, if started sufficiently close to a transition point, the Newton-Raphson method will converge quadratically to that transition point.

Note that H(r) is in general not invertible due to the presence of zero eigenvlaues corresponding to rigid translations. However it follows from conservation of energy and angular momentum that g(r) has no components in the null space of H and hence the system of equations is consistent.

Another standard multivariate optimization technique commonly used for potential energy function minimization is the Polak-Ribier variation of the non-linear conjugate gradient (CG) minimization algorithm. This iterative technique performs successive 1-dimensional minimizations over a step direction defined by a linear combination of the gradient at the current point and the step direction at the previous point. Here we use the same convergence criterion as with the Newton-Raphson method above(|g(r)| < 10^-8).

The CG technique has been successfully used by many investigators as a local optimization algorithm for the LJ PES. Our experience with LJ clusters has been that CG generally works well even when started from a random point quite far from a local minimum. Occasionally the CG algorithm will converge from such a random start to a saddle point rather than a local minimum, so it is necessary to inspect the Hessian at convergence. Also, we have observed that very occasionally the CG technique will fracture the cluster (i.e. create more than one sub-cluster for the conformation), usually resulting in a configuration of a single isolated atom and a cluster of N-1 atoms. This condition is easily detected by examining the individual pair potential terms in the summation forming U(r).