1549380323-Statistical Mechanics Theory and Molecular Simulation

Integrating the equations of motion 99

prove the ergodicity or lack thereof in a system with many degrees of freedom. The
ergodic hypothesis tends to break down locally rather than globally.The virial theorem
tells us that the average energy in a given mode iskTat equilibrium if the system
has been able to equipartition the energy. Instantaneously, however, the energy of a
mode or degree of freedom fluctuates. Thus, if some particular mode has a high barrier
to surmount, a very long time will be needed for a fluctuation to occur that amasses
sufficient energy in this mode to promote barrier-crossing. Biological macromolecules
such as proteins and polypeptides exemplify this problem, as important conformations
are often separated by barriers in the space of the backbone dihedral angles or other
collective variables in the system. Many other types of systems have severe ergodicity
problems that render them challenging to treat via numerical simulation, and one
must always bear such problems in mind when applying numerical methods such as
molecular dynamics. Keeping such caveats in mind, we begin with a discussion of
numerical integrators.

3.8 Integrating the equations of motion: Finite difference methods

3.8.1 The Verlet algorithm

There are three principal aspects to a molecular dynamics calculation: 1) the model
describing the interparticle interactions; 2) the calculation of energies and forces from
the model, which should be done accurately and efficiently; 3) the algorithm used to
integrate the equations of motion. Each of these can strongly influence the quality of
the calculation and its ability to sample a sufficient number of microstates to obtain
reliable averages. We will start by considering the problem of devisinga numerical
integrator or solver for the equations of motion. Later in this chapter, we will consider
different types of models for physical systems. Technical aspects of force calculations
are provided in Appendix B.
By far the simplest way to obtain a numerical integration scheme is touse a Taylor
series. In this approach, the position of a particle at a timet+∆tis expressed in terms
of its position, velocity, and acceleration at timetaccording to:

ri(t+ ∆t)≈ri(t) + ∆tr ̇i(t) +

1

2

∆t^2 ̈ri(t), (3.8.1)

where all terms higher than second order in ∆thave been dropped. Sincer ̇i(t) =vi(t)
and ̈ri(t) =Fi(t)/miby Newton’s second law, eqn. (3.8.1) can be written as

ri(t+ ∆t)≈ri(t) + ∆tvi(t) +

∆t^2

2 mi

Fi(t). (3.8.2)

Note that the shorthand notation for the forceFi(t) is used in place of the full expres-
sion,Fi(r 1 (t),...,rN(t)). A velocity-independent scheme can be obtained by writing a
similar expansion forri(t−∆t):

ri(t−∆t) =ri(t)−∆tvi(t) +

∆t^2

2 mi

Fi(t). (3.8.3)

Adding eqns. (3.8.2) and (3.8.3), one obtains