• What is molecular dynamics?
• Introduction tto a few integration schemes.

 

Motions of subunits (normally atoms) of one or more "molecules".

 

Two kinds of motions:

 
• Motions without disturbence to equilibrium state.

 
Example: Still water in a cup.
 
• Motions that changes macroscopic (global) state.

 
Example: Airflow.
 
 

Definition:

Molecular dynamics is a method in which the motion of each individual atom or molecule is computed according to Newton's second law.

Macroscopic quantities are extracted from the microscopic trajectories of particles.
 
 

Applications:

 
• The thermodynamic properties of gas, liquid, and solid.

• Phase transitions.

• Internal motions of macromolecules (e.g., polymers, DNA, proteins)

• Hydrodynamical fluid flow; plasma and electrons, transport phenomena, etc.

 
 

Equations of motions in molecular dynamics:

Newton's second law:
 
F = ma = m d²x/dt²
 
Its solution:
 
dx/dt = int F/m

x = int int F/m

 
which is an integration process.
 
For system of N atoms:
 
 
m_i d²x_i/dt² = sum _j F_ij
 
The integration is more complicated. Special computer numeric methods have been developed.
 
 
 

Solving differential equations is one of the major topics in numerical analysis. Several algorithms particularly suitable for molecular dynamics will be introduced.

 

Runge-Kutta method
 

This is a very straightforward but less accurate method.

 

Method
 

Our task is to solve the differential equation
 
 

dx/dt = f(t, y), x(t0)= x0
 
That is, given the function f(t,y), find the function x(t) such that it passes through the point (t0, x0), whose derivative is the given function. Any set of ordinary differential equations can be cast in this form if we allow x to be a vector of certain dimensions.

For example, a second-order differential equation can be rewritten as a set of two first-order differential equations. Thus it is suffice to just discuss how to solve the above equation, keeping in mind that x could mean
(x1, x2, ... , xd).

Clearly, the most obvious scheme to solve the above equation is to replace the differential by finite differences:
 
 

dt -> h

dx -> x(t+h) - x(t)
 
 

We then get the Euler method or first-order Runge-Kutta formula:
 
 
x(t+h) = x(t) + h f(t, y(t)) + O(h²)
 
The term first order refers to the fact that the equation is accurate to first order in the small step size h, thus the (local) truncation error is of order h². The Euler method is not recommended for practical use, because it is less accurate in comparison to other methods and it is not very stable.
 

Improvement of accuracy
 

The accuracy of the approximation can be improved by evaluating the function f at two points, once at the starting point, and once at the midpoint. This lead to the second-order Runge-Kutta or midpoint method.

 

k1 = h f(t, x(t))

k2 = h f(t+h/2, x(t)+k1/2)

x(t + h) = x(t) + k2 + O(h³)
 

However, the most popular Runge-Kutta formula is the following fourth-order one:
 
 
k1 = h f(t, x(t))

k2 = h f(t+h/2, x(t)+k1/2)

k3 = h f(t+h/2, x(t)+k2/2)

k4 = h f(t+h, x(t)+k3)

x(t + h) = x(t) + k1/6 + k2/3 + k3/3 + k4/6 + O(h⁵)
 

The fourth-order Runge-Kutta method is suitable for many applications for systems involving few degrees of freedom. But it was almost never used for molecular dynamics involving many degrees of freedom. Reasons for this:
 
• computationally demanding
 
four evaluations of the right-hand side of the equations
 
• lacks the time-reversal symmetry of the Newton equations.
 
 
 

Verlet algorithm
 

This algorithm is particularly suited for molecular dynamics. It has been widely used in many areas from simulation of liquids and solids to biological molecules.

Method

Our problem is to solve the set of Newton's equations:
 
 

m_i d² x_i / dt² = F_i
 
where x_i is the position vector of the i-th particle, m_i is the mass of this particle, and F_i is the total force acting on this particle. We'll discuss the form of the forces later on.

To get approximate difference formula for derivatives, let's look at the Taylor expansion of the function x(t):
 
 

x(t + h) = x(t) + h x'(t) + h²/2 x''(t) + h³/6 x'''(t)
                + O(h⁴)
 
Replacing h by -h, we get:
 
 
x(t - h) = x(t) - h x'(t) + h²/2 x''(t) - h³/6 x'''(t)
                 + O(h⁴)
 
By subtracting the 1st equation with the 2nd, we have:
 
 
x'(t) = [ x(t + h) - x(t - h) ]/2h + O(h²)
 
By adding the 1st equation with the 2nd, we have:
 
 
x''(t) = [ x(t + h) -2x(t) + x(t - h) ]/h² + O(h²)
 

The same formula applies if x is a vector. Replacing the second derivative by a difference using the above approximation, Newton's equation can be written as:
 
 

x_i(t + h) = - x_i(t - h) + 2x_i(t) + h2/m_i F_i(t)
                 + O(h⁴)
 
This equation was first used in 1791 by Delambre, and rediscovered many times, most recently by Verlet in 1960s for molecular dynamics. It is also called Stoermer algorithm.

It has the property that it is exactly time-reversible. That is, if you let the system develop according to this equation for some number of steps, you can back to your starting point exactly by reversing the time. This is the property of the original Newton's differential equation. If you were using the fourth-order Runge-Kutta, you do not have time-reversal symmetry exactly. This property is very important for long-time simulation.
 

The velocity of the particles can be computed from
 

v_i(t) ~ [ x_i(t + h) - x_i(t - h) ]/2h
 
 
Initial condition used in Verlet algorithm
 

One problem with Verlet algorithm is that how we get started, since one needs the values at two previous times.

