Equations of motion for atomic systems

Molecular Dynamics

The equations of motion

Molecular dynamics requires a technique for the solution of the equations of motion for atomic systems. The most fundamental form of these equations of motion is the Lagrangian.

where the Lagrangian is defined in terms of the kinetic, K and potential, V energies, L = K - V. L is a function of generalized coordinates q_k and their time derivatives dq_k/dt If we consider a system of atoms, with Cartesian coordinates r_i such that K = 1/2m_i(¶r_i/¶t) Eqn. (1) becomes

where m_i is the mass of atom i and

is the force on that atom. These equations also apply to the center of mass motion of a molecule with f_i representing the total force on molecule i.

The generalized momentum p_k conjugate to q_k is defined as

An alternative formulation uses Hamilton’s equations

For Cartesian coordinates Hamilton’s equations become

The center of mass trajectory for the system is calculated by solving 3N second order differential equations (LaGrangian) or 6N first order coupled differential equations (Hamiltonian). If the potential function V depends only on the motion of the particles relative to one another and there is no external field applied then V, H, and L are independent of translation and rotation about the center of mass. In other words the total linear and angular momenta of the center of mass are conserved quantities for a completely isolated set of interacting molecules or if we are modeling a single macromolecule. However, if we use rectilinear periodic boundary conditions then angular momentum is not conserved while linear momentum is conserved.

Conservation and boundary conditions

Linear momentum is conserved in periodic boundary conditions.

Angular momentum is conserved for spherical boundary conditions. Energy is conserved in MD simulations (dH/dt = 0).

Time reversal applies to MD simulations.

Solution of the differential equations of motion by finite difference methods

To propagate molecular motion we need a means of solving the above equations of motion for a system of many particles. Coupled linear differential equations of the type described above for the motion of various masses in a forcefield can be solved using finite difference methods. The equations are solved step-by-step in discrete time intervals dt. Finite difference methods use calculation of the velocity (i.e. ¶r/¶t) times the time step to produce a new set of positions. Then the new positions are used to reevaluate the velocities using the equations of motion. This procedure is repeated for each step of the simulation.

The choice of time step is clearly of critical importance to the success of the method. The time step must be short relative to the length of time it takes for a particle to travel its own length. Another way to pose the criterion is that the time step should be about 10 times shorter than the period of the highest frequency vibration in the simulation. However, the configuration space sampled during the simulation will be greater if the time step is longer, so in the interest of efficiency of calculation it is desirable to make the time step as long as possible.

There are several different techniques for propagating the motion of the particles in a simulation. Allen & Tildesley discuss two methods.

Verlet algorithm

basic
leapfrog
velocity Verlet

Gear predictor-corrector

The Verlet algorithm

Basic Verlet

The Verlet method is a direct solution of the second order differential equations. In the Verlet method the velocities are eliminated by comparison of two expansions about the position at time t. The Taylor series expansion about +dt and -dt are summed to give the expression

r(t + dt) = r(t) + dt v(t) + (1/2)dt² a(t) + …

r(t - dt) = r(t) - dt v(t) + (1/2)dt² a(t) + …

-------------------------------

r(t + dt) = 2r(t) + r(t + dt) + dt² a(t) + …

This is equation is correct except for errors of the order of dt⁴. The computed velocity (used to estimate the kinetic energy) is subject to errors of the order of dt². The velocity is computed by v(t) = [r(t + dt) – r(t – dt)]2dt, on the fly, in this method.

Notice that the property that dynamics must be symmetric with respect to time reversal is used in the definitions required for the Verlet algorithm.

Leapfrog Verlet

The Verlet algorithm may introduce numerical imprecision since numbers of the order of dt² are added to numbers of the order dt⁰ (» 1). For this reason the leapfrog Verlet method is used

r(t + dt) = r(t) + dt v(t + 1/2dt)

v(t + 1/2dt) = v(t - 1/2dt)r(t) + dt a(t)

The velocity equation is executed first and generates a new mid-step velocity. This velocity is then used to calculate the new position. The velocity is calculated from

v(t) = (1/2)v(t + 1/2dt) +(1/2)v(t - 1/2dt)

This leapfrog method also has the advantage that temperature scaling by velocity scaling is feasible.

Velocity Verlet

The handling of the velocity (and therefore the calculation of the kinetic energy) is not ideal in either of the above forms of the Verlet algorithm. The velocity Verlet algorithm stores positions, velocities, and accelerations.

r(t + dt) = r(t) + dt v(t) + (1/2)dt² a(t)

v(t + dt) = v(t) + (1/2)dt[a(t) + a(t + dt)]

The above velocity Verlet approach can be shown to be equivalent to basic Verlet algorithm by eliminating the velocities. The equations are implemented in two stages. First, the new positions at time t + dt are calculated. Then the velocities at mid-step are calculated using

v(t + 1/2dt) = v(t) + (1/2)dt a(t)

The forces and acceleration at time t + dt are then calculated and then the new velocity is calculated.

v(t + dt) = v(t + 1/2dt) + (1/2)dt a(t + dt)

The Gear predictor-corrector method

If the classical trajectory is continuous, then an estimate of the positions, velocities, accelerations etc. may be obtained by a Taylor series expansion about time t:

The p superscript refers to predicted values. The variables are

r = position

v = velocity (dr/dt)

a = acceleration (d²r/dt²)

b = third derivative of position with respect to time (d³r/dt³)

The equations of motion are introduced by calculating the acceleration, a due to the force, F. The force in turn is calculated from the potential function V(r^p) at the new positions, r^p so that the correct acceleration is a^c = F/m = (-dV(r^p)/dr)/m. The predicted positions and velocities must be corrected. The correction term is proportional to the difference between the predicted and correct acceleration

The corrector step is

The positions, velocities etc. are now a better approximation to the true solutions of the equations of motion. The coefficients c_j are discussed in Appendix E of Allen and Tildesley. For the second order scheme shown above they are

c₀= 1/6

c₀= 5/6

c₀= 1

c₀= 1/3

A second order scheme means that the second order equations of motion are being used for the equations of motion (Lagrange equations of motion). There is a first order scheme in the event that the Hamiltonian equations of motion are to be used.

Accuracy of calculated average energy

Comparison of the Verlet algorithm and the Gear predictor-corrector have been carried out as a function of the time step, dt. The Gear-predictor method can be made more accurate by including more correction steps. The accuracy of both methods increases as dt².