Numerical approaches
So far we have discussed the mathematical theories behind a molecular dynamics simulation. However, in order to be able to implement such a method into a program code, several numerical approaches have to be made. In this section the integrator schemes, periodic boundary conditions, the truncating of potentials and the conservation of linear momentum.
Integration schemes
Several numerical methods exist to solve the Newtonian equations of motion, the most popular being the Verlet algorithm. Often even an adapted version of this algorithm is used, called the leapfrog algorithm. Both algorithms are discussed briefly in this section.
To derive the Verlet algorithm we start with the Taylor series of the position of a particle around time t and with time step Δt
 | (1) |
Of course the first order derivative of r(t) with respect to t is equal to the velocity. Similarly the second order derivative equals the acceleration, which can be substituted by the force by using Newton’s second law. This reduces (1) to
 | (2) |
and similarly we have
 | (3) |
Summing these two equations and neglecting the higher order terms (Δt4, Δr4) gives us
 | (4) |
where r(t+Δt) is the new position of the particle, r(t) its present position, r(t−Δt) the position on the previous time step, F(t) the resulting force on the particle, and m the mass of the particle. Note that the Verlet algorithm does not use the velocity of a particle to compute its new position.
The difference between the Verlet algorithm and the leapfrog algorithm is that the latter evaluates the velocities at half-integer time steps, which are artificial points in time, exactly half way between two time steps. The algorithm uses these half-integer time steps velocities to compute the new positions. To derive the leapfrog algorithm from the Verlet algorithm let us first define the velocity at half-integer time steps
 | (5) |
and
 | (6) |
Using the above equations and the Verlet algorithm as expressed by (4), we can obtain an expression the velocity on time t+Δt/2
 | (7) |
and therefore it is from (6) that we immediately obtain an expression for the new positions
 | (8) |
It must be understood that instead of having the velocities at the same time as the positions there now is a difference of Δt/2 between the positions and velocities. Still, during the simulation the positions and velocities are treated as if they were at the same time. The elimination of the velocities from the leapfrog algorithm shows that the method is algebraically equivalent to Verlet’s algorithm, and, hence, it gives identical trajectories [1]. There are some advantages in using the leapfrog scheme. Numerical benefits derive from the fact that at no stage the difference of two relatively large quantities to obtain a small one (this happens in (4) in the Verlet algorithm). This adaptation minimizes the numerical error [2].
The steps to take in solving the equations of motions using the leapfrog algorithm are summarized in the list below. It is assumed that the initial positions, r(t), and velocities, v(t−Δt/2), are known.
- Calculate F(t).
- Calculate new velocities v(t+Δt/2) using v(t−Δt/2) and F(t).
- Calculate new positions r(t−Δt/2) using r(t) and v(t−Δt/2).
- Repeat these steps until a certain stop criterium.
The time steps associated with the Verlet or the leapfrog algorithm have to be chosen in such a way that the integration scheme resembles the motion of the particles correctly. When the time step is chosen too small the movement of the particles is very small and the simulation can take a very long time. However, more drastically, is the problem when the time steps are chosen too large, resulting in instabilities, such as extreme collisions. In the figure below this effect is illustrated schematically. For instance, if one is interested in the vibrational motions inside a molecule, a suitable time step is around 1 fs, but for more global motions a time step of several femtoseconds is not unusual.
 |
| The choice of an appropriate time step is important to prevent motions that develop too slow (left) or motions leading to instabilities (middle). With the appropriate time step collisions occur smoothly (right). |
Periodic boundary conditions
In order to conserve the macroscopic behavior of the system under investigation, it is important that the boundaries of the system are treated in a correct way. Several methods to treat the boundaries exist, but the most common are harmonic boundary conditions and periodic boundary conditions. Since it is most common to perform molecular dynamics simulations with the latter type, there are discussed in more detail.
When using periodic boundary conditions we divide the simulation universe into boxes, which are all exact replicas of the central box to form an infinite lattice. In our simulations we only look at the molecules in the central, or original, box. If, in the course of the simulation, a particle moves in the central box, its periodic image in each of the neighboring boxes moves in exactly the same way. Thus, if a particle leaves the central box, one of its images will enter trough the opposite face. Hence, there are in fact no walls at the boundaries of the central box. The number of particles in the central box is, therefore, conserved. The box simply forms a convenient axis system for measuring the coordinates of the N particles [2]. Fortunately it is not needed to store the coordinates of all images of the particles, which would be an infinite set, since when a particle leaves the central box by crossing a boundary, attention is switched to the image, which has just entered the central box.
Simply said, periodic boundary conditions mean that the system is its own neighbor in all dimensions. One can see this either as an universe consisting of an infinite amount of similar systems, or that the universe is the system by itself.
 |
