Integrators in Molecular Dynamics

Potential energy is a function of the atomic positions (3N) of all the atoms in the system. Due to the complicated nature of this function, there is no analytical solution to the equations of motion they must be solved numerically. Numerous numerical algorithms have been developed for integrating 

All the integration algorithms assume that the positions, velocities, and accelerations can be approximated by a Taylor series expansion  

---

Tabor, M. (1989). Chaos and Integrability in Molecular Dynamics (Wiley, New York). 

T. P. Straatsma and J. A. McCammon, Ghem. Phys. Lett., 167, 252 (1990). Free Energy Thermodynamic Integrations in Molecular Dynamics Simulations Using a Non-iterative Method to Include Electronic Polarization. 

Integrators in molecular dynamics simulations are supposed to be accurate, i.e., they should enforce the exact trajectory being followed as closely as possible. They should provide stability, meaning that the constants of motion, e.g., the total energy in the microcanonical ensemble, are preserved. Nevertheless, the integrators should be efficient, which means that a minimum number of force calculations are needed in order to save computer time. The best numerical methods are based on... 

In this paper, we discuss semi-implicit/implicit integration methods for highly oscillatory Hamiltonian systems. Such systems arise, for example, in molecular dynamics [1] and in the finite dimensional truncation of Hamiltonian partial differential equations. Classical discretization methods, such as the Verlet method [19], require step-sizes k smaller than the period e of the fast oscillations. Then these methods find pointwise accurate approximate solutions. But the time-step restriction implies an enormous computational burden. Furthermore, in many cases the high-frequency responses are of little or no interest. Consequently, various researchers have considered the use of scini-implicit/implicit methods, e.g. [6, 11, 9, 16, 18, 12, 13, 8, 17, 3]. 

In molecular dynamics, successive configurations of the system are generated by integrating Newton s laws of motion. The result is a trajectory that specifies how the positions and velocities of the particles in the system vary with time. Newton s laws of motion can be stated as follows ... 

There are many algorithms for integrating the equations of motion using finite difference methods, several of which are commonly used in molecular dynamics calculations. All algorithms assume that the positions and dynamic properties (velocities, accelerations, etc.) can be approximated as Taylor series expansions ... 

Fincham D and Heyes D M 1982. Integration Algorithms in Molecular Dynamics. CCP5 Quarterly 6A 10. 

The above formula for Z, the NPT partition function, was first reported by Guggenheim [74], who wrote the expression down by analogy rather than based on a detailed derivation. While this form of the partition function is thought to be broadly valid and is widely applied (for example in molecular dynamics simulation [6]), it introduces the conceptual difficulty that the meaning of the discrete volumes Vi is not clear. Discrete energy states arise naturally from quantum statistics. Yet it is not necessarily obvious what discrete volumes to sum over in Equation (12.50). In fact for most applications it makes sense to replace the discrete sum with a continuous volume integral. Yet doing so results in a partition function that has units of volume, which is inappropriate for a partition function that formally should be unitless. 

The motions of proteins are usually simulated in aqueous solvent. The water molecules can be represented either explicitly or implicitly. To include water molecules explicitly implies more time-consuming calculations, because the interactions of each protein atom with the water atoms and the water molecules with each other are computed at each integration time step. The most expensive part of the energy and force calculations is the nonbonded interactions because these scale as 77 where N is the number of atoms in the system. Therefore, it is common to neglect nonbonded interactions between atoms separated by more than a defined cut-off ( 10 A). This cut-off is questionable for electrostatic interactions because of their 1/r dependence. Therefore, in molecular dynamics simulations, a Particle Mesh Ewald method is usually used to approximate the long-range electrostatic interactions (71, 72). 

Molecular dynamics (MD) is a long standing molecular modelling technique [23], In molecular dynamics we solve by numerical integration the equations of motion of a system of interacting particles with or without periodic boundary conditions. 

A variety of algorithms have been used for integrating the equations of motion in molecular dynamics simulations of macromolecules. Most widely employed are the algorithms due to Gear91 and Verlet.90 The algorithm introduced by Verlet in his initial studies of the dynamics of Lennard-Jones fluids is derived from the two Taylor expansions,... 

The plain Lennard-Jones potential is shifted up, so that its minimum located at 21/6ct has value 0, and set to zero beyond that point. The advantage of including the r 6 contribution instead of merely using the purely repulsive r 12 is that Eq. 4 is exactly zero beyond rcut and merges smoothly to this value at rCM. The use of a smooth hard core in molecular dynamics simulations is advantageous since the force is the derivative of the potential therefore the latter should be differentiable. In fact, the derivative must also be bounded to ensure numerical stability of the discrete integrator. 

The implementation of thermodynamic integration in molecular simulation calculations is fairly straightforward. For each of a range of discrete values of the coupling parameter X between 0 and 1, a molecular dynamics or Monte Carlo simulation needs to be carried out. From each of these simulations, the ensemble average 33 f(p, qN X)/dX) is evaluated. The free energy difference is then found from... 

In the previous chapter, we discussed the growth of error in numerical methods for differential equations. We saw that if the time interval is fixed, the error obeys the power law relationship with stepsize that is predicted by the convergence theory. We also saw that this did not contradict the exponential growth in the error with time (when the stepsize is fixed). The latter issue casts doubt on the reUance on the convergence order as a means for assessing the suitability of an integrator for molecular dynamics. 

---