Perturbation Theory - Modeling the Structure of Liquids: The Integral Equation Approach - Liquid-State Physical Chemistry: Fundamentals, Modeling, and Applications (2013)

Liquid-State Physical Chemistry: Fundamentals, Modeling, and Applications (2013)

7. Modeling the Structure of Liquids: The Integral Equation Approach

7.4. Perturbation Theory*

The EoS as obtained from the pressure and compressibility equations using the PY (and simulations) results for the hard-sphere liquid can be used as a basis for rather more accurate estimates of the EoS for more realistic potentials. This approach is usually addressed as the perturbation approach, and is based on the use of a reference potential for which the results are known and a perturbation part for which the effects are to be assessed. The starting point is the Gibbs–Bogoliubov inequality for approximate solutions, which we discuss first. Thereafter, we indicate two possible divisions of the potential in a reference and a perturbation part.

7.4.1 The Gibbs–Bogoliubov Inequality

Suppose we have a Hamilton expression consisting of a reference part c7-math-50470 and a perturbation part c7-math-5048 and we write, using λ as a bookkeeping label,

(7.49) c7-math-0049

We see that λ = 0 means that we have the reference potential, while λ = 1 implies that we have the full energy expression. The configurational partition function Q becomes

(7.50) c7-math-0050

where Q0 is the configurational partition function of the reference system and ⟨···⟩0 denotes the canonical average in that reference system. The configurational Helmholtz energy is accordingly

(7.51) c7-math-0051

The first derivative

(7.52) c7-math-0052

and higher derivatives can be used to construct a Taylor series2) which reads

(7.53) c7-math-0053

and is usually addressed as the λ-expansion. Again, ⟨f0 means the average over the reference ensemble using c7-math-5021 only, but later we will use also ⟨f1 for the average over the complete ensemble using c7-math-5049. Since c7-math-5022, the λ2 term is always positive so it appears that c7-math-5023 is an upper bound for F. In fact, we can much further progress if we note that for q ≥ 0, qlnqq + 1 ≥ 0 and write q = f/g. For arbitrary, normalized distributions functions f and g with ∫fdr = ∫gdr = 1 we then have

(7.54) c7-math-0054

since g(r) ≥ 0, being a distribution function. Identifying the reference system (λ = 0) with c7-math-5024 and the actual system (λ > 0) with c7-math-5025 we obtain

(7.55) c7-math-0055

Reversing the roles of f and g, the result is

(7.56) c7-math-0056

In total we thus have, what is called the Gibbs–Bogoliubov inequality [16],

(7.57) c7-math-0057

Figure 7.4 shows the relationship between the various terms, from which it becomes clear that knowledge of c7-math-5039 and c7-math-5040 provides bounds on FλF0. Because c7-math-5041 is much more easily calculated than c7-math-5042, the upper bound is usually considered. If c7-math-5043 contains parameters, c7-math-5044 can be minimized with respect to these parameters to obtain the optimum result. Some consideration makes it clear that the “flatter” c7-math-5046, the smaller the perturbation contribution will be. In case we have a pairwise additive potential Φ = ½Σijϕij, which can be split as ϕ(r) = ϕ0(r) + ϕ1(r) with ϕ0(r) a reference part and ϕ1(r) the perturbation, c7-math-5026 can be simplified and to first order

(7.58) c7-math-0058

Figure 7.4 Terms in the Bogoliubov inequality for c7-math-5037 (upper) and c7-math-5038 (lower). Arrow → indicates the effect of F′′, while arrow c7-fig-5001 represents the effect of minimization.


Higher-order terms involve three-point and higher order correlations functions, also for pair potentials. Since a great deal of information is known about the hard-sphere system, this is a good reference system; however, potentials like the LJ do not contain a hard core so that an effective value for the hard-sphere diameter must be chosen. In the sequel, we discuss two ways in which to choose the reference part c7-math-5050 and a perturbation part c7-math-5051, so that use can be made of hard-sphere results. A more complete description of perturbation theory is given by Hansen and McDonald [7].

Problem 7.7

Verify Eq. (7.58).

Problem 7.8: The van der Waals EoS

For the vdW fluid originally an intuitive “derivation” was given leading to [P + a/V2](Vb) = NkT, where a and b are parameters while P, T, and V represent the pressure, temperature, and volume, respectively, for N molecules. Using the perturbation approach, the vdW equation can be derived so that the approximations involved become clear. Use the separation c7-math-5027 and consider the hard-sphere model with radius σ as the reference system.

a) Show that, if c7-math-5028 is small enough, then c7-math-5029.

b) Show, using Eq. (7.58) and approximating gHS(r) by gHS(r) = 0 for r < σ and gHS(r) = 1 for r > σ, that c7-math-5030 where the minus sign is included to make a > 0.

