In an isotropic ¬‚uid all six o¬-diagonal elements have the same correlation

function and one can best use the average of the correlation functions of all

o¬-diagonal elements of σ.

The determination of viscosity through the Green-Kubo relation requires

accurate determination of the (integral) of the correlation function of a heav-

ily ¬‚uctuating quantity. Hess (2002a,b) concluded that it is more e¬cient

to use non-equilibrium molecular dynamics (NEMD) to determine viscosity

coe¬cients (see Section 18.5).

18.4 The Einstein relation

All Green-Kubo relations contain the integral of an autocorrelation function

of a ¬‚uctuating observable f (t):

∞

f (0)f („ ) d„. (18.53)

0

Numerical evaluations of such integrals are di¬cult if no knowledge on the

analytical form of the tail of the correlation function is available. The statis-

tics on the tail and “ consequently “ on the integral, is often poor.

An alternative is to monitor not f (t), but its integral F (t):

t

F (t) = f (t ) dt , (18.54)

0

18.5 Non-equilibrium molecular dynamics 519

and observe the behavior of F 2 (t) for large t. The following Einstein

relation is valid:

∞

d2

lim F (t) = 2 f (0)f („ ) d„. (18.55)

t’∞ dt 0

This means that F 2 (t) , plotted versus t, should approach a straight line.

A common application is the determination of the single-particle di¬usion

constant

∞

D= v(0)v(„ ) d„, (18.56)

0

by observing the mean-squared displacement x2 (t) , which approaches 2Dt

for times much longer than the correlation time of the velocity.

Proof We prove (18.55). Consider

t

d2

F (t) = 2 F (t)f (t) = 2 f (t )f (t) dt .

dt 0

By substituting „ = t ’ t the right-hand side rewrites to

t

f (t ’ „ )f (t) d„.

2

0

Since the ensemble average does not depend on the time origin, the integrand

is a function of „ only and is equal to the autocorrelation function of f (t).

Obviously, the limit for t ’ ∞ yields (18.55).

18.5 Non-equilibrium molecular dynamics

When (small) external “forces” are arti¬cially exerted on the particles in a

molecular-dynamics simulation, the system is brought (slightly) out of equi-

librium. Following an initial relaxation the system will reach a steady state

in which the response to the disturbance can be measured. In such non-

equilibrium molecular dynamics (NEMD) methods the “forces” are chosen

to represent the gradient that is appropriate for the transport property of

interest. When the system under study is periodic, it is consistent to apply

a gradient with the same periodicity. This implies that only spatial Fourier

components at wave vectors which are integer multiples of the reciprocal ba-

sic vectors (2π/lx , 2π/ly , 2π/lz ) can be applied. This limitation, of course, is

a consequence of periodicity and the long-wavelength limit must be obtained

by extrapolation of the observed box-size dependence. We now consider a

few examples (Berendsen, 1991b).

520 Linear response theory

18.5.1 Viscosity

In order to measure viscosity, we wish to impose a sinusoidal shear rate over

the system, i.e., we wish to exert an acceleration on each particle. This is

accomplished by adding at every time step ”t a velocity increment ”vx to

every particle

”vix = A”t cos kyi . (18.57)

Here, A is a (small) amplitude of the acceleration and k = 2π/ly is the

smallest wave vector ¬tting in the box. Any multiple of the smallest wave

vector can also be used. When m is the mass of each particle and ρ is the

number density, the external force per unit volume will be:

ext

fx (y) = mρA cos ky. (18.58)

Since there is no pressure gradient in the x-direction, the system will react

according to the Navier“Stokes equation:

‚ 2 ux

‚ux

mρ =· + mρA cos ky. (18.59)

‚y 2

‚t

The steady-state solution, which is approached exponentially with a time

constant equal to mρ/·k2 (Hess, 2002b), is

mρ

ux (y) = A cos ky. (18.60)

·k2

Thus, the viscosity coe¬cient · is found by monitoring the gradient in

the y-direction of the velocities vx of the particles. Velocity gradients in

non-equilibrium molecular dynamics simulations are determined by a least-

squares ¬t of the gradient to the particle velocities. A periodic gradient can

be measured by Fourier analysis of the velocity distribution.

18.5.2 Di¬usion

Self-di¬usion coe¬cients can easily be measured from an equilibrium sim-

ulation by monitoring the mean-square displacement of the particle as a

function of time and applying the Einstein relation. The di¬usion coe¬-

cient measured this way corresponds to the special case of a tracer di¬usion

coe¬cient of a single tagged particle that has the same interactions as the

other particles. In general, the tagged particles can be of a di¬erent type

and can occur in any mole fraction in the system. If not dilute, the di¬us-

ion ¬‚ux is in¬‚uenced by the hydrodynamic interaction between the moving

particles.

18.5 Non-equilibrium molecular dynamics 521

Consider a binary mixture of two particle types 1 and 2, with mole frac-

tions x1 = x and x2 = 1 ’ x. Assume that the mixture behaves ideally:

μi = μ0 + RT ln xi . (18.61)

i

Now we can derive the following equation:

D

u1 ’ u2 = ’ ∇x. (18.62)

x(1 ’ x)

In an NEMD simulation we apply two accelerations a1 and a2 to each of

the particles of type 1 and 2, respectively. This is done by increasing the

velocities (in a given direction) every step by a1 ”t for species 1 and by a2 ”t

for species 2. The total force on the system must be kept zero in order to

avoid acceleration of the center of mass:

M1 xa1 + M2 (1 ’ x)a2 = 0, (18.63)

where M1 and M2 are the molar masses of species 1 and 2, respectively. The

balance between driving force and frictional force is reached when

RT

M2 a2 = ’ x(u1 ’ u2 ), (18.64)

D

or, equivalently,

RT

(1 ’ x)(u1 ’ u2 ).

M1 a1 = (18.65)

D