The most simple way for implementing periodic boundary conditions is to use the cubic cell, but different shapes are possible (such as the hexagonal prism cell or the rhombic dodecahedron cell). In the figure above the two-dimensional representation of the periodic boundary conditions for a cubic cell is depicted.
Not only the movement of the particles is limited by these periodic boundary conditions, also the interactions between the particles are treated in the same way. It is, therefore, important to ensure that the size of the box is chosen in such a way that the particle cannot ‘sense’ its own image in the neighboring box, and, hence, influences its own behavior.
Furthermore, when calculating the interaction between a certain particle and any other particle, one has to make sure to take the correct image of the other particle, in such a way that the smallest distance between the particle and any of its images is calculated. This prerequisite is known as the minimum image convention. In the minimum image convention, each particle ’sees’ at most just one image of every other particle in the system (which is repeated infinitely via the periodic boundary
conditions), and the interaction is calculated with the closest particle or image.
Using periodic boundary conditions and the minimum image convention ensures us the macroscopic behavior of the system is maintained, if we have chosen the size of the box properly. One major drawback, however, is that fluctuations having a wavelength larger than the length of the box cannot be achieved, since they do not fit in the periodic lattice. However, periodic boundary conditions model most systems well.
Truncating potentials
The most time consuming part in a molecular dynamics simulation with N particles is the calculation of the non-bonded interactions. The number of intramolecular interactions (for instance bonding) is proportional to the number of particles in the system. In contrast, the number of non-bonded interactions increases as the square of the number of particles, since in general the non-bonded interactions are calculated between every pair of particles in the system, and is therefore of order N2.
However, if we take the Lennard-Jones potential, we can see it falls off very rapidly when the separation between two particles increases. At a distance of already at 2.5 σij, the potential is almost zero. It is easily understood that when calculating the potential at distances much larger than 2.5 σij computational power is wasted, for the type of systems we are investigating. Therefore, it is clever to introduce a method by which we do not have to take the contribution beyond, for instance 2.5 σij, into account.
The most popular way to achieve this is by introducing a so-called cutoff distance. When this cutoff is employed, all the interactions between all pairs of particles that are further apart than a specified cutoff value are set to zero, taking into account the restrictions applied by the periodic boundary conditions.
However, the implementation of cutoffs also introduces new problems, such as a discontinuity in the potential, since the potential energy instantaneously changes to zero at the cutoff distance. As a consequence the total energy cannot be conserved, which is required in most simulations. There are several ways to deal with this problem, such as the use of shift functions.
The shift function simply shifts the potential in such a way that the discontinuity vanishes. This is done by adding a constant term from the potential. While this approach solves the problem with the energy conservation it does not solve all problems, since the force (being the first derivative of the potential energy function) still will have a discontinuity, which can make the system instable, i.e.
result in unphysical behavior. However, since the contribution at the cutoff distance is negligible, the discontinuity is assumed to have no influence.
By itself, the use of a cutoff and a shift function will not dramatically reduce the time to calculate the non-bonded interactions. This is because of the simple fact that we have to calculate the distance between every pair of particles in the system to determine whether they lie within the cutoff distance or not. This is almost as time consuming as calculating all the non-bonded interactions themselves.
To finally solve the problem of calculating the N·(N−1)/2 distances, the non-bonded nearest neighbor list was proposed by Verlet [3]. In this approach all particles within the cutoff distance of one certain particle, together with all particles that are just a little bit further away (at a distance called the Verlet radius) are stored in a list. Now only the non-bonded interactions between the particle and its neighbors
stored in the list have to be calculated. The number of non-bonded interactions to be calculated therefore is no longer equal to N·(N−1)/2 but scales on the average linearly with N.
However, the non-bonded nearest neighbor list has to be updated with regular intervals, since two particles are not likely to remain in close proximity during the entire simulation. As mentioned earlier, not only the particles within the cutoff distance are stored in the list, but also particles within the Verlet radius. This is done to make sure no particles will suddenly move into the cutoff distance without having the knowledge that they are approaching and therefore introduce unwanted effects into the simulation. This extra shell surrounding the cutoff radius has to be chosen of such a thickness that no particle approaches closer than the cutoff distance before the list is updated. It is clear that the size of this shell and the update frequency are strongly correlated.
There still is a small problem when updating the list: to determine which particles are in close proximity all inter-particle distances have to be calculated. For a large system this still is very time consuming. Therefore, it is normal to divide a cubic or rhombic simulation box in a number of cells. The length of each of these cells should be longer than the cutoff distance. If the system is divided into
M3 cells, there will be an average of N/M3 particles in each cell. To find the new neighbors only the cell in which the particle is positioned and its surrounding cells (26 in a 3D-system) are subject to the calculation of the distances. It is then necessary to consider just 27N/M3 particles on average rather than N, which clearly is a reduction of computational time.
Linear momentum conservation
In most cases we do not want an external force acting on our system, and thereby influences its behavior. Therefore, we require that the external force equals zero and, hence, that dp/dt=0, from which it immediately follows that the linear momentum should remain constant throughout a simulation. For convenience we require that the linear momentum equals
 | (9) |
This requirement is known as the conservation of linear momentum.
The choice of keeping the linear momentum of the system at zero is arbitrary, but for our purposes it is very useful, since if we allow the system to have a nonzero, but constant, linear momentum, the center of mass of the system moves with a constant speed. When subsequently calculating properties using the velocities of the particles we have to subtract the velocity of the center of mass to obtain the correct properties. When setting the linear momentum of the system at zero, and, therefore, the center of mass does not move, we can use the velocities of the particles immediately to obtain properties, without doing some correction first. This is best illustrated by an example. If we would fix all particles at a certain lattice (so they cannot move with respect to each other), but if we set the linear momentum at a nonzero, but constant value, the system would move as a rigid body (so all particles
have exactly the same velocity vector, corresponding to the vector implied by the nonzero linear momentum). However, if would like to calculate the temperature of the system (which depends on the kinetic energy only, and, hence, on the velocities of the particles), we would not obtain a temperature of zero, which would be correct, since the particles cannot move with respect to each other, but a temperature
imposed by the nonzero linear momentum. This is clearly not what we want. Therefore we require that the linear momentum of the system should be zero. One can compare the example with an ice cube flying trough the air, although the ice cube has a velocity as rigid body, its temperature remains constant.
[1] D. Frenkel and B. Smit, Understanding Molecular Simulation - From Algorithms to Applications, Academic Press, San Diego, USA (1996).
[2] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, 1st edition, Oxford University Press, Oxford, England (1987).
[3] L. Verlet, Computer ‘Experiments’ on Classical Fluids - Equilibrium Correlation Function, Phys. Rev., 165, 201–204 (1967).