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 0 and a perturbation part and we write, using λ as a bookkeeping label,
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
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
The first derivative
and higher derivatives can be used to construct a Taylor series2) which reads
and is usually addressed as the λ-expansion. Again, 〈f〉0 means the average over the reference ensemble using only, but later we will use also 〈f〉1 for the average over the complete ensemble using . Since , the λ2 term is always positive so it appears that is an upper bound for F. In fact, we can much further progress if we note that for q ≥ 0, qlnq − q + 1 ≥ 0 and write q = f/g. For arbitrary, normalized distributions functions f and g with ∫fdr = ∫gdr = 1 we then have
since g(r) ≥ 0, being a distribution function. Identifying the reference system (λ = 0) with and the actual system (λ > 0) with we obtain
Reversing the roles of f and g, the result is
In total we thus have, what is called the Gibbs–Bogoliubov inequality ,
Figure 7.4 shows the relationship between the various terms, from which it becomes clear that knowledge of and provides bounds on Fλ − F0. Because is much more easily calculated than , the upper bound is usually considered. If contains parameters, can be minimized with respect to these parameters to obtain the optimum result. Some consideration makes it clear that the “flatter” , 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, can be simplified and to first order
Figure 7.4 Terms in the Bogoliubov inequality for (upper) and (lower). Arrow → indicates the effect of F′′, while arrow 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 and a perturbation part , so that use can be made of hard-sphere results. A more complete description of perturbation theory is given by Hansen and McDonald .
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](V − b) = 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 and consider the hard-sphere model with radius σ as the reference system.
a) Show that, if is small enough, then .
b) Show, using Eq. (7.58) and approximating gHS(r) by gHS(r) = 0 for r < σ and gHS(r) = 1 for r > σ, that 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 = (V − Nb)N with excluded volume b = ½·4πσ3/3, and show that this approximation leads to a pressure contribution βP = ρ/(1 − bρ).
d) Show that the final expression for the pressure P becomes βP = ρ /(1 − bρ) − baρ2 or [P + a/V2](V − b) = 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  (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
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
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 . 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
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 . 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 . 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  (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
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 . The effective hard-sphere diameter d, in this case a function of T and ρ, was chosen according to
where r0 is the equilibrium distance of the potential ϕ. An analytical expression fitted to simulation results is given by Verlet and Weis . 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  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. .
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.