## Christopher G. Gray, Keith E. Gubbins, and Christopher G. Joslin

Print publication date: 2011

Print ISBN-13: 9780198556213

Published to Oxford Scholarship Online: January 2012

DOI: 10.1093/acprof:oso/9780198556213.001.0001

Show Summary Details
Page of

PRINTED FROM OXFORD SCHOLARSHIP ONLINE (www.oxfordscholarship.com). (c) Copyright Oxford University Press, 2017. All Rights Reserved. Under the terms of the licence agreement, an individual user may print out a PDF of a single chapter of a monograph in OSO for personal use (for details see http://www.oxfordscholarship.com/page/privacy-policy). Subscriber: null; date: 26 February 2017

# THERMODYNAMIC PROPERTIES OF PURE FLUIDS

Chapter:
(p.627) 6 THERMODYNAMIC PROPERTIES OF PURE FLUIDS
Source:
Theory of Molecular Fluids
Publisher:
Oxford University Press
DOI:10.1093/acprof:oso/9780198556213.003.0006

# Abstract and Keywords

This chapter gives an account of the statistical thermodynamics of pure fluids composed of non-spherical molecules. The early part of the chapter gives the derivation of equations for the various thermodynamic functions in terms of correlation functions and intermolecular forces, and includes discussion of quantum corrections and virial coefficients. The later parts of the chapter describe the application of molecular theory (perturbation theory, theory of hard non-spherical bodies, associating fluid theory) to the calculation of these properties, and compares these results with those from molecular simulations and experiments. The thermodynamics of nano-scale systems, in which some macroscopic laws and concepts may break down, are discussed at the end of the chapter. The mathematics of convex body geometry is given in an appendix.

[Statistical thermodynamics] differs from classical thermodynamics in that the thermodynamic functions for the assemblies and phases with which we deal are not left unspecified, or to be derived solely from measurements, but are always constructed a priori by the application to particular molecular models of the fundamental theorems of statistical mechanics with which we start.

Ralph Fowler and E. A. Guggenheim, Introduction to Statistical

Thermodynamics (1949)

In this chapter we consider the thermodynamic properties of a pure isotropic, homogeneous fluid in the classical, rigid molecule approximation (quantum corrections are briefly discussed in § 6.9). In the canonical ensemble the fundamental relation is that for the configurational free energy (the configurational part, in the classical limit, of the total free energy), given by Eq. (3.11) of Vol. 1,

$Display mathematics$
(6.1)

where Q c is the configurational partition function, given by (3.82):

$Display mathematics$
(6.2)

Here β = 1/kT, Ω = 4π or 8π2 for linear and nonlinear molecules, respectively, and Z c is the configurational phase integral. Expressions for the other thermodynamic properties in terms of the partition function can be derived from (6.1) using the usual relations of classical thermodynamics, e.g.

$Display mathematics$
(6.3)

$Display mathematics$
(6.4)

(p.628)

$Display mathematics$
(6.5)

From (6.1) and (6.3)–6.5 it is easy to write down the relations for the pressure p, configurational energy U c, and configurational heat capacity C vc in terms of the partition function Q c. If, in addition to the classical and rigid-molecule approximations we assume the potentials to be pairwise additive (Eq. (1.13)), then we can also relate these properties to the pair correlation function g(r ω 1 ω 2). These relations are derived in §§6.16.4 below. Alternative expressions for the pressure and chemical potential are provided by the compressibility relation (§ 6.3), and test particle expressions (also referred to as the potential distribution theorem, § 6.5), respectively. These expressions do not depend on either the rigid molecule or pairwise additivity approximations.

The energy, pressure, compressibility, free energy, and chemical potential equations ((6.9a), (6.15), (6.16), (6.24), (6.35), (6.37), and (6.49)) provide alternative routes to the thermodynamic properties. Four of these can be used to check the consistency of any theory of g(r ω 1 ω 2) by comparing results obtained by the four methods; these should be identical if the theory is exact, so that the degree of inconsistency gives a rough measure of the inexactness of the theory. Consistency is not a sufficient condition for exactness, of course; thus if we use a perturbation expansion for g(r ω 1 ω 2) identical results will be obtained by each of the four routes regardless of the order at which the series is truncated

In addition to the general expressions for p, U, χ, and A in terms of the partition function, intermolecular potential, and pair correlation function, it is often useful to express these properties in terms of the spherical harmonic expansion coefficients for the potential and the pair correlation function. Moreover, if the pair potential is of the site-site type some of these properties can be related to the simpler site-site correlation functions, rather than the full pair correlation function (in the case of χ such a relation exists irrespective of the type of potential). These expressions are therefore given in §§6.6 and 6.7. In §6.8 we note a rigorous inequality that exists for the free energy of a molecular fluid. Quantum corrections are briefly discussed in §6.9, and virial coefficients are considered in §6.10. Comparisons between theory and experiment are given in §6.11. The scaled particle approach for fluids of hard convex molecules, and its empirical extension to non-convex molecules, is described in §6.12, and theories for associating fluids are discussed in § 6.13.

# 6.1 The energy equation

The relation between the configurational energy and the partition function is obtained from (6.4) and (6.1) as

$Display mathematics$

(p.629) On performing the differentiation, we get:

$Display mathematics$
(6.6)

Using the pairwise additivity approximation and noting that the sum of the u(ij) ≡ u(r ij ω i ω j) contains ½N(N – 1) terms, each of which gives the same result on averaging, we have U c = ½N(N – 1) 〈u(12)〉, or

$Display mathematics$
(6.7)

where (3.105) has been used in carrying out the last step. This integral can be further simplified by noting that the integrand is a function of r 12 only, F(r 12) ≡ 〈u(12)g(12)〉ω 1 ω 2. On changing integration variables from (r 1, r 2) to (r 1, r 12), where r 12 = r 2r 1, we have

$Display mathematics$
(6.8a)

Expressing dr 12 in polar coordinates, $d r 12 = r 12 2 d r 12 d ω 12$, and carrying out the integrations over the angles ω 12 = θ 12 ϕ 12 gives

$Display mathematics$
(6.8b)

From (6.7) and (6.8b) we obtain the energy equation,

$Display mathematics$
(6.9a)

which provides a route to the thermodynamic properties if g(r ω 1 ω 2) is known. For example, it can be used to derive an equation for the pressure by first using U c = ∂(βA c)/∂β to derive an expression for A c, followed by the use of p = –(∂A c/∂V).

(p.630) For brevity (6.9a) is often written (since the integrand is independent of ω)

$Display mathematics$
(6.9b)

Similar relations hold for the pressure, free energy, and other equations given in subsequent sections.

# 6.2 The pressure equation

From (6.1) and (6.3) the pressure is given by

$Display mathematics$
(6.10)

or

$Display mathematics$
(6.11)

The evaluation of the volume derivative is complicated by the fact that V appears as the upper integration limit for the integrations over the r i. Such integrals have already been considered in §3.4.1. The differentiation can be performed by a change of integration variables suggested by H. S. Green,1

$Display mathematics$
(6.12)

so that (6.11) gives

$Display mathematics$
(6.13)

Assuming pairwise additivity of the potentials, we have

$Display mathematics$
(6.14)

If we substitute (6.14) into (6.13), note that each of the N(N – 1) terms in the ij sum gives the same result on averaging, and transform back to the original r i variables, we get p = ρkTN(N – 1)/6Vr 12 u′(12)〉, or

$Display mathematics$

(p.631) or

$Display mathematics$
(6.15)

where u′ = ∂u/∂r, and in the last step we have used (6.8a). Equation (6.15) is the pressure equation. The pressure equation provides a second route to the thermodynamic properties. Like the energy equation, (6.9a), it assumes rigid molecules and pairwise additivity. The pressure equation can also be derived from the virial theorem (see Eq. (E.14) of Vol. 1).

We note that in regions of r for which ∂u/∂r 〉 0 (attractive region) the contribution of the intermolecular potential term on the right-hand side of (6.15) will result in a decrease in pressure, while the reverse is true where ∂u/∂r 〈 0 (repulsive region).

# 6.3 The compressibility equation

The compressibility equation provides a third route to the thermodynamic properties, and has been written down in Eq. (3.113). It is most easily derived2 by differentiating the grand canonical expression for the first order distribution function f(r 1 ω 1) with respect to the chemical potential. Such derivatives have been considered in §§3.4.1 and 3.4.2, and (3.238) gives the general result of such a differentiation for the distribution function of order h. Putting h = 1 in (3.238) gives

$Display mathematics$
(6.16)

which is the compressibility equation. Here χ = ρ –1(∂ρ/∂p)T is the isothermal compressibility and g(r) is the centres pair correlation function. This expression does not rest on any assumptions concerning the form of the intermolecular potential energy (e.g. pairwise additivity or rigid molecules), and is therefore more general than the pressure equation. It holds for either isotropic or anisotropic fluids (e.g. liquid crystals), provided that they are homogeneous. Equation (6.16) relates χ to the centres correlation function, in contrast to the energy and pressure equations, where the full pair correlation function g(r ω 1 ω 2) is needed.

As shown in Appendix 3E, it is also possible to express χ in terms of the direct correlation function, c(r ω 1 ω 2), or the centres direct correlation function c(r) = 〈c(r ω 1 ω 2)ω 1 ω 2,

(p.632)

$Display mathematics$
(6.17)

which holds for an isotropic fluid.

# 6.4 The free energy equation

In order to relate the free energy to the pair correlation function g(r ω 1 ω 2) we use the coupling constant, or generalized charging process, approach.3 The intermole-cular potential energy is written

$Display mathematics$
(6.18)

so that λ = 0 corresponds to the reference system and λ = 1 to the real system. The configurational free energy A cλ for this system is

$Display mathematics$
(6.19)

where

$Display mathematics$
(6.20)

The difference AA 0 in free energies between the real and reference systems is given by

$Display mathematics$
(6.21)

Combining (6.19) and (6.20) and differentiating gives

$Display mathematics$
(6.22)

where 〈…〉λ indicates an ensemble average for a system with potential u 0 + λu p. For pairwise additive u p, (6.22) becomes (cf the derivation of (6.9a))

$Display mathematics$
(6.23)

where g λ is the pair correlation function corresponding to the system with pair potential u 0 + λu p. From (6.21), (6.22), and (6.23) we therefore have

$Display mathematics$
(6.24)

(p.633) which is the free energy equation. It provides a fourth route from the pair correlation function to the thermodynamic properties. For example, using p = –(∂A/∂V) and (6.24) yields a fourth expression for the pressure, in addition to those obtained from the energy, pressure, and compressibility equations. Equation (6.24) can be used in computer simulation studies to calculate free energy changes.4 The Monte Carlo values of (AA 0) given in Tables 4.1 and 4.2 were calculated in this way, as were the MSA and GMF values of (AA 0) obtained in §§5.4.4 and 5.4.9.

# 6.5 The test particle expressions for thermodynamic properties

## 6.5.1 The Widom method5

It is possible to relate the configurational thermodynamic properties to the interaction of a single molecule in the fluid with all of the other surrounding molecules.6 We call this single, specified molecule the ‘test particle’. These expressions are useful in both theoretical (see, for example, §6.12) and computer simulation (see later in this section) studies, and provide an alternative formalism to the expressions obtained in the previous four sections.

We first derive the expression for the chemical potential. Since we shall want to use these expressions for mixtures it is convenient to carry out the derivation for a mixture in which there are N A molecules of componentA, N B of component B, …, and N R of component R ; the total number of molecules is N = ∑β N β. The configurational part of the chemical potential µ αc of component α is given by

$Display mathematics$
(6.25)

where N′ means all N βα, and Q c for mixtures is (cf. (3.82))

$Display mathematics$
(6.26)

where x N = x 1 x 2x N with x i = r iωi, and dx i = dr ii. Since the N β are large we can approximate (6.25) by

$Display mathematics$
(6.27)

where Q c(N α) is given by (6.26) and Q c(N α – 1) is the corresponding partition function for the system of N – 1 molecules (one α molecule removed),

$Display mathematics$
(6.28)

(p.634) Here the product is over all components except α, u N–1 is the intermolecular potential energy for the (N – 1 ) molecule system, and we label the missing α molecule as molecule 1. The quantities Z c(N α) andZ c(N α–1) are the corresponding phase integrals given by (3.53). We now introduce u tα(x 1;x N–1), the intermolecular potential energy due to the interaction of the test particle of component α (molecule 1) with all other molecules in the system, by writing

$Display mathematics$
(6.29)

The potentials u and u N–1 can be expressed in terms of the coordinates of the molecules relative to those of the test particle, x 21 x 31x N1, where x j1 = r j1 ω j1 denotes positional (r j1 = r jr 1) and orientational ω j1 coordinates of molecule j relative to molecule 1; these coordinates are then independent of the coordinate x 1. We can therefore write Z c(N α) as

$Display mathematics$
(6.30)

where 〈…〉N–1 means an ensemble average over the (N – 1) molecule system with intermolecular potential energy u N–1 (x 2x N). From (6.27) and (6.30) we have7 , 8

$Display mathematics$
(6.31)

where ρ α = N α/V is the number density of α molecules. The quantity 〈exp[–βu ]〉 can be interpreted as follows. The test particle is placed with some fixed coordinates x 1 = r 1 ω 1 in the fluid of N – 1 molecules, and exp[–βu ] is then obtained as a function of the coordinates of all of the N – 1 molecules; a Boltzmann-weighted average over the coordinates x 2x N is then taken. Alternatively, we could imagine fixing the N – 1 molecules and moving the test particle freely among them; an unweighted average would then be taken over all possible positions and orientations of the test particle. In either case the test particle acts as a fictitious or ‘ghost’ particle, since it has no influence on the other N – 1 molecules.

(p.635) It is often more convenient to work with residual thermodynamic properties rather than configurational ones. The residual chemical potential μ αr is given by $μ α r = μ α − μ α id = μ α c − μ α c id$, where the superscript id signifies the ideal gas value at the same temperature, density, and composition as the real mixture. For an ideal gas u N = 0 and $Q c id = V N / ∏ β N β !$ from (6.26), so that, using Stirling’s approximation for ln N β!,

$Display mathematics$
(6.32)

we have

$Display mathematics$
(6.33)

so that8

$Display mathematics$
(6.34)

From this result and (6.31) we have

$Display mathematics$
(6.35)

An alternative but closely related expression for μ αr can be obtained by noting that

$Display mathematics$
(6.36)

where 〈…〉N is an ensemble average over the N molecule system for which the intermolecular potential energy is u N. From (6.35) and (6.36) we have9

$Display mathematics$
(6.37)

Note that in (6.37) the average is taken over a system in which the test particle does interact with the other N – 1 molecules, i.e. the test particle is not a ‘ghost’. We also note that both (6.35) and (6.37) are quite general. They do not rely on the assumption of pairwise additivity, nor on the rigid molecule approximation provided we suitably interpret the ensemble average, and are not restricted to uniform fluids7 They provide yet another route (the μ-route) to the thermodynamics; note that μ α is not expressed in terms of g(12), however10

Both (6.35) and (6.37) have been used in computer simulation studies to obtain the chemical potential11, 12 In using (6.35) it is common to introduce a fixed lattice of test particles having various orientations, and to obtain the interaction energy of the test particles with the real molecules; the average value of exp[–αu ] is then found. The major contribution to this average is from (p.636) configurations in which the interaction energy is attractive (u negative), and the method works well for low and moderate densities. At high densities (much above ρσ3 = 0.5 for Lennard-Jones (LJ) molecules12) the test particles overlap real molecules for the majority of configurations, leading to predominantly repulsive configurations; at such densities the usual Monte Carlo and molecular dynamics procedures do not effectively sample the regions of phase space most important in determining μ αr, and the statistical errors become large. One way of overcoming this problem is to bias the sampling so as to effectively cover the desired regions of phase space12 Similar problems are encountered in using (6.37) in simulations. In that case the dominant contributions to the average of exp[βu ] are from configurations in which U is repulsive (i.e. positive); since the test particle is a real molecule moving under the influence of forces from its neighbours, it tends to spend much of its time in attractive configurations, again leading to large statistical uncertainties in μ αr as calculated from the simulation. The situation can be illustrated schematically by rewriting (6.35) and (6.37) in the forms (for a pure fluid):

$Display mathematics$
(6.38)

$Display mathematics$
(6.39)

where f(u t) is the probability density for the test particle energy u t in the N – 1 molecule system, and g(u t) is the corresponding probability density in the N molecule system. As seen in Fig. 6.1, at high densities it is not enough to know f(u t) or g(u t) in the region of the maximum, since there are substantial contributions to the integrand for regions of u t where f or g is small.

FIG. 6.1 Probability density distributions needed in the test particle method for a typical liquid state condition (schematic): (a) moderate density, (b) and (c) high density. The function f is a broad distribution and extends to large U t values; the width of this distribution is compressed here for convenience of plotting.

(p.637)

FIG. 6.2 Comparison of experiment (lines) and Monte Carlo results (points) for the coexisting liquid density (ρ = ρσ 3) and vapour pressure (ρ = ρσ 3/ε) for diatomic Lennard-Jones fluids with bond-length = /σ = 0.630. For Cl2 the site-site parameters ε/k = 185 K, σ = 3.538 Å while for Br2 ε/k = 257.2 K, σ = 3.538 Å. These parameters gave the best fit to the available coexistence data. Units of P are the same in both figures. (From, Powles,11 based on simulation data of Romano and Singer.11)

An example of the use of such computer simulation results for μ r is shown in Fig. 6.2 for two-centre LJ models of chlorine and bromine. In this case μ r is calculated from (6.35), and the vapour pressure and coexisting liquid density are calculated by solving μ 1 = μ g and p 1 = p g = p. The gas-phase properties μ g and p g are obtained from the virial expansion. Since the conventional Monte Carlo sampling procedure is used it is only possible to obtain the coexistence properties for the liquid at rather high temperatures (low densities).

It should be noted that for fluids of hard bodies u can only take the values 0 or ∞. The average of exp(–βu ) that appears in (6.35) is still well behaved under these conditions (see §6.12 for applications to hard bodies). However, the average of exp(βu ) in (6.37) is indeterminate, so that this expression is not useful for such fluids.

It is possible to derive equations for other properties in terms of u (1) = u (x 1;x 2x N) by methods that closely follow that used above for the chemical potential6 For example, the configurational energy is u c = 〈u NN. If the inter-molecular potential energy is pairwise additive we can write

$Display mathematics$
(6.40)

and it is then easy to show that13

(p.638)

$Display mathematics$

or, using (6.27), (6.33), and (6.35),

$Display mathematics$
(6.41)

By similar arguments, we can derive the following equation for the pressure, starting from (6.13)and (6.14):

$Display mathematics$
(6.42)

where

$Display mathematics$
(6.43)

Equation (6.42) has been used to determine pressures for hard core systems,14 including the tensorial components of the pressure for inhomogeneous systems,15 and for continuous potentials16

The surface tension can also be expressed in terms of a test particle equation,17 and the pair correlation function g αβ(r ω 1 ω 2) can be expressed in terms of $u t α β ( 2 ) ( x 1 x 2 ; x 3 ⋯ x N )$, the intermolecular potential energy of interaction of two test particles of species α and β with all of the remaining molecules of the system.6,7

## 6.5.2 The Kirkwood method

Starting from (6.27) it is possible to derive a further expression for the chemical potential, due originally to Kirkwood,3 , 18 that involves the pair correlation function g(12) ≡ g(r ω 1 ω 2). In this approach molecule 1 of species α is coupled to the remaining N – 1 molecules by a coupling parameter λ. Equation (6.29) is now replaced by19

$Display mathematics$
(6.44)

where 0 ≤ λ ≤ 1 and x N–1 represents the configuration of all the other (N – 1) molecules. Thus λ = 0 and λ = 1 correspond to the N – 1 and N molecule systems, respectively. We have Z c(λ = 1) = Z c(N α) andZ c(λ = 0) = VΩα Z c(N α – 1), so that from (6.27) and (6.34) we obtain

$Display mathematics$
(6.45)

(p.639) where

$Display mathematics$
(6.46)

From (6.44)–6.46 we get

$Display mathematics$
(6.47)

We now introduce the approximation of pairwise additivity. We label the molecules of speciesA as 1,2,…, N A, those of B as 1,2,…, N B, etc., so that

$Display mathematics$
(6.48)

Substituting this in (6.47), noting that the sum over j contains (N βδ αβ) terms (each giving the same result on integration), and using the definition of g αβ(12) given in (3.252), we get

$Display mathematics$
(6.49)

which is the Kirkwood result. Here g αβ(12; λ) is the pair correlation function for β molecules around a central α1 molecule when the α1 molecule is coupled to the extent λ. This equation is less general than (6.35) and (6.37) of the previous subsection, since it is based on the assumptions of pairwise additivity and rigid molecules. However, it is a useful additional route to the chemical potential, and has been widely used in the theory of solutions20

# 6.6 Thermodynamic properties in terms of spherical harmonic expansion coefficients

The energy, pressure, compressibility, and free energy equations each involve integrals over functions of intermolecular distances and orientations. It is often convenient to expand these functions in spherical harmonics, and thus obtain an expression for the property in terms of the harmonic coefficients of g(r ω 1 ω 2) and u(r ω 1 ω 2). The property harmonic expansions are most readily derived using the generalized Parseval theorem; this is derived in Appendix B of Vol. 1. We consider here only the case of linear molecules, and employ space-fixed axes. The corresponding expressions for nonlinear molecules, and for the intermolecular r-frame, or for the space-fixed and k-frame expansions in k-space, can be obtained by using the appropriate form of the generalized Parseval theorem given in Appendix B.

(p.640) For linear molecules and space-fixed axes, the generalized Parseval theorem is

$Display mathematics$
(6.50)

where A(r ω 1 ω 2) and B(r ω 1 ω 2) are real and A( 1 2 ;r) is the space-fixed harmonic coefficient defined by

$Display mathematics$
(6.51)

where ω iθ i ϕ i is the direction of the symmetry axis of molecule i, ωθϕ is the direction of r, and C is a Clebsch-Gordan coefficient (App. A).

We can now apply (6.50) to the energy equation (6.9b) by taking A = u and B = g, so that

$Display mathematics$
(6.52)

where we have made use of the fact that u( 1 2 ;r) and g( 1 2 ;r) are real (see §§ 2.3 and 3.2.1) and depend only on r = |r|, and dr(…) here can be replaced by 4π r 2 dr(…). We note that the spherical harmonic expansion has reduced the original sevenfold integral to the sum of one-dimensional integrals in (6.52). Similarly, applying (6.50) to (6.15), (6.16), (6.17), and (6.24) gives

$Display mathematics$
(6.53)

$Display mathematics$
(6.54)

$Display mathematics$
(6.55)

$Display mathematics$
(6.56)

where u′( 1 2 ;r) = du( 1 2 ;r)/dr and h(000;r) = g(000;r) – (4π)3/2. In deriving (6.54) and (6.55) we note that B = 1, so that B( 1 2 ;r) = (4π)3/2 δ 10 δ 20 δ 0, and only the (000) harmonic of B is nonvanishing.

It is also possible to write down the corresponding expansions for U c, p,χ, and A in terms of the intermolecular-frame harmonic cosoft12sofefficients, by using the appropriate form of the generalized Parseval theorem, (B.85), in place of (6.50).

(p.641) Thus, for the configurational energy we obtain21

$Display mathematics$
(6.57)

where A( 1 2 m;r) is defined by

$Display mathematics$
(6.58)

where $ω i '$ refers to the r-frame and m ≡ –m. The space-fixed and r-frame coefficients of A are related by (3.147) and (3.148).

The intermolecular harmonics have been extensively used in computer simulation studies22 However, if the intermolecular potential consists of a small number of harmonic terms, the space-fixed harmonic expansion gives a simpler form for the properties. For example, if the anisotropic potential is u a = u QQ, the quadrupole-quadrupole interaction, then there is only one nonvanishing u( 1 2 ;r) harmonic, namely u(224; r), (see (2.177) and (2.24)):

$Display mathematics$
(6.59)

Thus, the expression (6.52) for the anisotropic contribution to U c, for example, reduces to a single term,

$Display mathematics$
(6.60)

so that only one space-fixed harmonic of g, namely g(224; r), need be known. In the intermolecular-frame expansion, (6.57), for the same intermolecular potential one would need five harmonics u( 1 2 m;r) and g( 1 2 m;r), viz. 1 2 m = 222,221,220,221,222. Also, since ru′ = –5u for the u QQ interaction we see that $p a = 5 3 ρ U c a / N$ for this case.

The convergence of the spherical harmonic expansion using the intermolecular frame has been studied23 for homonuclear diatomic molecules for an atom-atom LJ potential model together with molecular dynamics simulation results for the g( 1 2 m;r). It is found that for moderate elongations, 〈 0.45 (with = /σ where is the bond length and σ is the atom diameter), the series converges quite rapidly, but that the convergence quickly gets poorer as increases above this value. Typical results are shown in Table 6.1. The convergence has been found22 to be rapid for central LJ + quadrupole-quadrupole (QQ) potentials. For potentials that are quite hard, and elongated, the convergence worsens considerably22, 24 Of course, if the potential consists of only a few harmonic terms (e.g. LJ + QQ) speed of convergence is irrelevant for many properties, such as the energy, pressure, and G 2, which involve only particular low order harmonics. (p.642)

Table 6.1 Convergence of the spherical harmonic expansion for homonuclear diatomic molecules using an atom−atom Lennard-Jones model

U*/N

p *

*

T*

ρ *

Q*

Number of harmonics

Series

Exact

Series

Exact

0.425

1.543

0.50

0

4

−13.06

−13.11

−2.41

−2.34

0.425

1.543

0.50

0

8

−13.18

−13.11

−2.84

−2.34

0.5471

1.29

0.56

0.50

−25

−16.74

−14.45

−1.11

−0.61

0.5471

1.29

0.52

1.00

25

−13.18

−14.13

−0.08

−1.38

0.5471

2.20

0.52

2.00

25

−13.07

−16.91

−0.63

−0.67

(*) Reduced quantities are =/σ; T = kT/ε; ρ = ρσ 3;Q = Q/(εσ 5)½; U = U/ε; p = /ε, where σ and ε are the atom-atom LJ parameters. Here is bond length.

Four terms include coefficients 1 2 m = (000), (200), (020) and (220); eight terms include in addition to these (400), (040), (420), and (240); twenty-five terms include the additional harmonics (440), (600), (620), (260), (640), (460), (660), (800), (820), (280), (840), (480), (860), (680), (880), (10,0,0), and (0,10,0). From ref. 23.

# 6.7 Thermodynamic properties in terms of the site-site correlation functions

In this section we derive relations for the thermodynamic properties U c, p, χ, and A c in terms of site-site correlation functions. For the compressibility such relations can be derived for an arbitrary intermolecular potential, and are as general as the compressibility equation itself. The energy, pressure, and free energy equations can only be transformed to relations involving the g αβ correlation functions if the intermolecular potential is a sum of site-site potentials, as in (2.5).

## 6.7.1 Compressibility

We first consider the compressibility, χ. We have seen in § 3.1.6 that h αβ(k), the Fourier transform of h αβ(r) = g αβ(r) – 1, is given by

$Display mathematics$
(6.61)

where h(k ω 1 ω 2) is the transform of h(r ω 1 ω 2) and r ciα is the position of site α with respect to the molecular centre of molecule i. Taking the limit k → 0 in this relation gives (we use the notation$f ˜ ( k )$ in place of f(k) when there is a possibility of ambiguity, as in (6.62))

$Display mathematics$
(6.62)

and combining this result with (3.113) gives the compressibility equation

(p.643)

$Display mathematics$
(6.63)

or

$Display mathematics$
(6.64)

Hence χ can be obtained from any one of the site-site correlation functions g αβ (equivalently, χ is unaffected by the choice of molecular centres used for g(r) in (6.16), as expected physically).

It is possible to derive a similar expression for χ in terms of the reference interaction site model (RISM) site-site direct correlation functions c αβ, as follows. The (matrix) OZ relation between h αβ(k) and c αβ(k) is (5.238),

$Display mathematics$
(6.65)

where ω(k) is the matrix with elements ω αα(k) = sin(kℓ αα)/kℓ αα. In the limit k → 0, ω approaches the matrix with elements ω αα = 1 for all α, α′.

Hence from (6.65) we have:

$Display mathematics$
(6.66)

where

$Display mathematics$
(6.67)

and $c ˜ α β ( 0 ) ≡ c α β ( k = 0 )$. Hence we have

$Display mathematics$
(6.68)

(p.644) Thus the compressibility equation can be written

$Display mathematics$
(6.69)

or

$Display mathematics$
(6.70)

The number of $c ˜ α β$ in the sum over αβ sites in (6.69) and (6.70) corresponds to the number of $h ˜ α β$ in (6.65). In the event that only one site A is chosen in each molecule, e.g. at the molecular centre, (6.70) reduces to

$Display mathematics$

as one expects from the atomic liquid case, (5.20). It should be stressed that c AA(r) defined by the OZ equation (6.65) is not equal to c(r), where c(r) ≡ 〈c(r ω 1 ω 2)〉ω 1 ω 2 is the centres direct correlation function defined by (3.122), although their r-integrals are equal. In contrast, for this case h AA(r) is equal to 〈h(r ω 1 ω 2)〉ω 1 ω 2 see also §§ 3.1.5 and 5.5.2).

