by the intramolecular kinetic energy. Of the total of six degrees of freedom

three are internal: the bond vibration and two rotational degrees of free-

dom, together good for an average of 3 kB T kinetic energy. The harmonic

2

bond vibrational degree of freedom has an average kinetic energy of 1 kB T ,

2

which is equal to the average potential energy 2 k(d ’ d0 )

1 2 . Therefore the

contribution to the virial is 1 F d = 1 k d(d ’ d0 ) = 1 k (d ’ d0 )2 = 1 kB T ,

2 2 2 2

which exactly cancels the average kinetic energy of the bond vibration. But

how about the contribution of the two rotational degrees of freedom? They

have an average kinetic energy of kB T , but where is the virial compensation?

The answer is that rotation involves centrifugal forces on the bond, which is

then stretched, causing a compensating elastic force in the bond direction.

That force causes a contribution to the intramolecular virial that exactly

compensates the rotational kinetic energy (see exercise 17.8).

492 Review of statistical mechanics

With that problem solved, the next question is what happens if the bond

length is treated as a constraint? In that case there is no kinetic energy in

the bond vibration and their is no contribution to the virial due to vibra-

tional motion. But there still is rotational kinetic energy and that is still

compensated by the contribution to the virial of the constraint force. So

when constraints are used in an atomic description, the constraint forces

must be computed and accounted for in the atomic virial. Constraint forces

are the forces that compensate components of other forces in the constraint

directions; they follow directly from the constraint computation (see Section

15.8 on page 417).

17.8 Liouville equations in phase space

A classical system of particles that evolves under Hamiltonian equations

of motion (see Section 15.3 on page 399) follows a trajectory in the 2n-

dimensional space of all its (generalized) coordinates qi , i = 1, . . . n and con-

jugate momenta pi i = 1, . . . n, called phase space. Here n is the number of

degrees of freedom, which equals 3— the number of particles minus the num-

ber of constrained degrees of freedom, if any. A particular con¬guration of

all coordinates and momenta at time t de¬nes a point in phase space. Often

it is convenient to denote the 2n-dimensional vector (q1 , . . . qn , p1 , . . . pn )T

by the vector z (z ∈ R2n ).28 We consider for the time being systems that

obey a ¬rst-order di¬erential equation, so that the evolution of a point in

phase space is a deterministic initial value problem.

Consider evolution under the Hamilton equations (15.9):

‚H(q, p)

qi =

™ , (17.148)

‚pi

‚H(q, p)

pi = ’

™ . (17.149)

‚qi

In symplectic notation, and introducing the matrix

01

L0 = , (17.150)

’1 0

where 0 is a n — n all-zero matrix, and 1 is a n — n diagonal unit matrix,

28 This notation is referred to as the symplectic notation. A symplectic mapping is an area (and

volume) conserving mapping, such as the transformation of z for systems obeying Hamilton™s

equations. A symplectic algorithm is an algorithm to solve the time evolution that conserves

area and volume in phase space.

17.8 Liouville equations in phase space 493

Hamilton™s equations can be written as

‚H

™

z = L0 . (17.151)

‚z

Another notation is writing the operation on z on the right-hand side of

(17.151) as an operator :

ˆ

z = iLz.

™ (17.152)

ˆ

L is the Liouville operator, which we will de¬ne below in the context of the

rate of change of density in phase space.

The rate equation (17.152) is more general than its Hamiltonian form of

(17.151), and may for example contain the e¬ects of external time-dependent

¬elds. The operator is then called the generalized Liouville operator.

The rate equation can be formally integrated to yield a time evolution

equation for the point z in phase space:

t

ˆ

iL(„ ) d„ ] z(0),

z(t) = exp[ (17.153)

0

which for time-independent operators reduces to

z(t) = eiLt z(0).

ˆ

(17.154)

These equations are formally elegant but do not help solving the equations

of motion. Exponential operators and the integral over a time-dependent

operator must all be carefully de¬ned to be meaningful (see Section 14.6 on

exponential operators on page 385).

In statistical mechanics we deal not with single trajectories, but with

distribution functions f (z, t) in phase space. We are interested in the time

evolution of such distribution functions and in the equilibrium distributions

corresponding to the various ensembles. When we know the distribution

function, we can determine the observable ensemble average of any variable

that is known as a function of the point in phase space:

A (t) = A(z) f (z, t) dz. (17.155)

In order to ¬nd a di¬erential equation for the time evolution of f (z, t) we

try to ¬nd the time derivative ‚f /‚t. Since f concerns a probability distri-

bution, its integrated density must be conserved, and any change in density

™

must result from in- or outgoing ¬‚ux J = f z. The conservation law in 2n

dimensions is similar to the continuity equation in ¬‚uid dynamics (see (9.3)

494 Review of statistical mechanics

on page 282):

‚f (z, t)

= ’∇ · J = ’∇ · (f z) = ’z · ∇f ’ f ∇ · z,

™ ™ ™ (17.156)

‚t

where ∇ stands for the gradient operator in z: (‚/‚z1 , . . . , ‚/‚z2n ). Writing

the time derivative as the material or Lagrangian derivative D/Dt, i.e., the

derivative seen when traveling with the ¬‚ow (see page 282), we see that

Df def ‚f

+ z · ∇f = ’f ∇ · z.

™ ™

= (17.157)

Dt ‚t

This equation is often referred to as the generalized Liouville equation. For a

Hamiltonian system the term on the right-hand side is zero, as we see using

(17.151):

2n

‚2H

‚H

∇ · z = ∇ · L0

™ = L0ij

‚z ‚zi ‚zj

i,j=1

n

‚2H ‚2H

’

= = 0. (17.158)

‚qi ‚pj ‚pi ‚qj

i=1