fusion coe¬cient is easily found from either of these equations. The steady

state is reached exponentially with a time constant equal to M D/RT . The

amplitudes a1 and a2 should be chosen large enough for a measurable e¬ect

and small enough for a negligible disturbance of the system; in practice the

imposed velocity di¬erences should not exceed about 10% of the thermal

velocities.

18.5.3 Thermal conductivity

The thermal conductivity coe¬cient Λ can be measured by imposing a ther-

mal ¬‚ux Jq with the system™s periodicity:

Jq (y) = A cos ky. (18.66)

This is accomplished by scaling the velocities of all particles every time step

by a factor ». The kinetic energy Ekin per unit volume changes by a factor

522 Linear response theory

»2 , such that

(»2 ’ 1) 3

”Ekin

Jq = = ρkB T. (18.67)

”t ”t 2

The external heat ¬‚ow causes a temperature change; the temperature T (y)

will obey the following equation:

‚2T

‚T

ρcv = Jq + Λ 2 . (18.68)

‚t ‚y

The steady-state solution then is

A

T (y) = T0 + cos ky. (18.69)

Λk 2

Thus, by monitoring the Fourier coe¬cient of the temperature at wave vector

k in the y-direction, the thermal conductivity coe¬cient Λ is found. Also

here, the amplitude of the heat ¬‚ux should be chosen large enough for a

measurable e¬ect and small enough for a negligible disturbance. In practice

a temperature amplitude of 10 K is appropriate.

Exercises

18.1 When the delta response of a linear system equals an exponential

decay with time constant „c , compute the response to a step function

and the frequency response (18.9) of this system.

18.2 Verify the validity of he Kramers“Kronig relations for the frequency

response of the previous exercise.

18.3 Prove (18.15) by following the same reasoning as in the proof given

for (18.16).

18.4 Prove (18.23) by showing through partial integration that

‚Y ‚H

e’βH(z) ”zi e’βH(z) Y ”zi

dz = β dz. (E18.1)

‚zi ‚zi

When the total dipole ¬‚uctuation M (0) · M (t) appears to decay

18.5

exponentially to zero in a simulation with “tin-foil” boundary con-

ditions, how then do the real and imaginary parts of the dielectric

constant depend on ω? Consider µr (ω) = µr ’ iµr . Plot µr versus µr

and show that this Cole“Cole plot is a semicircle.

18.6 Derive (18.62).

19

Splines for everything

19.1 Introduction

In numerical simulations one often encounters functions that are only given

at discrete points, while values and derivatives are required for other values

of the argument. Examples are:

(i) the reconstruction of an interaction curve based on discrete points,

for example, obtained from extensive quantum calculations;

(ii) recovery of function values and derivatives for arbitrary arguments

from tabulated values, for example, for potentials and forces in MD

simulations;

(iii) the estimation of a de¬nite or inde¬nite integral of a function based

on a discrete set of derivatives, for example, the free energy in ther-

modynamic integration methods;

(iv) the construction of a density distribution from a number of discrete

events, for example, a radial distribution function.

In all these cases one looks for an interpolation scheme to construct a

complete curve from discrete points or nodes. Considerations that in¬‚uence

the construction process are (i) the smoothness of the curve, (ii) the accuracy

of the data points, and (iii) the complexity of the construction process.

Smoothness is not a well-de¬ned property, but it has to do with two

aspects: the number of continuous derivatives (continuous at the nodes),

and the integrated curvature C, which can be de¬ned as the integral over

the square of the second derivative over the relevant interval:

b

{f (x)}2 dx,

C[f ] = (19.1)

a

523

524 Splines for everything

or, in the multidimensional case:

{∇2 f (r)}2 dr.

C[f ] = (19.2)

V

The term “curvature” is used very loosely here.1 If data points are not in-

¬nitely accurate, there is no reason why the curve should go exactly through

the data points. Any curve from which the data points deviate in a statis-

tically acceptable manner, is acceptable from the point of view of ¬tting to

the data. In order to choose from the “ generally in¬nite number of “ ac-

ceptable solutions, one has to apply additional criteria, as compliance with

additional theoretical requirements, minimal curvature, minimal complexity

of the curve speci¬cation, or maximum “uncertainty” (“entropy”) from an

information-theoretical point of view. Finally, one should always choose the

simplest procedure within the range of acceptable methods.

The use we wish to make of the constructed function may prescribe the

number of derivatives that are continuous at the nodes. For example, in an

MD simulation using an algorithm as the Verlet or leap-frog scheme, the

¬rst error term in the prediction of the coordinate is of the order of (”t)4 .

The accuracy depends on the cancellation of terms in (”t)3 , which involve

the ¬rst derivatives of the forces. For use in such algorithms we wish the

derivative of the force to be continuous, implying that the second deriva-

tive of the potential function should be continuous. Thus, if potential and

forces are to be derived from tabulated functions, the interpolation proce-

dure should not only yield continuous potentials, but also continuous ¬rst

and second derivatives of the potential. This, in fact, is accomplished by

cubic spline interpolation. Using cubic spline interpolation, far less tabu-

lated values are needed for the same accuracy than when a simple linear

interpolation scheme would have been used.

We ¬rst restrict our considerations to the one-dimensional case of polyno-

mial splines, which consist of piecewise polynomials for each interval. These

local polynomials are by far to be preferred to global polynomial ¬ts, which

are ill-behaved and tend to produce oscillating solutions. The polynomial

splines cannot be used when the function is not single-valued or when the

x-coordinates of the data points cannot be ordered (x0 ¤ x1 ¤ · · · ¤ xn );

one then needs to use parametricn splines, where both x and y (and any

1 It would be better to call the curvature, as de¬ned here, the total energy of curvature. Curvature

is de¬ned as the change of tangential angle per unit of length along the curve, which is the

inverse of the radius of the circle that ¬ts the local curved line segment. This curvature is a

local property, invariant for orientation of the line segment. For a function y(x) the curvature

can be expressed as y (1 + y )’3/2 , which can be approximated by y (x). If an elastic rod is

bent, the elastic energy per unit of rod length is proportional to the square of the curvature.

See Bronstein and Semendjajew (1989) for de¬nitions of curvature.

19.2 Cubic splines through points 525

n

t

yn