For homonuclear diatomic AA molecules with two sites per molecule, (6.64) and (6.70) reduce to

$Display mathematics$
(6.71)

which can also be derived directly from (5.251).

## 6.7.2 Internal energy 25–27

We now turn to the properties U, A, and p for the particular case of a site-site pair potential, (2.5). The energy equation, (6.9a), becomes

$Display mathematics$
(6.72)

where Ω = dω and r αβ is the distance between site α in molecule 1 and site β in molecule 2. The integral in this equation can be written more simply in terms of the site-site correlation function g αβ(r αβ), which is proportional to the probability density that sites α and β (on different molecules) are at a distance r αβ apart, irrespective of the molecular orientations (see § 3.1.6). To carry out this simplification we transform the integration variables in (6.72) from (r 1 r 2 ω 1 ω 2) to (p.645)

FIG. 6.3 Vectors r 1 and r 1α locate the centre and site α of molecule 1, respectively; XYZ are the space-fixed axes, and xyz are the body-fixed axes (kept parallel in the transformation of centres).

(r 1α r 2β ω 1 ω 2), where r i is the position of the molecular centre and r , is the position of site α in molecule i (see Fig. 6.3); ω i = ϕ i θ i χ i (or θ i ϕ i in the case of linear molecules) is the molecular orientation of the body-fixed axes xyz with respect to space-fixed axes XYZ, and is unchanged in this transformation. The probability of finding molecules 1 and 2 in the element of configuration space (dr 1 dr 2 dω 1 dω 2) at (r 1 r 2 ω 1 ω 2) is P(r 1 r 2 ω 1 ω 2) dr 1 dr 2 dω 1 dω 2, where P(r 1 r 2 ω 1 ω 2) is the specific pair distribution function (see §3.1.3). This probability is equal to that of finding the variables (r 1α r 2β ω 1 ω 2) in their corresponding ranges, i.e.

