Integrating the equations of motion 101is a fundamental symmetry of Hamilton’s equations that should be preserved by a nu-
merical integrator. The second is the symplectic property of eqn.(1.6.29); we will dis-
cuss the importance of the symplectic property for numerical stability in Section 3.13.
While there are classes of integrators that purport to be more accurate than
the simple second-order Verlet and velocity Verlet algorithms, for example predictor-
corrector methods, we note here that many of these methods are neither symplectic
nor time-reversible and, therefore, lead to significant drifts in thetotal energy when
used. In choosing a numerical integration method, one shouldalwaysexamine the
properties of the integrator and verify its suitability for a given problem.
3.8.3 Choosing the initial conditions
At this point, it is worth saying a few words about how the initial conditions for a
molecular dynamics calculation are chosen. Indeed, setting up an initial condition can,
depending on the complexity of the system, be a nontrivial problem.For a simple
liquid, one might start with initial coordinates corresponding to the solid phase of the
substance and then simply melt the solid structure under thermodynamic conditions
appropriate to the liquid. Alternatively, one can begin with random initial coordinates,
requiring only that the distance between particles be large enough to avoid strong re-
pulsive forces initially. For a molecular liquid, initial bond lengths and bend angles may
be dictated by holonomic constraints or may simply be chosen to be equilibrium values.
For more complex systems such as molecular crystals or biological macromolecules, it
is usually necessary to obtain initial coordinates from an experimental X-ray crystal
structure. Many such crystal structures are deposited into structure databases such
as the Cambridge Structure Database, the Inorganic Crystal Structure Database, or
the Protein Data Bank. When using experimental structures, it might be necessary to
supply missing information, such as the coordinates of hydrogen atoms, that cannot
be experimentally resolved. For biological systems, it is often necessary to solvate the
macromolecule in a bath of water molecules. For this purpose, one might take coordi-
nates from a large, well-equilibrated pure water simulation, place themacromolecule
into the water bath, and then remove waters that are within a certain distance (e.g. 1.8
̊A) of any atom in the macromolecule, being careful to retain crystallographic waters
bound within the molecule. After such a procedure, it is necessary to re-establish equi-
librium, which typically involves adjusting the energy to give a certain temperature
and the volume (Chapter 4) to give a certain pressure (Chapter 5).
Once initial coordinates are specified, it remains to set the initial velocities. This is
generally done by “sampling” the velocities from a Maxwell–Boltzmann distribution,
taking care to ensure that the sampled velocities are consistent with any constraints
imposed on the system. We will treat the problem of sampling a distribution more
generally in Chapter 7; however, here we provide a simple algorithm for obtaining an
initial set of velocities.
The Maxwell–Boltzmann distribution for the velocityvof a particle of massmat
temperatureTis
f(v) =( m
2 πkT) 1 / 2
e−mv(^2) / 2 kT
. (3.8.10)
The distributionf(v) is an example of aGaussian probability distribution. More gen-