8.2 Molecular dynamics at constant energy 201
conserved too, so the time averages of physical quantities obtained by this type of
simulation are equivalent to averages in the microcanonical or (NVE) ensemble. In
this section we describe the microcanonical MD method in more detail.
The algorithm of a standard MD simulation for studying systems in equilibrium
is the following:
- Initialise;
- Start simulation and let the system reach equilibrium;
- Continue simulation and store results.
We will now describe these main steps in more detail.
Initialise:The number of particles and the form of the interaction are specified.
The temperature is usually of greater interest than the total energy of the system
and is therefore usually specified as an input parameter. We shall see below how
the system can be pushed toward the desired temperature.
The particles are assigned positions and momenta. If a Lennard–Jones potential is
used, the positions are usually chosen as the sites of a Bravais-fcc lattice, which is the
ground state configuration of the noble gases like argon (although the Lennard–Jones
system is hexagonal close-packed in the ground state [6]). The fcc lattice contains
four particles per unit cell, and for a cubic volume the system contains therefore 4M^3
particles,M=1, 2,...This is the reason why MD simulations with Lennard–Jones
interactions are often carried out with particle numbers 108, 256, 500, 864, ....
The velocities are drawn from a Maxwell distribution with the specified temper-
ature. This is done by drawing thex,yandzvelocity components for each particle
from a Gaussian distribution; for thex-component of the velocity this distribution
is exp[−mv^2 x/( 2 kBT)]. InAppendix B3it is described how random numbers with
a Gaussian distribution can be generated. After generating the momenta, the total
momentum is made equal to zero by calculating the average momentump ̄per
particle, and then subtracting an amountp ̄from all the individual momentapi.
Start simulation and let the system reach equilibrium:The particles being
released from fcc lattice positions, the system is generally not in equilibrium and
during the initial phase of the simulation it is given the opportunity to relax. We now
describe how the integration of the equations of motion is carried out and how the
forces are evaluated. Finally we shall explain how in this initial phase the desired
temperature is arrived at.
Numerical algorithms for molecular dynamics will be considered in detail in
Section 8.4. Suffice it here to mention briefly the most widely used algorithm
which is simple and reliable at the same time – the Verlet algorithm (see also
Appendix A7.1). The standard form of the Verlet algorithm for the integration of
the equation of motion of a single particle subject to a forceFdepending only on
the position of the particle reads
r(t+h)= 2 r(t)−r(t−h)+h^2 F[r(t)]/m (8.4)