$Display mathematics$
(6.73)

Converting these specific distribution functions to (generic) correlation functions g (see (3.91) and (3.105)) we have, for the isotropic, homogeneous case,

$Display mathematics$
(6.74)

where g(r ω 1 ω 2) is the usual angular pair correlation function and g(r αβ ω 1 ω 2) is proportional to the probability density that a pair of molecules are at orientations ω 1 and ω 2, with site-site separation r αβ = r r . The site-site correlation function is related to g(r αβ ω 1 ω 2) by (cf. (3.126))

$Display mathematics$
(6.75)

Since (6.74) is true for each αβ term in the sum in (6.72) we can substitute (6.74) in (6.72); note that u αβ (r αβ) is independent of ω 1 and ω 1, and use (6.75) to get

$Display mathematics$
(6.76)

(p.646) Changing variables from r 1α r 2β to r 1α r αβ, integrating over r 1α, and using polar coordinates for r αβ (see (6.8a) and (6.8b)) gives

$Display mathematics$
(6.77)

Thus a knowledge of the site-site correlation functions is sufficient to calculate the configurational energy, provided that the pair potential is of the site-site form. For the special case of homonuclear diatomic (AA) molecules, (6.77) becomes

$Display mathematics$
(6.78)

For the case of atomic liquids, having one site per molecule, (6.77) reduces to the standard result

