1549380323-Statistical Mechanics Theory and Molecular Simulation

104 Microcanonical ensemble

X=σ

√

−2 lnξ 2 ′cos 2πξ 1

Y=σ

√

−2 lnξ 2 ′sin 2πξ 1. (3.8.20)

Thus, we obtain two Gaussian random numbersXandY. This algorithm for gener-
ating Gaussian random numbers is known asBox–Muller sampling. By applying the
algorithm to the Maxwell–Boltzmann distribution in eqn. (3.8.10), 3Ninitial velocity
vector components can be generated. Note, however, that if there are any constraints
on the system, the velocities must be projected back to the surface of constraint af-
ter the sampling is complete in order to ensure that the first time derivatives of the
constraint conditions are also satisfied. Moreover, for systems in∑ which the total force
N
i=1Fi= 0, the center-of-mass velocity

vcm=

∑N

∑i=1mivi

N

i=1mi

(3.8.21)

is a constant of the motion. Therefore, it is often useful to choose the initial velocities
in such a way thatvcm= 0 in order to avoid an overall drift of the system in space.
Once the initial conditions are specified, all information needed to start a simulation
is available, and an algorithm such as the Verlet or velocity Verlet algorithm can be
used to integrate the equations of motion.

3.9 Systems subject to holonomic constraints

In Section 1.9, we discussed the formulation of classical mechanics for a system sub-
ject to a set of holonomic constraints, that is, constraints which depend only on the
positions of the particles and possibly time:

σk(r 1 ,...,rN,t) = 0 k= 1,...,Nc. (3.9.1)

For the present discussion, we shall consider only time-independent constraints. In this
case, according to eqn. (1.9.11), the equations of motion can be expressed as

mi ̈ri=Fi+

∑Nc

k=1

λk∇iσk, (3.9.2)

whereλkis a set of Lagrange multipliers for enforcing the constraints. Although it is
possible to obtain an exact expression for the Lagrange multipliers using Gauss’s prin-
ciple of least constraint, the numerical integration of the equations of motion obtained
by substituting the exact expression forλkinto eqn. (3.9.2) would not exactly pre-
serve the constraint condition due to numerical errors, which would lead to unwanted
instabilities and artifacts in a simulation. In addition, Gauss’s equations of motion are
complicated non-Hamiltonian equations that cannot be treated using simple methods
such as the Verlet and velocity Verlet algorithms. These problems can be circumvented
by introducing a scheme for computing the multipliers “on the fly” in a simulation
in such a way that the constraint conditions are exactly satisfied within a particular
chosen numerical integration scheme. This is the approach we will now describe.