the derivative of the potential, and therefore the potential is no longer con-

servative. As illustrated in Fig. 6.5a, the derivative of a truncated potential

contains a delta function at the cut-o¬ radius. Incorporation of this unphys-

ical delta function into the force leads to intolerable artifacts. The use of

a truncated force without delta function (Fig. 6.5b) implies that the e¬ec-

tive potential function in the simulation is not the truncated real potential

but a shifted potential obtained by integration of the truncated force. One

may then expect that equilibrium ensembles generated by dynamic trajec-

tories di¬er from those generated by Monte Carlo simulations based on the

truncated potential.

But even truncated forces with shifted potentials generate artifacts. When

particles di¬use through the limit of the interaction range they encounter

a sudden force change leading to extra noise, to heating artifacts and to

artifacts in the density distributions. The discontinuity of the force causes

errors when higher-order integration algorithms are used that rely on the

existence of force derivatives. Such e¬ects can be avoided by shifting the

force function to zero at the cut-o¬ distance (Fig. 6.5c), but this has an even

more severe in¬‚uence on the e¬ective potential, which now deviates from

the exact potential over a wide range. Several kinds of switching functions,

switching the force smoothly o¬ at the cut-o¬ radius, are used in practical

MD algorithms. The user of such programs, which seemingly run error-free

even for short cut-o¬s, should be aware of the possible inadequacies of the

e¬ective potentials.

Of course, the error due to neglect of the long-range interaction beyond

the cut-o¬ radius rc can be reduced by increasing rc . This goes at the

3

expense of longer pair lists (which scale with rc ) and consequently increased

computational e¬ort. In addition, rc should not increase beyond half the

smallest box size in order to restrict interactions to the nearest image in

periodic systems. One should seek an optimum, weighing computational

e¬ort against required precision. But what is the error caused by the use of

a modi¬ed interaction function?

Let us consider a homogeneous ¬‚uid with an interatomic dispersion in-

teraction v disp = ’C6 r’6 and compute the correction terms to energy and

pressure for three commonly used short-range interactions v sr : truncated

potential, truncated force, shifted force. The average number density is ρ

and the radial distribution function is g(r): given the presence of a par-

ticle at the origin, the probability of ¬nding another particle in a volume

element dr equals ρ g(r) dr. Noting that the correction ”v(r) involves the

full dispersion interaction minus the employed short-range interaction, the

6.3 Force ¬eld descriptions 161

correction to the potential energy and therefore to the internal energy u per

particle is

∞

1

”v(r)4πr2 g(r) dr,

”u = ρ (6.33)

20

and the pressure correction is

∞

2π d”v(r)

”P = ’ ρ2 r3 g(r) dr, (6.34)

3 dr

0

as is easily derived from the virial expression for the pressure (see Section

17.7.2 on page 484)

1

r ij · F ij ,

P V = N kb T + (6.35)

3

i j>i

with r ij = r i ’ r j and F ij is the force exerted by j on i:

dv(r) r ij

F ij = ’ . (6.36)

dr rij

r=rij

We obtain the following results, assuming that g(r) = 1 for r >= rc and

that the number of particles within rc 1:

(i) Truncated potential

v sr = ’C6 r’6 , (6.37)

2π ’3

”u = ’ ρC6 rc , (6.38)

3

’3

”P = ’2πρ2 C6 rc . (6.39)

The pressure correction consists for two-thirds of a contribution from

the missing tail and for one third of a contribution due to the discon-

tinuity of the potential at the cut-o¬ radius.

(ii) Truncated force

v sr = ’C6 r’6 + C6 rc ,

’6

(6.40)

4π ’3

”u = ’ ρC6 rc , (6.41)

3

4π 2 ’3

”P = ’ ρ C6 rc . (6.42)

3

(iii) Shifted force

v sr = ’C6 r’6 ’ 6C6 rrc + 7C6 rc ,

’7 ’6

(6.43)

7π ’3

”u = ’ ρC6 rc , (6.44)

3

162 Molecular dynamics

7π 2 ’3

”P = ’ ρ C6 rc . (6.45)

3

In order to avoid di¬culties caused by the discontinuity in the potential,

even for MC simulations where the pressure is a¬ected, the potential should

never be truncated, neither in MD nor in MC. In the derivation of the cor-

rections we have neglected details of the radial distribution function in the

integrals, which is justi¬ed when rc extends beyond the region where g(r)

di¬ers from 1. Thus ”u depends only “ and linearly “ on density, with

the consequences that the change in Helmholtz free energy equals ”u, that

there is no change in entropy, and that ”P = ρ”u. The e¬ect of the cor-

rection is appreciable, especially a¬ecting vapor pressure and the location

of the critical point.7 For example, for the Lennard“Jones liquid argon 8 at

85 K, not far from the boiling point (87.3 K), the energy correction for a

shifted force with rc = 1 nm is ’0.958 kJ/mol, and the pressure correc-

tion is ’20.32 kJ mol’1 nm’3 = ’337 bar! In Table 6.2 the corrections in

the thermodynamic properties at this particular density-temperature state

point are compared to the thermodynamic values themselves. The e¬ects

are non-negligible, even for this rather long cut-o¬ radius of about 3 σ.

For water with9 C6 = 0.002617 and at room-temperature liquid density the

shifted force correction for 1 nm cut-o¬ is ’0.642 kJ/mol for the energy

and ’357 bar for the pressure. Such corrections are essential. Still, they

are usually not applied, and the models are parameterized to ¬t empirical

data using MD with a given cut-o¬ method. It is clear that the model

parameters are then e¬ective parameters that incorporate the e¬ect of re-

strictions of the simulations; such parameters must be readjusted when other

cut-o¬s or long-range methods are applied. While such e¬ective potentials

are convenient (commonly they do not only imply the e¬ects of long-range

treatment but also of lack of polarizability and other contributions that are

not pair-additive, neglect of speci¬c interaction terms and neglect of quan-

tum e¬ects), they tend to restrict the generality and transferability of force

¬elds.

At this point it is worthwhile to remark that for dispersion (and other

power-law interactions) long-range contributions can be evaluated under pe-

7 See for the Lennard“Jones equation-of-state and a discussion on truncation of LJ interactions:

Nicolas et al. (1979), Smit (1992), Johnson et al. (1993) and Frenkel and Smit (2002). The

critical temperature for the truncated force shifts as much as 5% downward when the cut-o¬

is decreased from 3.5 to 3σ.

σ = 0.34 nm; µ = 119.8 kB = 0.9961 kJ/mol; C6 = 4µσ 6 = 0.006155 kJ mol’1 nm6 . Density at

8

85 K and 1 bar: 35.243 mol/dm3 (number density 21.224 nm’3 ).

9 Value for the SPC (Berendsen et al., 1981) and SPC/E model (Berendsen et al., 1987), which

is the value recommended by Zeiss and Meath (1975), based on experimental data.

6.3 Force ¬eld descriptions 163

Table 6.2 Corrections to thermodynamic properties of liquid argon at 85 K

and a density of 35.243 mol/dm’3 , if “measured” by isobaric/isochoric