$Display mathematics$
(6.79)

## 6.7.3 Pressure25–28

For a site-site potential the pressure equation, (6.15), becomes

$Display mathematics$

or, using (6.74)

$Display mathematics$
(6.80)

The derivative ∂u αβ (r αβ)/∂r 12 is at fixed molecular orientations and fixed ω 12, and is given by

$Display mathematics$
(6.81)

To evaluate (∂r αβ/dr 12) we note that (see Fig. 3.8)

$Display mathematics$
(6.82)

where r cαβ = r c2βr c1α. Thus we have

$Display mathematics$
(6.83)

where $r ^ 12 = r 12 / r 12$ is the unit vector in the r 12 direction. Differentiating (6.83), noting that r cαβ and r cαβ are fixed (since ω 1 and ω 2 are fixed), we have

$Display mathematics$

(p.647) Thus, rearranging and using (6.82), we have

$Display mathematics$
(6.84)

where $r ^ 12 = r 12 / r 12$ and $r ^ α β = r α β / r α β$ are the unit vectors along r 12 and r αβ, respectively, and cos $γ α β = r ^ 12 ⋅ r ^ α β$ (i.e. γ αβ is the angle between r 12 and r αβ). Substituting (6.81) and (6.84) into (6.80), and using (6.75), we get

