Numerical techniques: Molecular Dynamics
I. The idea of molecular dynamics
Molecular dynamics (MD) is a method for computing equilibrium and kinetic properties of
molecular systems, which obey the laws of classical physics. In general, any molecular process
with the characteristic timescale > 0 ≈0.2 ps can be described using classical physics.
Therefore, except for the fluctuations in bond lengths and angles all other motions in molecular
systems can be treated as classical. As a consequence the motions of atoms can be described by
Newton equations of motions. The main advantage of MD is its ability to produce “true”
microscopic dynamics governed by the underlying energy landscape and interatomic forces. The
implementation of MD is similar to actual experiment and involves the following stages:
1. Setting the parameters describing the conditions of MD simulations, such as temperature,
number of atoms etc.
2. Initialization, which includes reading in the coordinates of atoms and generation of initial
distribution of velocities
3. Computation of forces
4. Integration of Newton’s equation of motion
5. Repeating steps 3 and 4 until MD simulations are done for the desired timescale
6. Computation of averages
The most typical mistakes in MD are associated with inadequate sampling of molecular
conformations, unphysical choice of initial conditions, and wrong choice of variables describing
the system. In what follows we consider the stages in MD simulations.
II. Initial conformation
After setting the values of parameters for MD simulations the initial structure is read in. It is
important that this conformation does not include atom overlaps, unusual local conformations
(i.e., unphysical values of dihedral angles), which may result in large forces causing the
instability of the integrator. Another consideration to keep in mind is that the initial structure must
be consistent with the purpose of simulations. For example, if the goal is to study the kinetics of
relaxation to the native state from unfolded conformations, the initial structure must resemble
random coil (unfolded) states. If the purpose is to probe equilibrium properties near the native
state, then the initial structure must be native-like. In this case, the PDB structure may serve as
starting conformation.
III. Initial distribution of velocities
The initial distribution of velocities can be drawn from the Maxwell-Boltzmann distribution
1
mvi2,
m 2 2 k BT
P(vi , )
,
2k T e
B
(1)
1
where vi,a is the (=x,y,z) component of the velocity of the atom i, m is a atom mass, T is a
temperature, and kB is a Boltzmann constant. The distribution in Eq. (1) may also be used to
define the instantaneous temperature T(t) using the equipartion theorem as
2
mv
2
1
k BT ,
2
(2)
which relates the average kinetic energy with the temperature (<…> brackets indicate the
ensemble average). Eq. (2) can be directly obtained from Eq. (1). Because the ensemble average
corresponds to the average over all velocities of atoms, the instantaneous temperature T(t) is
defined as
k B T (t )
1
Nf
mv ,
2
i,
(3)
i,
where Nf is the number of degrees of freedom. Therefore, Eq. (3) allows us to compute the
instantaneous temperature from the distribution of velocities. It is also clear that for a given
realization of the velocity distribution, T(t) is not strictly equal to T. Indeed, although the
velocities are generated using the distribution in Eq. (1) at the temperature T, the molecular
system in MD simulations contains only finite number of atoms and the actual (instantaneous)
temperature T(t) will deviate from T. To keep the temperature constant, one can rescale the
velocities according to
vi’ ,
T
vi , .
T (t )
(4)
It is straightforward to show that the instantaneous temperature after rescaling T’(t)≡T. If
rescaling is not carried out, the relative fluctuations of temperature in the system of N atoms is
given by
T (t ) T
T (t )
2
(t ) T (t )
T (t )
1
2 2
~N
1
2
.
(5)
Thus, for the system of about 1000 atoms the fluctuations in T(t) are about 3%.
IV. Force computation
Because the forces in Newton equations must be computed between all pairs of atoms, their
computation is numerically expensive. As a result, the number of force evaluations scales with
the number of atoms N as N2. (There exist efficient numerical methods, which reduce this scaling
to almost linear. These methods will be considered in Lecture 12.) To illustrate the computation
of forces consider a typical interatomic Lennard-Jones potential
2
12
6
V LJ (rij ) h 12 2 6
r
rij
ij
.
(6)
Here rij ri r j is the distance between atoms i and j ( ri , ri are their position vectors,
respectively), the equilibrium distance between atoms is , and the energy of attractive
interatomic interactions is h. The potential in Eq. (6) reaches minimum at rij 21 / 6 . Because
the force acting on the atom i is the negative gradient of the interatomic potential energy, we
write
f i, x
xi x j
V
V rij
LJ LJ
12 h
x i
rij x i
rij2
rij
12
6
rij
(7)
The force fj,x has the opposite sign to fi,x. The y- and z-components of forces are computed in a
similar way.
V. Integration of Newton equations of motion
In this section the Verlet algorithm is derived. Using Taylor expansion write the position r(t+t)
and r(t-t) at times t+t and t-t as
r (t t ) r (t ) v(t )t
f (t ) 2 r ‘ ‘ ‘ (t ) 3
t
t o(t 4 )
2m
6
(8)
r (t t ) r (t ) v(t )t
f (t ) 2 r ‘ ‘ ‘ (t ) 3
t
t o(t 4 ) .
2m
6
(9)
In Eqs. (8,9) we omit the atom or component indexes, use Newton equation f (t ) mr (t ) and
denote v(t ) r (t ) . Adding Eqs. (8) and (9) yields
r (t t ) 2r (t ) r (t t )
f (t ) 2
t o(t 4 ) .
m
(10)
In order to propagate equations of motion the positions at t and t-t as well as the force at t must
be known. The accuracy of the algorithm is o(t4). The velocities, which are needed to compute
the kinetic energy and temperature, can be calculated using central finite divided difference
scheme as
v(t )
r (t t ) r (t t )
o(t 2 ) .
2t
(11)
Note that the accuracy of computations of velocities is lower than that of coordinates. The
implementation of Verlet algorithm is given in Box 2. Because we take into account only the
3
interactions and forces between the atoms in the system, the total energy and momentum must be
conserved. During MD simulations these two conditions must be explicitly checked.
VI. Requirements for MD integrator
There are several requirements that all MD algorithms must meet.
1. MD integration algorithms must be fast, although this requirement is not crucial, because
the time spent on propagation of equations of motion is typically far less than the time
needed to complete full evaluation of energy and forces.
2. The MD algorithm must be accurate. High accuracy would allow us to use larger
integration steps and thus reduce the amount of expensive evaluations of energy and
forces. However, any MD algorithm is inherently approximate and, therefore, with time
the deviation between the computed trajectory and actual, true trajectory always
develops. This issue is related to so called Lyapunov instabilities. Consider two, almost
identical initial conditions (1) and (2), which include coordinates and momenta. As both
molecular systems evolve with time, the difference between current conformations
r (t ) ri ,1 (t ) ri , 2 (t ) would grow exponentially as r (t ) ~ e t , where >0 is a
2
i
Lyapunov exponent and i is the atom index (for more information, see J. Comp.Phys.
151, 9 (1999)). Treating the inaccuracy developing in MD trajectories as a Lyapunov
instability, we conclude that the error in MD trajectory increases exponentially. However,
it is believed that MD trajectories after the onset of Lyapunov instability follow closely
“shadow” orbits (trajectories). The “shadow” orbits represent some “other” true MD
trajectories, which pass through slightly perturbed initial conditions as compared to those
that were actually used in simulations.
3. Because Newton equations are time reversible, so must be the MD algorithms. This
implies that the algorithms should not be changed upon the substitution t →-t.
4. MD algorithms must conserve energy and momentum. It is important to note that longterm energy drift develops in any MD algorithm. Therefore, the algorithm selection
should be also based on minimizing energy drift on a long time scale.
VII. Additional MD integrators
There are several other algorithms that may be used in MD. The Euler algorithm is based on
direct application of Taylor expansion, which implies using Eq. (8) with the terms up to t2 order.
The Euler algorithm is not time reversible and does not conserve energy. Similarly, all RungeKutta algorithms are not energy preserving and not time reversible and consequently cannot be
used for MD applications. It is also possible to derive a version of Verlet algorithm, which
explicitly propagates coordinates and velocities (a Leap-frog algorithm). Using Verlet algorithm
we write the velocities at the half integration steps
r (t t ) r (t )
t
)
2
t
r (t ) r (t t )
t
v(t )
.
2
t
v(t
(12)
(13)
4
From Eq. (12) we get
r (t t ) r (t ) v(t
t
)t
2
(14)
Using Verlet algorithm for coordinates in Eq. (12) we write
v(t
t
)
2
f (t ) 2
t r (t )
f (t )
t
m
v(t )
t
t
2
m
2r (t ) r (t t )
(15)
Therefore, to propagate equations of motion using Leap-frog algorithm one need f(t), r(t) and v(tt/2). The inherent problem of Leap-frog algorithm is that coordinates and velocities are not
computed synchronously. As a result it is difficult to compute the total energy of the system.
Leap-frog algorithm generates trajectories identical to Verlet algorithm.
The velocity form of Verlet algorithm (VFVA) computes the coordinates and velocities at the
same time. According to VFVA the position r(t+t is computed using Euler algorithm (Eq. (8)),
but the velocity is obtained from forces at t and t+t as
r (t t ) r (t ) v(t )t
v(t t ) v(t )
f (t ) 2
t
2m
(16)
f (t ) f (t t )
t .
2m
(17)
It can be shown that VFVA is equivalent to Verlet algorithm. Using Taylor expansion we get
r (t 2t ) r (t t ) v(t t )t
r (t ) r (t t ) v(t )t
f (t t ) 2
t
2m
(18)
f (t ) 2
t
2m
(19)
Adding Eqs. (18) and (19) we obtain
r (t 2t ) r (t ) 2r (t t ) (v(t t ) v(t ))t
f (t t ) f (t ) 2
t
2m
or taking into account Eq. (17)
r (t 2t ) 2r (t t ) r (t )
f (t t ) 2
t ,
m
which is the expression for Verlet algorithm. The use of VFVA offers important advantages over
standard form of Verlet algorithm. First, velocities and coordinates are propagated synchronously.
5
Second, computation of velocities and coordinates at t+t requires storing the positions and
velocities only at t, whereas in the original form of Verlet algorithm the positions at t-t must be
retained as well. The peculiarity of VFVA is that coordinates must be advanced first (Eq. (16)),
from which the forces at t+t are computed. After this the computation of velocities takes place.
Verlet algorithms fare very well on all points in section VI except for perhaps (2). Although they
do not show very high short-term accuracy, the long-term energy drift is small as compared to
most other algorithms. More importantly, Verlet algorithms are time reversible and energy
preserving. Almost all available molecular dynamics codes use, at least as an option, the velocity
form of Verlet algorithm.
Box 1 Computation of forces
subroutine Force
en = 0
do i=1,npart
f(i) = 0.0
enddo
do i=1,npart-1
do j=i+1,npart
xr = x(i)-x(j)
yr = y(i)-y(j)
zr = z(i)-z(j)
r2=xr**2+yr**2+zr**2
if(r2<rc2)then
r2i=1/r2
r6i=r2i**3
ff=12.0*r2i*r6i*(r6i-1.0)
f(i) = f(i)+ff*xr
f(j) = f(j) –ff*xr
en = en + r6i*(r6i-2.0)
endif
enddo
enddo
return
end subroutine Force
# initialization of energy and
# force array
# loop over all (i,j) pairs
#
# x,y,z distances between particles i and j
# applying the cut-off for distance r2
# computing force factor
# adding x-force acting on the particle i
# adding x-force acting on the particle j
# adding energy
Notes:
1. h and are set to unity
2. npart is the number of atoms
3. the algorithm computes only x-components of forces, but full Lennard-Jones energy
Box 2 Integration of equations of motions using Verlet algorithm
subroutine VerletIntegrator
sumv = 0
# initializing the total momentum
6
sumv2 = 0
do i=1,npart
xx=2*x(i)-xm(i)+f(i)*delta**2
vi=(xx-xm(i))/(2*delta)
sumv = sumv + vi
sumv2 = sum2 + vi**2
xm(i) = x(i)
x(i) = xx
enddo
temp = sum2/npart
etot = en + sum2/2
return
end subroutine VerletIntegrator
# and kinetic energy
# implementation of Verlet algorithm
# getting the velocity of i
# adding up momenta
# and energies
# updating coordinates
#
# computing instantaneous temperature
# and total energy (en is potential energy)
Notes:
1. mass and Boltzmann constant are set to unity
2. delta is the integration time step
3. for simplicity algorithm is presented only for x-coordinate
Homework:
1. Read online notes
2. Find the error in section IV.
7
