e’iω„ d„.

µr (ω) = 1 ’ P (0) · P („ ) (18.40)

3kB T d„

0

At zero frequency the static dielectric constant is obtained:

V

P2 .

µr (0) = 1 + (18.41)

3kB T

In simulations one monitors the total dipole moment M = i μi = V P and

computes (1/V ) M (0) · M (t) . This is generally not a trivial calculation

because M is a single quantity that ¬‚uctuates slowly and long simulations

are needed to obtain accurate converged values for the correlation function.

Matters are somewhat more complicated than sketched above.2 In fact,

(18.40) and (18.41) can only be trusted for very dilute systems in which

mutual dipole interactions can be ignored. The reason is that the local

electric ¬eld, to which the molecules respond, includes the ¬elds due to other

dipoles and the reaction ¬eld introduced by boundary conditions. Without

derivation we give the correct result (Neumann and Steinhauser, 1983) for

the relation between the frequency-dependent dielectric constant and the

correlation function of the total dipole moment M (t) for the case that a

reaction ¬eld is employed in the simulation (see Section 6.3.5 on page 164

and (13.82) on page 347 for a description of reaction ¬elds):

∞

2µRF + 1 1 d

e’iωt dt ,

[µr (ω) ’ 1] ’ M (0) · M (t)

=

2µRF + µr (ω) 3µ0 V kB T dt

0

(18.42)

2 The theory relating dipole ¬‚uctuations with dielectric constants goes back to Kirkwood (1939).

The theory for deriving the dielectric constant from dipole ¬‚uctuations in simulations with

various boundary conditions is most clearly given by Neumann (1983), with extension to the

frequency-dependent case by Neumann and Steinhauser (1983). The theory was tested on a

Stockmayer ¬‚uid (Lennard“Jones particles with dipole moments) by Neumann et al. (1984).

516 Linear response theory

’ f ” x” z

vx(y)

y

”z x

f ” x” z

z

”x

Figure 18.2 Two planes, moving relative to each other in a ¬‚uid, experience a

viscous drag force proportional to the velocity gradient.

where µRF is the relative dielectric constant used for the reaction ¬eld. For

the static dielectric constant µs = µr (0) it follows that

2µRF + 1 1

(µs ’ 1) M2 .

= (18.43)

2µRF + µr (ω) 3µ0 V kB T

Equations (18.42) and (18.43) are implicit equations for µr ; they reduce to

the simpler (18.40) and (18.41) when µRF = ∞, i.e., for conducting boundary

conditions. These are also valid for the use of complete lattice sums with

dipole correction (see the discussion of tin-foil boundary conditions on page

373). For simulations with µRF = 1, i.e., using a cuto¬ radius for Coulomb

interactions, the dipole ¬‚uctuations are quenched and become rather in-

sensitive to the the value of the dielectric constant. Such simulations are

therefore unsuitable to derive dielectric properties.

18.3.2 Viscosity

Consider an isotropic ¬‚uid between two plates (each in the xz-plane and

separated in the y-direction) that move with respect to each other in the x-

direction (Fig. 18.2), causing a laminar ¬‚ow with velocity gradient dvx /dy.

18.3 Relation to time correlation functions 517

On each plate the ¬‚uid will exert a drag force f per unit surface of the xz-

plane, proportional to the velocity gradient. The proportionality constant

is the viscosity coe¬cient ·:

‚vx

f =· . (18.44)

‚y

According to (17.113) on page 481 (see also Fig. 17.7), the force per unit

surface in a continuum is determined by the stress tensor σ:

dF = σ · dS. (18.45)

For example, on the lower plane in Fig. 18.2, with dS in the y-direction,

Fx = σxy Sy = σxy ”x”z. This phenomenological de¬nition agrees with the

de¬nition of · as given in the derivation of the Navier“Stokes equation in

Section 9.2, where · connects the o¬-diagonal elements of the stress tensor

with the velocity gradient as follows (see (9.11) and (9.12)):

‚u± ‚uβ

σ±β = · + (± = β). (18.46)

‚xβ ‚x±

This is the basic equation de¬ning ·. How the stress tensor can be deter-

mined from simulations is explained in Section 17.7.2 on page 484. The

o¬-diagonal elements of the average stress tensor over a given volume V are

equal to the negative o¬-diagonal elements of the pressure tensor,3 which is

“measured” by the virial:

1 1

Fi± xiβ = ’ σ±β (r) d3 r, ± = β.

P±β = (18.47)

V V V

i

The (Green“)Kubo relations for the viscosity coe¬cient can be found by

following the standard series of steps. We consider the xy-component of σ

without loss of generality.

(i) Apply a delta-disturbance g — δ(t) to ‚ux /‚y by imposing an additional

velocity vix = g — yi δ(t) to each particle in the considered volume.

After the delta pulse the ith particle is displaced in the x-direction

by ”xi = g — yi .

(ii) Compute ”H as a result of the disturbance, according to (18.24) and

using (18.47):

‚H

= ’g — Fxi yi = g — V σxy .

”H = ”xi (18.48)

‚xi

i i

3 The pressure tensor also contains a momentum transfer part due to particle velocities, but that

part is diagonal.

518 Linear response theory

(iii) De¬ne the response Y such that Y is proportional to ”H. This is

ful¬lled for

1

Y = σxy = ’ Fxi yi , (18.49)

V

i

for which ”H = g — V Y .

(iv) Find the δ-response function ¦(t) for Y from (18.25):

V

¦(t) = σxy (0)σxy (t) . (18.50)

kB T

The end result is

∞

V

σxy (0)σxy („ ) e’iω„ d„,

·(ω) = (18.51)

kB T 0

with the Kubo“Green relation for ω = 0:

∞

V

·0 = σxy (0)σxy („ ) d„. (18.52)

kB T 0