$Display mathematics$
(6.85)

where 〈…〉rαβ denotes a weighted average over the orientations keeping the site-site distance r αβ fixed, i.e.

$Display mathematics$
(6.86)

Changing variables in (6.85) from r 1α r 2β to r 1α r αβ, integrating over r 1α, and using polar coordinates for r αβ gives

$Display mathematics$
(6.87)

where $u α β ′ ( r α β ) = d u α β ( r α β ) / d r α β$. This is the pressure equation for a site-site potential. Equation (6.87) is formally of the site-site type, but since 〈r 12 cos γ αβr αβ requires the complete g(r ω 1 ω 2) for its evaluation, it is not completely analogous to the relations (6.64) and (6.77) and (6.98) for χ, U c, and A c. Thus p cannot be expressed directly in terms of g αβ even for a site-site potential29 For an atomic liquid, with one site per atom, (6.87) reduces to the usual result, (5.23). For diatomic AA 4 molecules (6.87) becomes

$Display mathematics$
(6.88)

We now consider the special case where the molecules are composed of hard spheres rigidly joined together (with sphere α of diameter σ αα, etc.). The site-site potential is therefore

$Display mathematics$
(6.89)

(p.648) where σ αβ = ½ (σ αα + σ ββ). It follows that

$Display mathematics$
(6.90)

and

$Display mathematics$
(6.91)

where δ(r αβσ αβ) is the Dirac delta function. Equation (6.91) follows because exp[–βu αβ(r αβ)] ≡ θ(r αβ) is the unit step function. The relation dθ(x)/dx = δ (x) is verified by noting that the defining conditions (δ(x) = 0 at x ≠ 0, δ(x) = ∞ at x = 0, and d(x) = 1) for the δ-function are satisfied (see Appendix B).

It follows from (6.91) that

$Display mathematics$

or

$Display mathematics$
(6.92)

Substituting this result into (6.85), and using the fact that y αβ (r αβ) = exp[βu αβ(r αβ)]g αβ(r αβ) is continuous, gives

$Display mathematics$
(6.93)

where $g α β ( σ α β + )$ is the value of the site-site correlation function at contact.

For a hard sphere atomic fluid we have 〈r 12 cos γ αβr αβ = σ = σ, and (6.93) reduces to the usual result, (5.24). For a fluid of hard homonuclear diatomic AA molecules, each of the terms in the summation over α and β in (6.93) is the same, and the pressure equation becomes

$Display mathematics$
(6.94)

Using (6.82) and (6.84) we find that

$Display mathematics$
(6.95)

where = /σ, = 2r c1A = 2r c2A being the site-site distance, and $θ 1 A ′$ and $θ 2 A ′$ are the polar angles corresponding to r c1A and r c2A with Z chosen to lie along r AA (see Fig. 6.4). (p.649)

FIG. 6.4 Coordinates for two hard AA molecules.

## 6.7.4 Free energy 28,30

In § 6.4 we have shown that (see (6.21), (6.22))

$Display mathematics$
(6.96)

where A is the Helmholtz free energy for the full system with potential energy u(r N ω N), A 0 is the value for a reference system with potential energy, u 0(r N ω N), u puu 0 is the perturbing energy, and 〈…〉λ is an average over configurations with Boltzmann weight corresponding to u λ = u 0 + λu p. If u p is pairwise additive, and the intermolecular pair potentials u p(ij) are of site-site form, then from (6.77) we have

$Display mathematics$
(6.97)

and substituting this result in (6.96) gives

$Display mathematics$
(6.98)

where $g α β λ ( r )$ is the site-site pair correlation function for the system with potential energy u λ = u 0 + λu p. If the reference system is the ideal gas, then $u α β p$ is replaced by the full site-site potential u αβ in the above expressions. For diatomic AA molecules (6.98) becomes

(p.650)

$Display mathematics$
(6.99)

The above method of introducing the charging parameter λ is not suitable for fluids composed of fused hard sphere molecules. For such fluids it is more convenient to use a charging parameter that scales both28 the site-site diameters σ αβ and the intramolecular centre-site distances r . Thus molecules characterized by the value λ of this parameter have site-site diameters $σ α α λ$ and intramolecular distances $r c α λ$ defined by