c) Next, Z0 has to be chosen. Approximate, with van der Waals, Z0 by Z0 = (VNb)N with excluded volume b = ½·4πσ3/3, and show that this approximation leads to a pressure contribution βP = ρ/(1 − ).

d) Show that the final expression for the pressure P becomes βP = ρ /(1 − ) − baρ2 or [P + a/V2](Vb) = NkT.

7.4.2 The Barker–Henderson Approach

For a continuous potential, like the LJ potential, a choice for the hard-sphere diameter is not unique. In order to obtain an optimum choice Barker and Henderson [17] (BH) defined a modified potential v(r;α,γ,d), containing the three parameters α, γ and d, corresponding to ϕ(r) by


Here, σ represents the point where ϕ(σ) = 0. The parameters α and γ represent the “inverse steepness” and the “depth” of the potential, respectively. Note that v(r;α,γ,d) is independent of d, providing freedom to choose d optimally. Moreover v(r;α,γ,d) reduces to ϕ(r) if α = γ = 1, and reduces to the hard-sphere potential ϕHS(r) with diameter d if α = γ = 0. Consequently, by varying α and γ we can change the potential from ϕHS(r) to ϕ(r).

The next step taken was to expand the Helmholtz energy F in a power series in both α and γ around α = γ = 0, leading to

(7.59) c7-math-0059

Because the derivatives are evaluated at α = γ = 0, they represent reference – that is, hard-sphere – quantities. While omitting the differentiations [17, 18], the final result is

(7.60) c7-math-0060

where f(r) = exp[−βϕ(r)] − 1 denotes the Mayer function. The subscript 0 for the Helmholtz function F0, pair correlation function g0(r) and compressibility (∂ρ/∂P)0 indicates that they are quantities for the reference system with diameter d. The second-order term in γ has been approximated, as the three-point and four-point correlation functions are not available. Furthermore, Barker and Henderson made plausible that the second-order term in α2 is negligible. Since d is still free, we can choose it as c7-math-5034. In this way, d becomes a temperature-dependent3) diameter, but for all temperatures and densities the first-order term in α, as well as the second-order term in αγ, vanishes identically. Taking α = γ = 1, the result for the Helmholtz energy is

(7.61) c7-math-0061

where the functions F0, g0(r) and (∂ρ/∂P)0 are to be evaluated for d as chosen. In spite of being approximate, the second-order term can be evaluated quite accurately [19]. The agreement with results from simulations is quite good. We only show in Figure 7.5 results for the EoS of a calculation (to second order) for the LJ potential and the density of the coexisting phases [17]. One observes an excellent agreement with the (essentially exact) results of computer simulations. For further comparison, the reader is referred elsewhere [18, 19].

Figure 7.5 (a) The EoS and (b) density for coexisting phases for the LJ potential from second-order BH theory. The curves are labeled with T* values, while the points indicate results from computer simulations and experimental data. The quantities ρ* and T* are given by ρσ3 and kT/ε, respectively.


7.4.3 The Weeks–Chandler–Andersen Approach

Another successful perturbation theory is due to Weeks, Chandler and Andersen [20] (WCA), who wrote the potential as ϕ(r) = ϕ0(r) + ϕ1(r) with



where ε is the depth of the potential and r0 is the equilibrium distance (Figure 7.6). This, in principle, is a good choice as it makes the perturbation ϕ1(r) as slowly varying (“flat”) as possible. WCA used the first-order perturbation expression

(7.62) c7-math-0062

with F0 = FHS and g0(r) = exp(−βϕ0)yHS(r), where yHS(r) is the hard-sphere cavity function. Other choices for g0(r) are discussed elsewhere [21]. The effective hard-sphere diameter d, in this case a function of T and ρ, was chosen according to

(7.63) c7-math-0063

where r0 is the equilibrium distance of the potential ϕ. An analytical expression fitted to simulation results is given by Verlet and Weis [22]. On the one hand, it is not immediately clear how to choose the hard-sphere reference diameter. One choice is Eq. (7.63), but on the other hand WCA approximates converges faster so that second-order terms are not required. Both, the BH and WCA theories, however, provide good approximations. Figure 7.7 shows a comparison [23] of results for the pair correlation function of the PY, BH, and WCA theories, from which a good agreement can be seen. It should be noted, however, that minor differences in g0(r) can result in relatively large differences in thermodynamic properties. The optimum choice for the reference system is therefore rather important.

Figure 7.6 Separation of the intermolecular potential in WCA approximation. (a) The original potential; (b) Separation in a reference part ϕ0(r) and a perturbation ϕ1(r). Modified from Ref. [31].


Figure 7.7 The pair correlation function g(r) for the LJ fluid. The curves marked - - -, –−− and … show results of the PY theory, the BH theory, and the WCA theory, respectively. The points refer to results of simulation studies. The quantities ρ* and T* are given by ρσ3 and kT/ε, respectively.