Usually, the second-order differential equation will have a unique solution if the variables and their first derivatives are given. In our case, if the initial positions x_i(0) and velocities v_i(0) are given, our solution x(t) is uniquely determined. One way to start the simulation is to use:
 
 

x_i(h) = x_i(0) + h v_i(0) + h²/2m_i F_i(0)
 
The subsequent values x_i(2h), x_i(3h), etc, are obtained from the Verlet difference equation.
 
 
 

Symplectic algorithms

What is symplectic algorithm?
 

In the Verlet algorithm, our basic dynamic variables are the coordinates of the particles. And we have a somewhat inconvenient situation that the initial conditions are specified as the initial positions and initial velocities, while the difference scheme involves position only. As you may have learned from mechanics, it is possible to reformulate the Newtonian mechanics so that both positions and momenta are the basic dynamic variables. This formulation is called Hamiltonian dynamics (which is equivalent to Newtonian dynamics). A Hamiltonian system is completely described by the Hamiltonian function (kinetic energy plus potential energy):
 
 

H(p,q) = K + V
 
together with Hamilton's equations of motion:
 
 
dp_i/dt = - dH/dq_i,    dq_i/dt = dH/dp_i
 
 
where H is considered to be a function of the generalized coordinates qi and generalized momentum p_i. The derivative on H is a partial derivative. The index i labels the degree of freedom, not the particle.
 

Algorithms for solving this set of  Hamilton's equations of motion are symplectic algorithms.
 

Systems for which a symplectic algorithm can be derived:
 

The symplectic algorithm can be derived if the Hamiltonian takes a simple form:
 
 

H = K + V = sum_i p_i²/2m_i + V(q₁, q₂, ...)
 
This is certainly the case if we use Cartesian coordinates and the interactions between particles depend only on the positions of the particles.

With this particular form of Hamiltonian, Hamilton's equations can be written as:
 
 

dp_i/dt = - dV/dq_i = F_i

dq_i/dt = dK/dp_i = p_i/m_i
 
 
This set of equationa reduce to Newton's equations if one eliminates the variable p_i.
 
 
 

Symplectic algorithms:
 
 

Algorithm A
 
 

We use an Euler type approximation for the momentum derivatives, and then a similar Euler algorithm for the position derivatives. However, for the positions, we use immediately the new result for the momenta. Thus, it is not exactly an Euler algorithm. For brevity, we'll use a vector notation:
 
 
p = ( p1, p2, ... )

q = ( q1, q2, ... )

F = ( F1, F2, ... )
 

and assume that all particles have the same mass m. With these notations and assumption, we have:
 
 

p(t + h) = p(t) + hF(q(t))

q(t + h) = q(t) + h/m p(t + h)
 

Clearly, the algorithm is accurate to first order in h.
 
 

Algorithm B
 
 

We do exactly the same thing as in algorithm A, except that we reverse the order of calculation. The new position is found first.
 
 
q(t + h) = q(t) + h/m p(t)

p(t + h) = p(t) + hF(q(t + h))
 

If you compare algorithms A and B, you will find that it is not exactly the same. It is very interesting to note that both algorithm A and B are in fact mathematically equivalent to the Verlet algorithm if you eliminate the momenta from the equations. However, they have better roundoff error control over the plain Verlet algorithm.
 
 

Algorithm C (second-order symplectic algorithm or velocity form of the Verlet algorithm):
 
 

This is the most widely used algorithm that we are really interested.

This algorithm is derived by dividing the step size h into two, for the first half step size h/2, we use the algorithm A, for the second half of the step, we use the algorithm B. In total, we have full step size h. That is, use the formula for algorithm A but set h/2. We get equations relating the quantities at time t to time t + h/2. Then we use the equations for algorithms B with step size again h/2, but the time is for t+h/2 to t + h. We get four equations:
 
 

p(t + h/2) = p(t) + h/2 F(q(t))

q(t + h/2) = q(t) + h/2m p(t + h/2)

q(t + h) = q(t + h/2) + h/2m p(t + h/2)

p(t + h) = p(t + h/2) + h/2 F(q(t + h))
 

Now we eliminate the intermediate quantities evaluated at time t + h/2, we get the final result for the symplectic algorithm:
 
 

q(t + h) = q(t) + h/m p(t) + h2/2m F(q(t))

p(t + h) = p(t) + h/2 [ F(q(t)) + F(q(t + h)) ]
 
 

Note that the coordinates q have to be evaluated first since the new values at time t+h are used to evaluate the forces.

The algorithm C is also second-order in step size h. From the derivation it is not obvious that it is necessarily superior than Verlet algorithm. However, the symplectic algorithm has several advantages:
 

• It is time-reversible
 
• The system can be started naturally, with initial position q0 and initial momentum p0 = mv0
 
• The symplectic algorithms A, B, and C have one important property that share with the original Hamiltonian system---they preserve Poincare invariants.
 
 
This last property of Hamiltonian system is usually discussed in advanced mechanics or dynamical system courses. It is about the phase space (p,q) flow property. A solution with initial condition (p0,q0) of the Hamiltonian system p=p(t,p0,q0), q = q(t,p0, q0) can be viewed as a map in the phase space from the point (p0,q0) to the point (p,q), where the time parameter t is taken as some constant. The map is also called a canonical transformation, or symplectic transformation.

One of the series of Poincare invariants is the phase space volume. Certain region of points of volume V0 can be mapped to a set of new points forming volume V. Invariance means V=V0. This is called Liouville theorem.

Not only is the volume in phase space invariant under the mapping of the dynamics, but also there are other lower dimensional objects (hypersurfaces) which are invariant. The symplectic algorithms preserve this invariant property exactly.

 
 

Homework