$Display mathematics$
(6.100)

so that λ = 0 corresponds to the ideal gas of point noninteracting molecules. By a derivation analogous to that of (6.24) in § 6.4 we find that

$Display mathematics$
(6.101)

where A 0 is the Helmholtz free energy of the ideal gas of point particles at the same temperature, density, and number of molecules as the fluid of fused hard spheres. We note that (6.101) cannot be reduced to a form involving site-site correlation functions, since the term $∂ u α β λ ( r α β λ ) / ∂ λ$ involves the molecular orientations ω 1 and ω 2 (see below).

We must now evaluate $∂ u α β λ / ∂ λ$. The potential $u α β λ$ depends on λ through both $σ α β λ$ and $r α β λ$, since (see (6.82))

$Display mathematics$
(6.102)

where r cαβ = r c2βr c1α. Thus we have

$Display mathematics$
(6.103)

From (6.92) we get

$Display mathematics$
(6.104)

(p.651) Similarly $∂ u α β λ / ∂ r α β λ$ is given by (6.92). To obtain $∂ r α β λ / ∂ λ$ we use (6.83),

$Display mathematics$

where (6.82) and (6.100) have been used. Thus we have

$Display mathematics$
(6.105)

where $r ^ α β = r α β / r α β$. Substituting these results into (6.103), noting that $∂ σ α β λ / ∂ λ = σ α β$ from (6.100), and substituting the resulting form for $∂ u α β λ / ∂ λ$ into (6.101) gives

$Display mathematics$
(6.106)

We note that the first term on the right-hand side of this expression can be written in terms of $g α β λ$, because the term $( ∂ u α β λ / ∂ σ α β λ ) ( ∂ σ α β λ / ∂ λ )$ in (6.103) is independent of ω 1, ω 2; the same is not true of the second term on the right in (6.106) since r caβ depends on ω 1 and ω 2. In (6.106) both $g α β λ$ and g λ(r ω 1 ω 2) refer to the fluid with potential u λ. By using the identity p = ∂A/∂V, together with the virial theorem, it is possible to derive28 the equation of state (6.93) from (6.106).

An alternative expression for the free energy can be obtained30 for fused hard sphere fluids by using a charging parameter λ′ that scales only the site diameters σ αβ, so that (6.100) is replaced by

$Display mathematics$
(6.107)

and the centre-site distances r cα (and the angles between the r cα) remain unchanged during the charging process and equal to the real values. In this case λ′ = 0 corresponds to an ideal gas composed of fused hard sphere molecules with infinitesimally small spheres. This gas, whose Helmholtz free energy is written A0′ ≡ A(λ′ = 0), differs from that used as the reference state in (6.106) in that the molecules possess rotational kinetic energy. By retracing the steps in deriving (6.106) it is clear that in this case the free energy will be given by

(p.652)

$Display mathematics$
(6.108)

It should be noted that $g α β λ ′ ( σ α β λ ′ + )$ in (6.108) is the site-site correlation function for molecules with variable site diameters $σ α β λ ′$, but with the fixed real values of the r ; in (6.106), by contrast, $g α β λ ( σ α β λ + )$ is the site-site correlation function when both the $σ α β λ$ and $r c α λ$ are variable. Equations (6.106) and (6.108) must give identical results if the exact correlation functions are used. If an approximate theory is used for the correlation functions they provide alternative routes to the thermodynamic properties which may give different results.

For a fluid of hard homonuclear diatomic AA molecules, (6.108) becomes

$Display mathematics$
(6.109)

# 6.8 A rigorous inequality for the free energy31–41

In this section we discuss an application of the classical Gibbs inequality31 for the configurational Helmholtz free energy. The quantum version of the inequality (due to Bogoliubov) can also be derived;32 35 the result then applies to the entire free energy, rather than just the configurational part.

If we write the total intermolecular potential energy as the sum of a reference and a perturbation part, u λ = u 0 + λu p, then the free energy difference between the real and reference systems is given by (6.96). Moreover, we have

$Display mathematics$
(6.110)

which is the analogue of (3.249) in the canonical ensemble. The right-hand side of (6.110) must be negative or zero. It follows that 〈u pλ is a decreasing function of λ, or a constant. Hence the maximum and minimum values of the integrand in (6.96) are 〈u pλ = 0 and 〈u pλ = 1, respectively. Multiplying these by the integration interval in (6.96) (unity), we arrive at bounds for the integral itself,

$Display mathematics$
(6.111)

where 〈u p〉 = 〈u pλ = 1. The lower bound 〈u p〉 in (6.111) is probably of little computational value, since it requires knowledge of the real fluid distribution functions. The upper bound

$Display mathematics$
(6.112)

(p.653) is valid for angle-independent or dependent u p, and is referred to as the Gibbs-Bogoliubov inequality.32 , 33 , 36 Atomic fluid applications are discussed by Mωnster36 When 〈u p0 is nonzero it can be used39 , 40 as a basis for a variational determination of AA 0.

If the fluid is molecular and the Pople reference system is chosen (§ 4.5), then 〈u p0 = 0 and

$Display mathematics$
(6.113)

i.e. the anisotropic perturbation rigorously lowers the free energy41 The inequality (6.113), which holds to all orders and for arbitrarily large Pople-type u p, is manifestly satisfied for multipolar perturbations in second order (see (4.49)). For this case we also have, from (6.111), that 〈u p〉 ≤ 0; this is expected, since the Boltzmann averaging places more weight on the negative parts of u p.

The quantity 〈u p0 equals the first-order perturbation term A 1 (see (4.5)), so that A 1 gives a rigorous upper bound for the free energy difference AA 0.

It is clear from (6.110) that (for β ≠ 0) the equalities in (6.111) hold if and only if u p is a constant, equal to its average value for all configurations. If 〈u p0 = 0 the only possible value is u p = 0. For β → 0 (infinite temperature) the situation is more subtle;41 we can have AA 0 = 0 or AA 0 = – ∞, depending on the nature of the potential.

General inequalities for the entropy or internal energy do not exist; it is possible that such inequalities can be derived for special potentials41

# 6.9 Quantum corrections42

When quantum corrections are small, as is expected for most liquids except He and H2, (see § 1.2.2 and Table 1.1 of Vol. 1), their effect on the thermody-namic properties can be treated by expanding the partition function in powers of ħ2. Details of this expansion are given in Appendix 3D of Vol. 1, and by Powles and Rickayzen42 The corresponding expansion for the Helmholtz free energy, A = –kT ln Q, gives the following results for molecules of the symmetries indicated:

Linear (I x = I yI):

$Display mathematics$
(6.114)

Spherical top (I x = I y = I zI):

$Display mathematics$
(6.115)

(p.654) Symmetrie top (I x = I yI⊥, I zI ):

$Display mathematics$
(6.116)

Asymmetric top (I x, I y, I z not equal):

$Display mathematics$
(6.117)

where A is the quantal and A C the classical value of A, and 〈…〉 is a classical configurational average. Each of these expressions contains three terms of order ħ 2; the first, involving 〈F 2〉, arises from translational (diffraction) effects, the second involves 〈τ 2〉 from rotational potential energy effects, and a final term (e.g. ħ 2 N/6I in (6.114)) arises from rotational kinetic energy effects (see § 1.2.2 for a qualitative discussion). Here 〈F 2〉 and 〈τ 2〉 are the classical mean squared force and mean squared torque (about the centre of mass) on a molecule, respectively, and can be obtained experimentally from the isotope separation factor (see § 1.2.2) or from spectral moments (infrared, Raman, or neutron), as discussed in Chapter 11. These quantities can also be calculated theoretically or by computer simulation (see Chapter 11). The quantities m and I α are the molecular mass and principal αα-component of the moment of inertia, respectively, and in (6.117)

$Display mathematics$
(6.118)

Also, τ α, is the principal α component of τ,$〈 τ ∥ 2 〉 ≡ 〈 τ z 2 〉$, and $〈 τ ⊥ 2 〉 = 〈 τ x 2 〉 + 〈 τ y 2 〉$ is the total mean squared torque perpendicular to the symmetry axis z.

For linear molecules, (6.114) can be rearranged to the dimensionless form

$Display mathematics$
(6.19a)

or, alternatively,

$Display mathematics$
(6.119b)

where F 2 = F 2/(ε/σ)2, τ 2 = τ 2/ε 2, and I = I/ 2 are the reduced quantities, with ε and σ being the usual energy and distance parameters in the isotropie part of the potential, and m the molecular mass; the quantities Λ = h/2πσ()½ and θ r = ħ 2/2Ik are the dimensionless de Broglie wavelength and the characteristic rotational temperature, respectively, introduced in § 1.2.2. The magnitude of these corrections for liquids with dipole-dipole and (p.655)

FIG. 6.5 Quantum correction to the free energy for a liquid of linear molecules with pair potential u LJ + u a, where u LJ is the Lennard-Jones potential and u a is either the dipole-dipole (curves labelled μ) or quadrupole-quadrupole (curves labelled Q) term; the state condition 0.001 and 0.07, is ρ = 0.800,T = 0.719. Solid and dashed curves are for I = I/ 2 = 0.001 and 0.07, respectively.

Table 6.2 Quantum corrections to configurational energy and pressure for ortho-baric liquids

Liquid

T(K)

$U − U c ℓ U c ℓ$

$p − p c ℓ ρ k T$

N2

73

−5

0.8

F2

59

−1.7

0.26

Cl2

250

−0.7

0.07

HCl

188

−2

0.16

$CH 4 †$

207

−4.4

0.09

† This state point is not orthobaric. The density is 24.9 mol −1. From ref. 42.

quadrupole-quadrupole forces is shown in Fig. 6.5, based on Monte Carlo simulation values43 , 44 of 〈F 2〉 and 〈τ 2〉. Values ofI for most linear molecules lie in the range 0.001 (HI) to 0.07 (I2), and results for these two extremes are shown in Fig. 6.5. The relative magnitudes of the three terms in (6.119a) are strongly dependent on the value ofI and on the strength of the anisotropic forces. Thus, when I ~ 0.001 it is clear that the rotational quantum corrections are the dominant ones unless the electrostatic forces are very weak. Conversely, for large I, e.g. I ~ 0.07, the translational correction dominates even for relatively strong electrostatic forces.

Rough estimates of the quantum corrections to the configurational energy and pressure have been made for several liquids by Powles and Rickayzen,42 and are shown in Table 6.2. These calculations are based on computer simulation values (p.656)

FIG. 6.6 Density of hydrogen from experiment48 (points), theory neglecting quantum effects (---), and theory including the O(ħ 2) quantum correction (–) as given by (6.114). (From ref. 45.)

of 〈τ 2〉 and 〈F 2〉 for site-site potential models; for nitrogen the quadrupole-quadrupole potential was also included. The ratio of the rotational to the trans-lational potential energy terms in (6.119a), i.e. 〈τ 2m/〈F 2〉I, is about 0.6 for N2 and F2, and about 4 for HCl. The results show that the corrections to U are small, but not negligible. For the pressures the quantum corrections are larger. Furthermore, the quantum corrections to the pressure may be a large fraction of the total pressure for liquids near their triple points, since the total pressure is then small. For orthobaric liquid N2 at 73 K, for example,p ≈ 1 bar, whereas pp cℓ = 140 bar. For compressed fluid CH4 at 207 K, 24.9 mol ℓ–1 (well away from the triple point), on the other hand,p ≈ 818 bar, while pp cℓ = 69 bar.

For hydrogen the quantum corrections are particularly large, so that the ħ 2 term alone in (6.114) is not sufficient at low temperatures, and a full quantal (p.657)

FIG. 6.7 Configurational internal energy of hydrogen. Key is as described in legend to Fig. 6.6. (From ref. 45.)

treatment becomes necessary45 , 46 Nevertheless, the ħ 2 correction significantly improves the accuracy of calculated thermodynamic data for the temperature range 100–200 K45 Figures 6.6 and 6.7 show such results, compared to a purely classical calculation, for equation of state data. The calculations are based on the pair potential model

$Display mathematics$
(6.120)

(p.658) where u LJ(r) is the Lennard-Jones model of (2.2), u QQ( 224) is the quadrupole-quadrupole potential given by (2.185), u dis ( 202 + 022 + 224) is the anisotropic dispersion potential as given by the 1 2 = 202, 022, and 224 terms (see (2.223) and also (2.228) and (2.232)), and u ov(202 + 022) is the anisotropic overlap potential as given by the 1 2 = 202 and 022 terms (see (2.247), (2.248)). Calculations are based on the Padé approximant of (4.47), with the anisotropic potential parameters Q = 0.65 × 10–26 esu, K = 0.125 (see (D.8)), and δ = 0.1 (values of Q and K are taken from sources listed in Appendix D of Vol. 1). Values of δ, ε, and σ are obtained by matching theory with experimental density and configurational energy data at higher temperatures, where quantum effects are small. It is seen from Figs. 6.6 and 6.7 that the quantum corrections increase at lower temperatures and higher pressures, as expected. Quantum corrections are appreciable at temperatures below 160 K, particularly at the higher pressures.

Quantum corrections to the critical constants have been calculated for LJ fluids47 The quantum corrections cause the critical temperature, pressure, and density to decrease.

# 6.10 Virial coefficients

The second virial coefficient is given by (3.268) and (3.269) as

$Display mathematics$
(6.121a)

$Display mathematics$
(6.121b)

where u′ = ∂u/∂r. These expressions are based on the assumptions that (a) the molecules are ‘rigid’ (see § 1.2.1), (b) quantum effects can be neglected (§ 1.2.2), and (c) the angle-averaged pair potential, 〈u(r ω 1 ω 2)〉ω 1 ω 2, decreases faster than r –3 for large r. Assumption (a) can lead to errors for molecules that have internal rotations or are otherwise flexible, while (b) may not hold for light gases at low temperature. Nevertheless, (a) and (b) lead to negligible errors for most small molecules that do not contain hydrogen. Condition (c) is necessary in order that the integrals appearing in (3.268) and (3.269) be independent of the shape and size of the system, and is required if the thermodynamic properties are to exist. Condition (c) will hold for all cases considered in this book; note that (c) is satisfied for a potential of the form u 0(r) + u μμ, where u 0 is shorter-ranged than r –3 and u μμ is the dipole-dipole potential (which varies as r –3), since 〈u μμω 1 ω 2 = 0.

Qualitatively B 2(T) behaves as shown in Fig. 6.8 as a function of temperature T. At low temperature the attractive forces dominate molecular behaviour. On (p.659)

FIG. 6.8 General form of the reduced second virial coefficient, $B 2 * = B 2 / b 0$ with b 0 = 2πσ 3/3, as a function of reduced temperature, T = kT/ε. The curve shown is for a classical LJ gas.

average the molecules spend more time in each others ‘spheres of influence’ than they would in the absence of interaction. Since the molecules then hit the walls less frequently,p is reduced below the ideal gas or kinetic value ρkT, so that B 2 is negative. At high temperatures the opposite effect occurs; the molecules bounce off each other, due to the repulsive forces, hit the walls more frequently, and hence increase p over ρkT, giving a positive B 2. The rapid decrease with temperature of B 2 (T) at very low temperatures is caused by more and more (exponentially more) molecules becoming mutually trapped in bound states. The slow decrease with temperature of B 2 (T) at very high temperatures is due to the slowly decreasing effective diameter of the molecules; as the collisions become more energetic, the molecules penetrate further into the core region of the potential.49

## 6.10.1 B2 and the pair potential

For the majority of cases, where the above conditions of rigidity and classical behaviour are fulfilled, (6.121a) shows that measurements of B 2 provide a valuable source of information about the pair potential. Complications due to the influence of three-body potentials or the need for an approximate theory, always a concern in interpreting liquid state data, are not present here. Information about the pair potential can be extracted from B 2 data either by inversion of the data, or by comparing calculated B 2 values with the data for various model potentials. The first procedure is the most desirable, and is possible in principle for atomic gases.

(p.660) Thus, if the molecules are spherical, and the pair potential has a monotonically repulsive portion u +(r +) for r = r +r m, and a monotonically attractive portion u (r ) for rr r m (cf. Fig. 2.4), where rm is the separation at the potential minimum, it is possible to rearrange (6.121b) to the form50 , 51

$Display mathematics$
(6.122)

where ϕu + ε is the pair potential referred to the minimum of the well as zero, ε is the well depth, ℒ[ ] is the Laplace transform, and:

$Display mathematics$
(6.123)

Thus, if accurate B 2 data are available over a range of temperatures, it should be possible to invert the Laplace transform in (6.122) to obtain Δ (ϕ), i.e. the potential itself for the purely repulsive region (ϕε), and the well-width function $( r + 3 + r − 3 )$ for ϕε. Such a procedure has been successfully used to obtain the repulsive part of the potential for helium;52 for this case data are available over a very wide range of temperatures and the attractive forces are weak. This inversion procedure has proved difficult for other atomic gases, since the inverse transform is highly sensitive to inaccuracies in the data, and because in many cases data are not available over a sufficient temperature range53 Less formal methods of inversion based on semiempirical techniques have been devised, and have proved successful for atomic gases51 , 53 , 54

For molecular gases it is not possible to write (6.121a) or (3.272) as a Laplace transform that could be used to obtain information about u(12), so that a formal inversion procedure does not exist in this case55 The use of B 2 data to investigate potentials has therefore been by the second of the methods mentioned earlier, i.e. comparisons of experiment with calculated B 2 values using various models. Unfortunately, it is found that B 2 is only sensitive to the form of the intermolecular potential energy function when data are available over a wide temperature range. When data are available over only a limited range B 2 is relatively insensitive to the form of the potential (in particular the anisotropic part), so that a variety of models often give an equally good fit to the data. Thus, for nitrogen and carbon dioxide it is found that the Kihara, LJ atom-atom (two sites), and Gaussian overlap potential models, with or without the addition of a quadrupole-quadrupole term, all give an equally good description of the available second virial coefficient data, although they do not work equally well for other properties (e.g. heat of sublimation, crystal lattice parameter)56 This situation is largely due to the insensitivity of the temperature-dependent pseudopotential u 0(r), defined by (3.273), to the nature of the anisotropic intermolecular forces; thus both multipolar and nonspherical shape forces tend to have a similar effect on u 0(r).57 Indeed, this insensitivity is (p.661)

FIG. 6.9 Second virial coefficient of HCl: • experiment;57 – – site–site Lennard-Jones (SSLJ) fitted to liquid; –– SSLJ + multipole terms fitted to B 2;…….. second model fitted to liquid. (From ref. 63.)

the basis for simple and widely used corresponding states methods for predicting second virial coefficients58

Second virial coefficient data have often been used to estimate values for anisotropic potential parameters, by comparison of the data with calculations for a particular model; multipole moments and anisotropic overlap parameters have often been obtained in this way, by using a generalized Stockmayer potential model59 Multipole moments so obtained are often unreliable because of uncertainties in the potential model used (see Appendix D of Volume 1). A more stringent test of the potential model is to compare theory and experiment for both the pressure and the dielectric second virial coefficients60 (see § 10.6, and also the collision-induced absorption virial coefficients, Chapter 11), and where possible the mean squared torque for the gas61 (see Chapter 11).

Although second virial data alone are of limited value in testing potential models, a model that fits other data should, of course, also fit the virial data if it is the correct pair potential. A test of this sort on two potential functions whose parameters were obtained by comparing computer simulation results with experimental thermodynamic data for liquid hydrogen chloride is shown in Fig. 6.9. Either of these potential functions can be made to fit the B 2 data if the site-site parameters ε αβ and σ αβ are suitably adjusted (here α, β are H or Cl); the solid line shows the result for the second potential, which includes multipole terms. When the parameters are adjusted to fit liquid data, however, the potentials no longer fit the virial data. The potential including multipole terms appears to be the better of the two. The failure of the liquid potentials to fit the gas data could be due to

(p.662) errors in the pair potential, the neglect of multibody potential terms in treating the liquid, or a combination of these two effects.

## 6.10.2 Calculations for simple model potentials

There have been many calculations of B 2, and in some cases higher virial coefficients, for a wide variety of simple model potentials of the types discussed in § 2.1; the earlier work has been reviewed elsewhere51 , 64 , 65 These have included calculations for the generalized Stockmayer model, site-site LJ models with and without point charges, the Kihara model, Gaussian Overlap model, etc. (These models are discussed in §2.1 of Vol. 1.) The addition of electrostatic forces causes B 2 to decrease at any given temperature. This is because such forces are, on the average, attractive, and cause a decrease in the frequency of molecules hitting the wall, and hence a decrease in the pressure.

Fluids of hard nonspherical molecules are of particular interest since they are frequently used as reference fluids in perturbation theories for both isotropic and anisotropic (e.g. liquid crystal) fluids. For a fluid of hard spheres of diameter a it is easy to show from (6.121a) that

$Display mathematics$
(6.124)

where v m = πσ 3/6 is the volume of a sphere. For a fluid of hard nonspherical molecules, each of volume v m, it is convenient to write

$Display mathematics$
(6.125)

where f B is a factor whose difference from unity indicates the degree of departure from hard sphere behaviour.

The simplest molecules to treat are convex ones, i.e. ones in which any line segment whose endpoints are inside the body lies wholly within the body (see Appendix 6A); examples of such bodies are spherocylinders and ellipsoids (a nonconvex example is two fused hard spheres). For such molecules f B is always greater than unity, so that convex nonspherical molecules always have second virial coefficients larger than those for spherical molecules of the same volume. The physical reason for f B 〉 1 is that, as the molecule rotates, the effective volume swept out is greater than that of the molecule, thus tending to increase B 2. Isihara66 has derived the following general expression for f B for the convex molecule case (see Appendix 6A.5 for derivation):

$Display mathematics$
(6.126)

where γ = r m s m/v m is the shape factor, and s m and s m are the mean radius of curvature and the surface area of the molecule, respectively (see Appendix 6A for an account of the geometry of hard convex bodies). From (6.125) and (6.126) we have, for hard convex molecules

$Display mathematics$
(6.127)

(p.663) We note that γ, and hence B 2, is independent of temperature as in the case of hard spheres. For spheres of radius a we have r m = a, s m = 4πa 2, v m = 4πa 3/3, γ = 3, and f B = 1. Formulae for r m, s m, v m, and f B for other shapes have been worked out;65 , 66 some of these are given in Table 6A.1.

For general nonconvex molecules there is no simple expression analogous to (6.127). For homonuclear (AA) hard sphere diatomic molecules (dumbbells), Isihara has derived an exact analytic expression for B 2.67 This result has been generalized to heteronuclear hard diatomics (AB), unlike pair diatomics (AB/CD), and symmetrical linear triatomics (ACA) by Wertheim.68 Numerical calculations of B 2 have been made for other fused hard sphere molecules, and som