THE JOURNAL OF CHEMICAL PHYSICS 142, 214110 (2015)

Computation of virial coefficients from integral equations Cheng Zhang, Chun-Liang Lai, and B. Montgomery Pettitt Department of Biochemistry and Molecular Biology, Sealy Center for Structural Biology and Molecular Biophysics, The University of Texas Medical Branch, Galveston, Texas 77555-0304, USA

(Received 4 February 2015; accepted 18 May 2015; published online 5 June 2015) A polynomial-time method of computing the virial coefficients from an integral equation framework is presented. The method computes the truncated density expansions of the correlation functions by series transformations, and then extracts the virial coefficients from the density components. As an application, the method was used in a hybrid-closure integral equation with a set of self-consistent conditions, which produced reasonably accurate virial coefficients for the hard-sphere fluid and Gaussian model in high dimensions. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4921790]

I. INTRODUCTION

The virial expansion1–7 of a fluid is a series for the pressure, P, in terms of the number density, ρ = N/V (N is the number of particles, V is the volume), βP =

+∞ 

Bn ρ n ,

(1)

n=1

where β = 1/(k BT) is the inverse temperature (k B is the Boltzmann constant), and Bn are the virial coefficients. The virial expansion is a unifying concept of liquid state theory,1–3,5,6 for it provides a simple equation of state for fluids. Accurate highorder virial coefficients are valuable, if the radius of convergence is sufficiently large. Even beyond the radius of convergence of the series, we can often use the virial coefficients to construct an equation of state by choosing a more suitable algebraic form with compatible low-density behavior.8 As is well appreciated, high-order virial coefficients are important measures for judging different liquid theories.8 There are many works in the literature8–16 for computing exact or approximate high-order virial coefficients, for liquids of various complexities (hard-spheres, water, polymers, proteins, etc.). While the low-order (n ≤ 4) virial coefficients can be computed accurately,12,17–20 the computation of higher-order ones remains challenging.8,21 This is so even for simplest systems such as the hard-sphere fluid,8,11,13,15,16,19,22 in which spheres interact via a repulsive hard-core potential [that is, the potential, φ(r), between two spheres of unit diameter separated by a distance r is zero if r > 1, or infinity if r < 1]. Despite its simplicity, the hard-sphere fluid is theoretically important, and the system and its virial coefficients have been intensively studied not only in two and three dimensions, but also in much higher dimensions.9,11,12,20,23,24 High-dimensional (particularly, hard-sphere) fluids are attractive theoretical models as the high dimensionality reduces the particle correlation and hence the complexity of the fluid. Studies on these fluids, and their virial coefficients,9,11,12,20,23,24 can show whether a theory or method that works for higher dimensions can be mapped or extrapolated to lower dimensions. For example, in the theory of integral 0021-9606/2015/142(21)/214110/12/$30.00

equations, the Percus-Yevick (PY) equation,25 which is excellent for the low-dimensional hard-sphere fluid, becomes less successful than the hypernetted-chain (HNC) equation26–28 in higher dimensions. One may use the tactic of studying the high-dimensional cases in hope of gaining some insight or techniques for lower (1-3) dimensional systems (e.g., by extrapolation along dimensionality24). In higher dimensions, reasonably accurate high-order virial coefficients for the hard-sphere fluid can be obtained15 from simulation by Mayer sampling.10,13,14,29,30 The resulting virial coefficients showed a regular pattern as the order n changes. Thus, we wish to consider whether they can be derived, at least approximately, by other theoretical means. One approach is to compute virial coefficients from integral equations, such as the PY25 or HNC.26–28 From such integral equations, one can readily obtain the density expansions of the correlation functions, and then extract the virial coefficients from the integrals of the density components.14 The high-order virial coefficients computed from the above procedure can be used to measure the accuracy of an integral equation. Particularly, for the three-dimensional hard-sphere fluid, many integral equations have produced similar low-order virial coefficients, as we shall see below. Thus, one has to know the virial coefficients to a sufficiently high order to judge which integral equation is better for a dense state. This approach also allows us to understand the observed pattern of the virial coefficients of the high-dimensional hard-sphere fluid through a sufficiently accurate integral equation. Despite the importance, it is generally nontrivial to compute the virial coefficients for a given integral equation theory. Although previous algorithms14 work reasonably for the PY equation, the HNC version was far more intensive computationally for high orders because of its exponential time complexity. Further, since neither the PY nor the HNC integral equation produces reliable virial coefficients for high orders or in higher dimensions, we wish to consider alternative routes. Here, we generalize and recast the previous algorithm of Shaul et al.14 to compute the virial coefficients from a larger class of integral equations. Our method has a polynomial time complexity, and thus computationally, it is particularly useful

142, 214110-1

© 2015 AIP Publishing LLC

214110-2

Zhang, Lai, and Pettitt

J. Chem. Phys. 142, 214110 (2015)

for high orders. As an application, we used the algorithm in a hybrid integral equation with a set of order-by-order self-consistent conditions. The method yields rather accurate virial coefficients of the high-dimensional hard-sphere fluid and Gaussian model. We describe the method in Sec. II, use it on the hard-sphere fluid and Gaussian model in Sec. III, and conclude the article in Sec. IV.

and the indirect correlation function as t(r) = ρ t 1(r) + ρ2t 2(r) + · · ·,

where c0(r) = f (r), and t 0(r) vanishes. The components of cl (r) and t l (r) are obtained below. Substituting Eqs. (7) and (8) into Eq. (3), and comparing the coefficients of ρl , yields, for l ≥ 1, t˜l (k) =

II. METHODS

l−1 

c˜l−1− j (k) [c˜j (k) + t˜j (k)].

(9)

j=0

A. Integral equations

In liquid state theory, the integral equation method5 is a way of computing the pair correlation function, g(r). The equation can be cast as a pair of coupled equations involving the direct correlation function, c(r), and the indirect correlation function, t(r) ≡ g(r) − c(r) − 1. For an isotropic fluid, as assumed here, the above correlation functions depend on the coordinates, r, only through the distance, r = |r|. The first of the coupled equations is the Ornstein-Zernike (OZ) relation t(r) = ρ (c ∗ h)(r) = ρ [c ∗ (c + t)](r),

(2)

where “*” denotes a spatial convolution and h(r) ≡ c(r) + t(r) = g(r) − 1 is the total correlation function. This relation is more convenient in Fourier space, as it assumes the product form t˜(k) = ρ c(k)[ ˜ c(k) ˜ + t˜(k)],

(3)

where the tilde indicates the Fourier transform. The second equation, or the closure, can be written in terms of the cavity distribution function, y(r) ≡ g(r) exp[ βφ(r)], c(r) = [1 + f (r)] y(r) − 1 − t(r),

(4)

where f (r) ≡ exp[− β φ(r)] − 1. The formally exact expression of y(r) is y(r) = exp[t(r) + B(r)], where B(r) is the bridge function, which cannot be written as a simple convolution expression in terms of pair functions. A closure without B(r) often approximates y(r) as a simple functional of t(r). The PY equation uses the following closure:25 y(r) ≈ 1 + t(r). Alternatively, the HNC equation assumes

(8)

(5) 26–28

y(r) ≈ exp t(r).

(6)

In either case, the closure along with the OZ relation, Eq. (2), determines the approximate correlation functions, c(r) and t(r). Below, we will illustrate the method with the PY and HNC equations, and the extension to other integral equations31–47 is straightforward. B. PY equation

We now compute the virial coefficients from the PY equation.14 First, we expand the direct correlation function as a series of the density, ρ, as c(r) = c0(r) + ρ c1(r) + ρ2c2(r) + · · ·,

(7)

Similarly, for closure Eq. (4), we have cl (r) = [1 + f (r)] yl (r) − t l (r) − δl0.

(10)

In the PY case, Eq. (5) means yl (r) = δl0 + t l (r), and thus, cl (r) = f (r) δl0 + f (r) t l (r).

(11)

Equations (9) and (11) can be applied in succession, F.T.

Eq. (9)

cl−1(r) −−−→ c˜l−1(k) −−−−−→ F.T.

Eq. (11)

t˜l (k) −−−→ t l (r) −−−−−−→ cl (r), to get higher-order density components of the correlation functions. For example, in the first iteration, Eq. (9) gives t˜1(k) = [c˜0(k)]2 = [ f˜(k)]2, i.e., t 1(r) = ( f ∗ f )(r); then, Eq. (11) gives c1(r) = ( f ∗ f )(r) · f (r). Note that the correlation functions on the right-hand side of Eq. (9) have orders lower than l; thus, they are always available when we compute t˜l (k). 1. Virial coefficients

We can compute the virial coefficients from the correlation functions, cl (r) and t l (r), in a few ways. First, in the compressibility route,  (c) ∂ρ ( βP ) = 1 − ρ c(r) dr, (12) where ∂ρ ≡ ∂/∂ ρ. Then, expanding both sides as series of ρ using Eqs. (1) and (7) yields  1 cn−2(r) dr. (13) Bn(c) = − n Next, in the virial route, we have, in D dimensions,  ρ2 β (v) (r · ∇φ) g(r) dr βP = ρ − 2D  ρ2 = ρ+ (r · ∇ f ) y(r) dr. 2D So, for an isotropic fluid, we have, for n ≥ 2,  1 Bn(v) = r f ′(r) yn−2(r) dr. 2D

(14)

(15)

For the hard-sphere fluid, f ′(r) = δ(r − 1), which is the δfunction peaked at r = 1, and Eq. (15) is reduced to Bn(v) = B2 yn−2(1),

 where B2 = r f ′(r) · 1 · dr/(2D) = π D/2/[D Γ(D/2)] is half of the volume of the unit sphere in D dimensions, and Γ(z)

214110-3

Zhang, Lai, and Pettitt

J. Chem. Phys. 142, 214110 (2015)

is the gamma function. Since the closure is approximate, Eqs. (13) and (15) generally yield different values for n ≥ 4. In addition to Eqs. (13) and (15), which apply to a general integral equation, there are two formulae specific to the PY equation. Baxter’s pressure formula48 can be written as  ρt(0) + ρ2c(0) ˜ dk (c) − βP = log[1 − ρ c(k)] ˜ (2π) D 2    [1 + t(r)]2(r · ∇ f ) 2 = ρ− ρ c(r) + dr, (16) 2D where we have used the virial theorem3,49 in the second step (see supplementary material 262). The second expression of Eq. (16) is valid even if ρ c(k) ˜ exceeds 1.0 for some k. For the hard-sphere fluid, it can be further simplified as − βP(c) = ρ[c(0) ˜ + c2(1−)ρB2] = ρ[c(0) ˜ + g 2(1+)ρB2], Taking the nth component of the density expansion yields a different expression for Bn . The resulting value, however, is the same as Bn(c).48 The second derivative of pressure can be written as (see supplementary material 262)  − β χex ≡ −ρ2∂ρ2( βPex) ≈ ρ2 [c(r) + h(r) t(r)] dr. Thus, Bn(PY, χ)

1 ≈ n(1 − n)

   n−2  c (r) + hn−2− j (r) t j (r) dr.  n−2  j=1  (17)

This value, referred to as the χ-route virial coefficient below, differs from Bn(c) and Bn(v). Since the above virial coefficients are computed from numerical integration, their precision can be improved by an extrapolation on the grid size (see supplementary material 1(B)62). C. HNC equation

We now turn to the HNC case. The OZ relation, Eq. (9), is still applicable. Below, we determine the density expansion of the HNC cavity distribution function, y(r) = exp t(r), y(r) =

+∞ 

yl (r) ρl ,

(18)

l=0

with y0(r) = 1, such that the components yl (r) can be used in Eq. (10). Previously, Shaul et al. proposed the following formula:14 l   1 t l (r) · · · t l k (r), yl (r) = k ! (l , ..., l ) 1 k=0 1 k l 1+···+l k =l

where the inner sum is carried over all possible k-tuples of positive integers (l 1, . . . , l k ), whose sum is l. The formula, however, becomes difficult to use for a high order, l, because of the exponentially growing number of the k-tuples. Below, we propose an algorithm with polynomial scaling of time.

First, differentiating Eq. (6) with respect to ρ yields ∂ρ y(r) = y(r) [∂ρt(r)]. Then, expanding both sides using Eqs. (8) and (18) gives, for l ≥ 1, yl (r) =

l 

yl− j (r) j t j (r)/l.

(19)

j=1

In Eq. (19), the yl− j (r) on the right-hand side have orders lower than l. So we can readily compute yl (r) using an increasing order of l [with y0(r) = 1]. For the first few orders, we have [dropping the “(r)”], y1 = y0 t 1 = t 1, y2 = ( y1 t 1 + 2 y0 t 2)/2 = 21 t 12 + t 2, and y3 = ( y2 t 1 + 2 y1t 2 + 3 y0 t 3)/3 = 61 t 13 + t 1t 2 + t 3, . . .. Once yl (r) is known, we can compute cl (r) from Eq. (10). Applying Eqs. (9), (10), and (19) in succession as F.T.

Eq. (9)

F.T.

cl−1(r) −−−→ c˜l−1(k) −−−−−→ t˜l (k) −−−→ Eq. (19)

Eq. (10)

t l (r) −−−−−−→ yl (r) −−−−−−→ cl (r), yields the correlation functions of any order l in polynomial time. We can generalize the above technique to a larger class of closures y(r) = Y [t (r)]32–35,39–45,47 (supplementary material 1(C)62), and to the Yvon-Born-Green (YBG)2,4,5,18,50 and Kirkwood2,4,51,52 integral equations (see supplementary material 1((D) and (E)), respectively62). 1. Virial coefficients

The compressibility- and virial-route virial coefficients, B(c) and B(v), can be still computed from Eqs. (13) and (15), respectively. Additionally, the virial coefficient can be derived from the HNC-specific expressions of the chemical potential, µ, and the free energy26,27,53–55 (see also supplementary material 262). For example, using the chemical potential formula with ∂ρ P = ρ ∂ρ µ yields    n−2 1−n 1  (HNC, µ) cn−2(r) − Bn = hn−2− j (r) t j (r) dr. n 2   j=1 The value, however, is the same as Bn(v)27,53 because of the virial theorem3,49 (see supplementary material 2 for a proof62). In terms of the graphical expansion,3,5,55 the integrands of Bn(HNC, µ) and Bn(v) differ only by the choices of the root vertices. For the hard-sphere fluid, we have an additional cavityroute virial coefficient52 (see also supplementary material 262), n−1 [log y(0)]n−1, (20) n where [log y(0)]n−1, the coefficient of ρn−1 in the density expansion of log y(0), can be computed from the components, yl (0) (l = 0, . . . , n − 1). This relation, however, demands a highly accurate y(0), which makes it applicable to very few cases, such as the HNC equation in high dimensions. Bn(y) =

D. Detailed self-consistent (DSC) equation

Neither the PY nor the HNC equation yields the same virial coefficients from the compressibility and virial routes.

214110-4

Zhang, Lai, and Pettitt

J. Chem. Phys. 142, 214110 (2015)

Here, we consider a variant on interpolation between the PY and HNC equations to force the consistency32–35,38–40,42,44 of all virial coefficients. We assume a form for the new cavity distribution function to be y(r) = y (PY)(r) + correction +∞  = 1 + t(r) + λ l wl (r) ρl ,

Bn(c) = (1 + κ n ) Bn(v),

l=2

where wl (r) is the lth order correction function, and the parameter λ l is to be determined. Particularly, we set wl (r) =

yl(HNC)(r)

−

yl(PY)(r)

correct. If we assume that the wl (r) given by Eq. (21) deviates from the correct functional form by a small amount ε, then the virial coefficients from the two routes should also differ by an amount proportional to ε, to the leading order. Instead of a strict self-consistency, we now set

with the parameter κ n of order ε. We recover Eq. (23) if κ n = 0. Equation (24) becomes λ n−2 = −

= [exp t(r)] l − t l (r), (21)

with [exp t(r)]l denoting the coefficient of ρl in the density expansion of exp t(r). Thus, w2(r) = [t 1(r)]2/2 = [( f ∗ f )(r)]2/2 and w3(r) = t 1(r) t 2(r) + [t 1(r)]3/6. Below, we determine the parameter λ l to force the equality of Eqs. (13) and (15) for each order l. From Eq. (4), we have c(r) = f (r) [1 + t(r)] + [1 + f (r)]

+∞ 

λ l wl (r)ρl .

l=2

(25)

Bn(c, uncorr) − (1 + κ n ) Bn(v, uncorr) ∆Bn(c) − (1 + κ n ) ∆Bn(v)

.

Now that Bn(c) , Bn(v), we shall use Bn(c) (which is usually more accurate than Bn(v), as shown in Sec. III A) as the best estimate, Bn(c) = Bn(c, uncorr) + λ n−2 ∆Bn(c). For simplicity, we will focus on a constant κ n that is turned on for n ≥ n0. Once the κ n and n0 are tuned for low-order virial coefficients, e.g., by matching the known values, we anticipate them to also improve higher-order virial coefficients.

So cl (r) = f (r) δl0 + f (r) t l (r) + λ l [1 + f (r)] wl (r). (22) Now, by demanding

set Bn(c) = Bn(v),

where

B(c, uncorr) − n ∆Bn(c)

− −

Bn(v, uncorr) , ∆Bn(v)

We now give an alternative to the DSC equation. Here, we

(23)

and using Eqs. (13), (15), and (22), we get, for n ≥ 4, λ n−2 =

2. Density-dependent parameter and λ-DSC equation

(24)

 1 ≡− f (r) t n−2(r) dr, n 1 Bn(v, uncorr) ≡ r · ∇ f (r) t n−2(r) dr, 2D  1 ∆Bn(c) ≡ − [ f (r) + 1] w n−2(r) dr, n Bn(c, uncorr)

and

 1 ≡ r · ∇ f (r) w n−2(r) dr. 2D The cases of λ l ≡ 0 and λ l ≡ 1 correspond to the PY and HNC closures, respectively. Note that Eq. (23) differs from the commonly used selfconsistent condition,32–35,39–47 which demands the same compressibility or pressure from Eqs. (12) and (14) at a certain density, in that the former is a set of order-by-order conditions instead of a single condition. For this reason, we refer to Eq. (23) as the DSC condition. Some generalizations are discussed in supplementary material 1((F) and (G)).62

y(r) = 1 + t(r) + λ [exp t(r) − 1 − t(r)], and expand the parameter λ as a power series of the density, ρ,33,37,39,42 λ = λ 2 + λ 3 ρ + λ 4 ρ2 + · · ·.

(26)

The coefficients λ l are to be determined to satisfy Eq. (23) following a procedure in supplementary material 1(F).62 We refer to this variant as the λ-DSC equation below, to distinguish it from the integral equation based on Eq. (21), which is simply referred to as the DSC equation. The λ-DSC equation is superior for low orders in low dimensions, whereas the DSC equation based on Eq. (21) is preferred for higher orders or in higher dimensions (supplementary material 1(F)62).

∆Bn(v)

1. Empirical correction

The DSC equation is still approximate, and below is an empirical correction. The DSC condition, Eq. (23), would yield the exact result if the functional form of wl (r) were

TABLE I. Integral equations based on the OZ relation. Closure PY25 HNC26–28 Hurst32 Rowlinson 133 Rowlinson 234 HC39,40 MS43 BPGG45 Verlet41 MP47 a RY44 Quadratic aη

y(r) = Y [t(r)] y = 1+t y = expt t = y m log y t = η log y + (1 −η) (y − 1) y = (1 −η)(1 + t) +η expt y = (1 + s t)1/s √ y = exp( 1 + 2t − 1) y = exp[(1 + st)1/s − 1] y = exp[t − At 2/(2 + Bt)] y = 1 + [exp(η t) − 1]/η y = 1 + [exp( f swt) − 1]/ f sw, with f sw = 1 − exp(−α r ) y = 1 + t +ηt 2/2

here is equivalent to 1/(a + 1) in the original paper.47

Integral equation

B5/B24

B6/B25

B7/B26

B 8/B 27

B 9/B 28

B10/B29

B11/B 210

B 12/B211

c v y

0.342 41 0.225 22 0.337 84

0.133 50 0.047 47 0.075 95

0.031 852 −0.006 357 −0.010 595

−0.021 721 −0.007 000 −0.012 001

−2.874 × 10−2 −4.118 × 10−3 −7.207 × 10−3

−2.907 × 10−2 −1.206 × 10−4 −2.143 × 10−4

−1.540 × 10−2 −6.823 × 10−4 −1.228 × 10−3

−1.074 × 10−2 5.239 × 10−4 9.526 × 10−4

−6.654 × 10−4 −5.123 × 10−4 −9.393 × 10−4

Kirkwood2,4,51,52

c v y

0.441 82 0.139 96 0.209 19

0.282 54 0.065 56 0.111 44

0.220 338 0.026 005 0.061 075

0.145 991 0.029 434 0.050 743

1.305 × 10−1 9.369 × 10−3 2.779 × 10−2

8.310 × 10−2 1.670 ×10−2 2.673 ×10−2

8.835 × 10−2 9.215 × 10−4 1.104 × 10−2

4.711 × 10−2 1.175 × 10−2 1.575 × 10−2

6.729 × 10−2 −3.971 × 10−3 2.345 × 10−3

PY25

c v y

0.296 88 0.250 00 0.078 13

0.121 09 0.085 94 −0.021 88

0.044 922 0.027 344 0.011 068

0.015 625 0.008 301 −0.004 325

5.188 × 10−3 2.441 × 10−3 1.984 × 10−3

1.663 × 10−3 7.019 × 10−4 −8.613 × 10−4

5.188 × 10−4 1.984 × 10−4 3.922 × 10−4

1.583 × 10−4 5.531 × 10−5 −1.772 × 10−4

4.745 × 10−5 1.526 × 10−5 8.146 × 10−5

HNC26–28

c v y

0.209 19 0.445 31 0.890 63

0.049 27 0.144 75 0.367 44

0.028 067 0.038 180 0.154 134

0.006 113 0.044 086 0.090 725

−1.459 × 10−3 6.461 × 10−3 1.077 × 10−2

4.033 × 10−3 −8.554 × 10−3 3.242 × 10−3

−1.010 × 10−3 1.127 × 10−2 1.800 × 10−2

−1.563 × 10−3 −1.949 × 10−3 −1.030 × 10−2

2.209 × 10−3 −7.102 × 10−3 −3.095 × 10−3

Hurst32 m = 0.417 18

c v

0.282 35 0.282 35

0.110 19 0.091 51

0.038 585 0.035 220

0.012 813 0.009 847

4.309 × 10−3 3.271 × 10−3

1.219 × 10−3 1.304 × 10−3

4.107 × 10−4 −1.341 × 10−4

8.969 × 10−5 2.582 × 10−4

2.133 × 10−5 −6.910 × 10−5

Rowlinson 133 η = 0.165 64

c v

0.282 35 0.282 35

0.113 24 0.078 81

0.041 426 0.031 550

0.014 013 0.008 704

4.942 × 10−3 1.856 × 10−3

1.432 × 10−3 1.610 × 10−3

4.915 × 10−4 −4.017 × 10−4

1.432 × 10−4 2.350 × 10−4

1.596 × 10−5 3.024 × 10−5

Rowlinson 234 η = 0.165 64

c v

0.282 35 0.282 35

0.109 20 0.095 68

0.036 963 0.041 070

0.011 460 0.014 171

3.466 × 10−3 5.922 × 10−3

7.813 × 10−4 2.681 × 10−3

2.042 × 10−4 5.186 × 10−4

9.288 × 10−7 5.440 × 10−4

−1.031 × 10−5 4.024 × 10−5

HC39,40 s = 0.834 36

c v

0.282 35 0.282 35

0.111 89 0.084 43

0.040 157 0.033 256

0.013 485 0.009 078

4.659 × 10−3 2.485 × 10−3

1.334 × 10−3 1.457 × 10−3

4.576 × 10−4 −3.026 × 10−4

1.177 × 10−4 2.630 × 10−4

1.877 × 10−5 −2.975 × 10−5

MS43

c v

0.296 88 0.250 00

0.111 35 0.126 63

0.038 090 0.023 846

0.011 537 0.021 918

4.055 × 10−3 −3.723 × 10−3

8.826 × 10−4 7.791 × 10−3

4.929 × 10−4 −5.373 × 10−3

−7.903 × 10−5 5.132 × 10−3

1.525 × 10−4 −4.482 × 10−3

BPGG45 s = 1.834 36

c v

0.282 35 0.282 35

0.103 76 0.118 38

0.033 542 0.036 488

0.010 794 0.014 865

3.426 × 10−3 3.981 × 10−3

9.949 × 10−4 1.749 × 10−3

2.638 × 10−4 4.766 × 10−5

4.516 × 10−5 2.452 × 10−4

2.045 × 10−5 −7.377 × 10−5

Verlet41 A = 0.834 36 B = 1.547 51

c v

0.282 35 0.282 35

0.106 59 0.106 59

0.034 177 0.046 552

0.010 100 0.015 639

2.993 × 10−3 5.844 × 10−3

7.601 × 10−4 1.814 × 10−3

2.436 × 10−4 2.341 × 10−4

2.875 × 10−5 2.884 × 10−4

1.430 × 10−5 −5.269 × 10−5

MP47 η = 0.165 64

c v

0.282 35 0.282 35

0.110 54 0.090 05

0.038 712 0.036 136

0.012 713 0.010 472

4.243 × 10−3 3.442 × 10−3

1.168 × 10−3 1.559 × 10−3

3.886 × 10−4 −1.293 × 10−4

8.324 × 10−5 2.951 × 10−4

1.508 × 10−5 −5.052 × 10−5

RY44 α = 0.170 15

c v

0.280 56 0.280 56

0.111 01 0.086 31

0.038 375 0.037 204

0.012 714 0.009 674

4.474 × 10−3 3.444 × 10−3

1.002 × 10−3 1.951 × 10−3

4.974 × 10−4 −7.341 × 10−4

2.546 × 10−5 6.736 × 10−4

1.021 × 10−6 −2.314 × 10−4

Quadratic η = 0.165 64

c v

0.282 35 0.282 35

0.110 81 0.088 94

0.039 034 0.035 331

0.012 917 0.009 974

4.359 × 10−3 3.135 × 10−3

1.220 × 10−3 1.466 × 10−3

4.118 × 10−4 −1.901 × 10−4

9.377 × 10−5 2.780 × 10−4

1.760 × 10−5 −5.299 × 10−5

0.286 950

0.110 251

0.038 882 1

0.013 022 9

4.1825 × 10−3

1.3097 × 10−3

4.030 × 10−4

1.206 × 10−4

3.77(32) × 10−5

Referencea

a The reference values were from, Ref. 13 except for

B12, which has been improved from the previous study.15 A few other improved or new values were 1.64(18) × 10−3 [214], 1.241(19) × 10−4 [311], 1.01(8) × 10−5 [410], 1.01(37) × 10−6 [411], 3.13(10) × 10−5 [511], and −2.34(20) × 10−5 [512] (the numbers in “[. . . ]” assume the format of D n ).

J. Chem. Phys. 142, 214110 (2015)

YBG2,4,5,18,50

Zhang, Lai, and Pettitt

B4/B23

214110-5

TABLE II. Virial coefficients of the three-dimensional hard-sphere fluid. In each cell, the numbers on the first, second, and third (if exists) rows are from the compressibility (c), virial (v), and cavity (y) routes, respectively. The cavity-route numbers were generally inaccurate, so we only list them for the first four equations. All integral equations produced the correct B3/B22 = 5/8 from the compressibility and virial routes, but not necessarily from the cavity route. The cavity-route values were 5/6 (YBG), 5/8 (Kirkwood), −1/12 (PY, MS, RY), 5/4 (HNC), and 0.137 52 (the rest).

214110-6

Zhang, Lai, and Pettitt

J. Chem. Phys. 142, 214110 (2015)

FIG. 1. Virial coefficients of the hard-sphere fluid from the DSC, PY, and HNC equations compared to those from Mayer sampling (reference) in (a) D = 2, (b) D = 7, (c) D = 10, and (d) D = 15 dimensions. Empty (thin) and filled (thick) symbols indicate positive and negative numbers, respectively. In (e) and (f), the (ref) relative error ε(B n ) ≡ B n /B n − 1 from the DSC equation (diamonds) was compared with that from the B4-SC equation (circles). For D = 10 [panel (e)], the errors of the two appeared to be similar in high orders as they became comparable to that of the Mayer sampling results.

III. RESULTS

In Sec. III A, the method described in Secs. II B and II C was used to compute the virial coefficients from many existing integral equations. We further show in Secs. III B and III C that the DSC equation described in Sec. II D yielded reasonably accurate virial coefficients for the hard-sphere fluid

and Gaussian model in high dimensions, although it is less successful in low dimensions. A. Comparison of integral equations

To show the generality of the method, we computed the virial coefficients from a few integral equations based on the

214110-7

Zhang, Lai, and Pettitt

OZ relation, Eq. (2), as listed in Table I, as well as those from the YBG2,4,5,18,50 and Kirkwood2,4,51,52 equations, using the techniques in supplementary material 1((C)–(E)).62 For an integral equation with a free parameter, one usually adjusts the parameter to achieve the consistency of Eqs. (12) and (14) at a certain density. Doing so for the low density limit yields the same B4 (and B5 for the Verlet equation41) from the compressibility and virial routes. The resulting parameters, as listed in Table II, may differ slightly from the literature values chosen for the consistency at another density. Note that the above consistency condition, referred to as the B4-SC condition below, is limited to the lowest order only, and thus, higher-order virial coefficients from the two routes would still differ. The virial coefficients of the integral equations are shown in Table II for the three-dimensional hard-sphere fluid after the above parameter adjustment. The low-order results agreed with the literature values. Except for the YBG, Kirkwood,

J. Chem. Phys. 142, 214110 (2015)

and HNC cases, most equations produced reasonably accurate compressibility-route virial coefficients, B(c), up to about the tenth order, although B11 and B12 were less satisfactory. This shows the importance of higher-order virial coefficients in differentiating integral equations.8 For higher orders, B(c) were more accurate than the virial-route results, B(v) [this prompted us to use the former in the modified DSC condition, Eq. (25)]. The cavity-route results, B(y), however, were usually much less accurate than B(c) and B(v). B. Virial coefficients of the hard-sphere fluid

We computed the virial coefficients of the hard-sphere fluid in various dimensions using the DSC and related equations in Sec. II D. The results are shown in Figs. 1–3, as well as Table III. The Mayer-sampling results13,15 were used as the references. Generally, the DSC equation worked the best in high dimensions, but less so in lower dimensions.

FIG. 2. Virial coefficients of the hard-sphere fluid from the DSC and λ-DSC integral equations (“ ,” “△,” “D,” and “×;” black dashed lines for the DSC equation with κ n = 0, blue dotted lines for that with κ n , 0, green dotted-dashed lines for the λ-DSC equation) compared to those from Mayer sampling (reference, “,” “▽,” “ ⃝,” and “+;” red solid lines). (a) D = 2, 3, and 4; (b) D = 5, 6, 7, and 8; (c) D = 9, 10, and 11; (d) D = 12, 13, and 14. Empty (thin) and filled (thick) symbols indicate positive and negative numbers, respectively.

214110-8

Zhang, Lai, and Pettitt

J. Chem. Phys. 142, 214110 (2015)

FIG. 3. Virial coefficients of the hardsphere fluid for dimensions D = 15 to 30 from the DSC integral equation with κ n ≡ 0 (“ ” and “△”) compared to those from Mayer sampling (reference, “” and “▽”). Empty (thin) and filled (thick) symbols indicate positive and negative numbers, respectively.

1. DSC equation

We first compared the DSC equation under Eq. (23) with the PY and HNC ones, as the former is an interpolation of the latter two. As shown in Fig. 1, the DSC equation was superior to the PY and HNC equations in most cases. The PY equation gave reasonably accurate Bn(c) and Bn(v) in D = 2, but not so in the higher dimensions. The Bn(PY, χ) from Eq. (17) was inferior in D = 2, although its error approached those of Bn(c) and Bn(v) in the higher dimensions. We omit the cavity-route Bn(y) from Eq. (20) in the PY case due to the poor accuracy. For the HNC equation, although Bn(c), Bn(v), and Bn(y) were all erratic in low dimensions, they improved quickly with D, with Bn(c) being the best among the three. The above comparison shows that the PY and HNC equations worked best for low and high dimensions, respectively. Thus, we expect the DSC equation to be closer to the PY one in low dimensions, but to the HNC one in high dimensions. This was indeed the case. We found that the λ n−2 in Eq. (24) increased with the dimension and order. In very high dimensions, λ n−2 was close to 1.0, the HNC limit.

As shown in Figs. 1(e) and 1(f), the order-by-order DSC equation based on Eq. (23) was more accurate than the singleparameter B4-SC Rowlinson 2 equation (cf. Table I) in high dimensions (note that the former is reduced to the latter, if λ n−2 ≡ η). 2. Improving the DSC equation

The DSC equation can be improved in two ways. Consider the modified DSC condition, Eq. (25), versus the strict one, Eq. (23). The latter is the special case of κ n ≡ 0. In high dimensions, it sufficed to use Eq. (23), as shown in Table III and Figs. 2 and 3, and the accuracy improved with increasing D. For D = 2, 3, or 6, more accurate virial coefficients can be achieved by Eq. (25) with nonzero κ n , as shown in Figs. 2(a) and 2(b). For D = 4 or 5, however, the modified DSC equation failed. Alternatively, in low dimensions, one can switch to the λDSC equation based on the density series of λ, Eq. (26). For D = 5, the λ-DSC equation yielded the correct signs of the first 12 virial coefficients, as shown in Fig. 2(b). For D = 6, only

214110-9

Zhang, Lai, and Pettitt

J. Chem. Phys. 142, 214110 (2015)

TABLE III. Virial coefficients of the hard-sphere fluid in high dimensions. In each cell, the upper and lower numbers were obtained from the DSC integral equation, Eq. (23), and Mayer sampling15 [algorithm G3, with the parameter M = 1 (or M = 2 for D = 30)], respectively. The former number was extrapolated to zero grid size (supplementary material 1(B)62). We only keep a maximum of 10 significant digits. B64/B263

B96/B295

D

B 32/B231

10

−1.108 603 771(2) −1.16(6)

−5.810 334 42(2) × 105

−1.144 362 730(5) × 1012

−3.919 593 67(2) × 1018

11

−1.310 853 792 −1.270(27)

−1.546 364 657 × 106

−7.854 488 089 × 1012

−7.370 834 712 × 1019

12

−1.119 249 796(2) −1.135 1(43)

−2.247 090 07(1) × 106 −2.36(13) × 106

−2.246 234 59(2) × 1013

−4.414 497 57(5) × 1020

13

−0.752 218 329 1 −0.756 9(19)

−2.105 998 963 × 106 −2.195(27) × 106

−3.410 012 696 × 1013

−1.155 110 104 × 1021

14

−0.422 622 529(1) −0.426 6(8)

−1.426 782 124(8) × 106 −1.445(13) × 106

−3.239 990 574(3) × 1013 −3.19(27) × 1013

−1.636 649 41(2) × 1021

15

−0.207 220 569 9 −0.208 50(8)

−7.568 079 484(2) × 105 −7.679(25) × 105

−2.161 072 905(3) × 1013 −2.176(47) × 1013

−1.458 667 826(6) × 1021 −1.77(15) × 1021

16

−0.091 491 486 1(3) −0.092 12(17)

−3.327 636 25(3) × 105 −3.358(13) × 105

−1.098 618 930(7) × 1013 −1.122(12) × 1013

−9.104 481 7(5) × 1020 −9.54(33) × 1020

17

−0.037 226 128 39 −0.037 331(11)

−1.264 478 996 × 105 −1.273 1(19) × 105

−4.520 299 46(1) × 1012 −4.594(15) × 1012

−4.305 514 16(4) × 1020 −4.347(46) × 1020

18

−0.014 202 636 7(7) −0.014 206(13)

−4.283 365 85(3) × 104 −4.309(6) × 104

−1.574 349 52(9) × 1012 −1.588(6) × 1012

−1.636 037 7(2) × 1020 −1.669(14) × 1020

19

−5.148 274 944 × 10−3 −5.155 4(10) × 10−3

−1.324 198 556(2) × 104 −1.328 0(10) × 104

−4.801 958 36(3) × 1011 −4.837(9) × 1011

−5.223 260 65(9) × 1019 −5.263(23) × 1019

20

−1.791 101 61(1) × 10−3 −1.792 3(7) × 10−3

−3.804 058 6(2) × 103 −3.821 2(35) × 103

−1.316 684 1(1) × 1011 −1.322 9(26) × 1011

−1.450 072 1(2) × 1019 −1.467(7) × 1019

21

−6.027 766 825(1) × 10−4 −6.038 3(34) × 10−4

−1.029 817 070(2) × 103 −1.031 6(9) × 103

−3.312 378 87(3) × 1010 −3.320 9(48) × 1010

−3.595 503 4(1) × 1018 −3.620(12) × 1018

22

−1.974 462 77(2) × 10−4 −1.975 8(10) × 10−4

−265.633 90(2) −265.90(17)

−7.768 860(1) × 109 −7.786(12) × 109

−8.132 020(2) × 1017 −8.184(20) × 1017

23

−6.325 876 896(2) × 10−5 −6.326 0(30) × 10−5

−65.860 111 1(3) −65.83(5)

−1.720 539 69(3) × 109 −1.723 5(20) × 109

−1.705 974 11(9) × 1017 −1.711 3(27) × 1017

24

−1.990 073 84(6) × 10−5 −1.990 6(6) × 10−5

−15.806 042(1) −15.823(11)

−3.634 870(1) × 108 −3.639 8(27) × 108

−3.364 252(2) × 1016 −3.375(9) × 1016

25

−6.166 835 490(2) × 10−6 −6.169 9(8) × 10−6

−3.692 688 13(2) −3.691 0(22)

−7.385 874 2(2) × 107 −7.380(8) × 107

−6.304 513 4(6) × 1015 −6.295(10) × 1015

26

−1.887 167 46(8) × 10−6 −1.887 8(6) × 10−6

−0.843 673 86(7) −0.843 7(5)

−1.453 087 1(4) × 107 −1.453 5(10) × 107

−1.132 401 5(1) × 1015 −1.133 5(17) × 1015

27

−5.715 054 215(4) × 10−7 −5.717 6(13) × 10−7

−0.189 211 501(2) −0.189 11(9)

−2.783 166 2(1) × 106 −2.785 7(18) × 106

−1.964 150 5(3) × 1014 −1.963 7(24) × 1014

28

−1.715 682 01(9) × 10−7 −1.715 37(25) × 10−7

−0.041 782 545(4) −0.041 795(18)

−5.212 393(3) × 105 −5.217 9(20) × 105

−3.306 362 0(4) × 1013 −3.305 9(37) × 1013

29

−5.112 987 071 × 10−8 −5.114 5(12) × 10−8

−9.107 961 141(3) × 10−3 −9.106 1(24) × 10−3

−9.582 179 92(3) × 104 −9.581 5(43) × 104

−5.436 905(1) × 1012 −5.438 4(43) × 1012

30

−1.514 410 5(1) × 10−8 −1.514 72(17) × 10−8

−1.963 941 0(7) × 10−3 −1.963 6(6) × 10−3

−1.733 619 45(5) × 104 −1.734 5(12) × 104

−8.743 38(2) × 1011 −8.746(11) × 1011

the λ-DSC equation produced the qualitatively correct density components of the direct correlation function, cl (r), as shown in Fig. 4. Unfortunately, the λ-DSC equation can be less stable than the DSC equation for higher orders or higher dimensions, as discussed in supplementary material 1(F).62 In future work, we may consider interpolations of the two equations for better accuracy.

B128/B2127

C. Virial coefficients of the Gaussian model

We applied the DSC equation to the Gaussian model.3,56,57 The model has an artificial potential such that the Mayer f function is given by f (r) = − exp(−r2). The main appeal of the model is that it offers many exactly computable virial coefficients.3,6,7,57 Such a calculation involves an enumeration of all

214110-10

Zhang, Lai, and Pettitt

J. Chem. Phys. 142, 214110 (2015)

FIG. 4. Density components c l (r ) of the direct correlation function of the sixth-dimensional hard-sphere fluid from the DSC, λ-DSC, PY, and HNC integral equations compared to those from Mayer sampling (reference). Note that c l (r ) is discontinuous at r = 1: c l (r ) = −t l (r ) for r < 1 and c l (r ) = yl (r ) − t l (r ) for r > 1. The r < 1 results and their continuous extensions to r > 1 (only for the Mayer sampling results, light red lines) use the left vertical scale; the r > 1 results and their continuous extensions to r < 1 (light red lines) use the right vertical scale.

TABLE IV. Virial coefficients of the Gaussian model. For D ≥ 6, each cell contains, in addition to the exact virial coefficient (the last number, 10 significant digits), results from the DSC integral equation with zero or nonzero κ n . D

B 9/B28

B10/B 29

B11/B210

κ na

B12/B211

Error

1 2 3 4 5

+0.239 049 430 7 −0.050 424 322 27 +0.015 197 607 20 −0.010 195 771 05 +0.022 768 548 66

+0.301 364 071 8 +0.017 748 037 02 −0.013 849 654 13 +8.798 988 664 × 10−3 −0.021 415 050 98

+0.098 334 123 20 +0.039 204 109 10 +1.334 755 917 ×10−3 −5.576 767 192 × 10−3 +0.021 195 589 66

−0.229 341 236 5 −0.013 244 042 30 +7.604 327 249 × 10−3 +2.318 690 261 × 10−3 −0.022 023 012 29

6

+0.031 0 +0.032 2 +0.031 815 511 25

−0.036 6 −0.037 80 −0.037 875 081 35

+0.045 1 +0.046 3 +0.046 983 841 60

−0.057 6 −0.058 8 −0.060 326 097 09

0 (0.047)n≥8

4.6% 2.6% (Exact)

7

+0.020 83 +0.020 982 +0.020 973 693 43

−0.026 25 −0.026 44 −0.026 483 325 18

+0.034 5 +0.034 71 +0.034 835 667 14

−0.046 86 −0.047 16 −0.047 394 686 09

0 (0.008 5)n≥8

1.12% 0.50% (Exact)

8

+0.010 165 +0.010 196 +0.010 190 170 60

−0.012 85 −0.012 883 −0.012 882 810 25

+0.016 94 +0.016 991 +0.016 998 429 37

−0.023 15 −0.023 215 −0.023 232 262 53

0 (0.003 25)n≥8

0.35% 0.13% (Exact)

9

+4.219 × 10−3 +4.222 5 × 10−3 +4.222 773 149 × 10−3

−5.197 × 10−3 −5.201 2 × 10−3 −5.202 299 404 × 10−3

+6.701 × 10−3 +6.707 × 10−3 +6.708 746 298 × 10−3

−8.973 × 10−3 −8.980 × 10−3 −8.983 225 704 × 10−3

0 (0.000 924)n≥7

0.119% 0.036% (Exact)

10

+1.599 1 × 10−3 +1.599 49 × 10−3 +1.599 635 591 × 10−3

−1.893 5 × 10−3 −1.894 0 × 10−3 −1.894 245 463 × 10−3

+2.355 7 × 10−3 +2.356 3 × 10−3 +2.356 666 276 × 10−3

−3.052 6 × 10−3 −3.053 4 × 10−3 −3.053 853 001 × 10−3

0 (0.000 29)n≥6

0.041% 0.014% (Exact)

a The

(Exact) (Exact) (Exact) (Exact) (Exact)

( ) (ref) (ref) nonzero κ n were optimized to minimize the maximal relative error of the first 12 virial coefficients: max max{log(B n /B n ), B n /B n − 1} . n≤12

214110-11

Zhang, Lai, and Pettitt

J. Chem. Phys. 142, 214110 (2015)

TABLE V. Higher-order virial coefficients of the Gaussian model. Each cell contains the virial coefficient from the DSC integral equation with κ n = 0 (first row), or that with the nonzero κ n , as listed in Table IV (second row), and the exact value (the third row, obtained from summing biconnected graphs with fewest numbers of edges, followed by an extrapolation along the number of edges assuming a geometric series). Dn

B n /B2n−1

Dn

B n /B 2n−1

Dn

B n /B2n−1

Dn

B n /B2n−1

813

+0.032 58 +0.032 67 +0.032 700 584 60

814

−0.047 02 −0.047 14 −0.047 190 578

815

+0.069 3 +0.069 50 +0.069 574 8

816

−0.104 1 −0.104 39 −0.104 499

913

+0.012 399 +0.012 409 +0.012 413 892 18

914

−0.017 60 −0.017 613 −0.017 619 937 78

915

+0.025 56 +0.025 582 +0.025 591 147 5

916

−0.037 88 −0.037 906 −0.037 918 635

+0.057 11 +0.057 161 +0.057 177 34

918

−0.087 47 −0.087 54 −0.087 562 90

919

+0.135 82 +0.135 93 +0.135 957 7

920

−0.213 5 −0.213 69 −0.213 724

1013

+4.093 1 × 10−3 +4.094 3 × 10−3 +4.094 811 210 × 10−3

1014

−5.650 × 10−3 −5.651 8 × 10−3 −5.652 506 747 × 10−3

1015

+7.997 × 10−3 +7.999 6 × 10−3 +8.000 503 367 × 10−3

1016

−0.011 569 −0.011 572 1 −0.011 573 166 09

1017

+0.017 059 +0.017 063 3 +0.017 064 514 4

1018

−0.025 583 −0.025 589 8 −0.025 591 150

1019

+0.038 950 +0.038 960 7 +0.038 962 13

1020

−0.060 11 −0.060 127 2 −0.060 128 3

917

cluster integrals that can be represented graphically by biconnected diagrams.3,5 The number of these diagrams, however, grows rapidly with the order, n. We enumerated biconnected diagrams using the program N,58 and obtained the exact virial coefficients up to the 12th order (the previous literature57 went up to n = 9). The exact values were used below as the references. We computed the virial coefficients of the Gaussian model using the DSC equation. Particularly, we wish to see whether the modified DSC condition, Eq. (25) with a nonzero κ n , improved over the κ n = 0 case [Eq. (23)]. In this model, κ 4 ≡ 0 is assumed, because it happens to yield the correct functional form of w2(r) ∝ exp(−r 2) and hence the correct B4. As shown in Table IV, a nonzero κ n , obtained from minimizing the largest deviation from the exact B5, . . . , B12, reduced the error considerably. The same parameter also improved the accuracy of the next few coefficients in comparison to the κ n = 0 case for D ≥ 8, as shown in Table V. For D ≤ 5, however, the method failed.

IV. CONCLUSIONS

In summary, we have shown an algorithm with polynomial time scaling (which differs from that by Shaul et al.14) to compute the virial coefficients from various integral equations. Using series transformations on an order-by-order basis with the DSC criteria for a constraint, we accurately predicted the virial coefficients of the high-dimensional hardsphere fluid and Gaussian model. Despite the success, the method met some difficulties in low dimensions. In the future, the method may be improved by choosing better closures and/or better correction functions for the cavity distribution function. Alternatively, one might consider to incorporate other self-consistent conditions,31,52 or explicit bridge functions.14,30,59

ACKNOWLEDGMENTS

It is a pleasure to thank Dr. Y. Mei, Dr. K. Dyer, Dr. J. Perkyns, and Dr. G. Lynch for helpful discussions. The authors gratefully acknowledge the financial support of the National Science Foundation (No. CHE-1152876) and the Robert A. Welch Foundation (No. H-0037). A portion of this work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. OCI-1053575. In particular, the calculations were performed on Stampede60 at the Texas Advanced Computing Center at the University of Texas at Austin and on SuperMIC61 at the High Performance Computing center at Louisiana State University. 1J.

E. Mayer and M. G. Mayer, Statistical Mechanics (John Wiley & Sons, New York, 1940). 2T. L. Hill, Statistical Mechanics (McGraw-Hill Book Company, New York, 1956). 3G. E. Uhlenbeck and G. W. Ford, in Studies in Statistical Mechanics, edited by J. de Boer and G. E. Uhlenbeck (North-Holland Publishing Company, Armsterdam, 1962), Vol. 1, p. 119. 4S. A. Rice and P. Gray, The Statistical Mechanics of Simple Liquids (Interscience Publishers, 1965). 5J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, 3rd ed. (Academic Press, Amsterdam, 2007). 6D. A. McQuarrie, Statistical Mechanics (Harper & Row, New York, 1976). 7L. E. Reichl, A Modern Course in Statistical Physics (University of Texas Press, Austin, 1980). 8R. J. Wheatley, Phys. Rev. Lett. 110, 200601 (2013). 9M. Adda-Bedia, E. Katzav, and D. Vella, J. Chem. Phys. 128, 184508 (2008); 129, 144506 (2008). 10K. M. Benjamin, J. K. Singh, A. J. Schultz, and D. A. Kofke, J. Phys. Chem. B 111, 11463 (2007); A. J. Schultz, N. S. Barlow, V. Chaudhary, and D. A. Kofke, Mol. Phys. 111, 535 (2013); A. J. Schultz and D. A. Kofke, ibid. 107, 2309 (2009); K. R. S. Shaul, A. J. Schultz, and D. A. Kofke, Mol. Simul. 36, 1282 (2010); J. K. Singh and D. A. Kofke, Phys. Rev. Lett. 92, 220601 (2004). 11N. Clisby and B. M. McCoy, J. Stat. Phys. 122, 15 (2006); 114, 1361 (2004). 12I. Lyberg, J. Stat. Phys. 119, 747 (2005). 13A. J. Schultz and D. A. Kofke, Phys. Rev. E 90, 023301 (2014).

214110-12 14K.

Zhang, Lai, and Pettitt

R. S. Shaul, A. J. Schultz, A. Perera, and D. A. Kofke, Mol. Phys. 109, 2395 (2011). 15C. Zhang and B. M. Pettitt, Mol. Phys. 112, 1427 (2014). 16M. Bishop, N. Clisby, and P. A. Whitlock, J. Chem. Phys. 128, 034506 (2008); N. Clisby and B. M. McCoy, Pramana 64, 775 (2005). 17J. D. van der Waals, Proc. Kon. Ak. Wet. 1, 138 (1899); J. J. van Laar, ibid. 1, 273 (1899); L. Boltzmann, ibid. 1, 398 (1899); G. S. Rushbrooke and P. Hutchinson, Physica 27, 647 (1961); J. S. Rowlinson, Mol. Phys. 7, 593 (1964); P. C. Hemmer, J. Chem. Phys. 42, 1116 (1965); S. Kim and D. Henderson, Phys. Lett. A 27, 378 (1968); J. E. Kilpatrick, in Advances in Chemical Physics, edited by I. Prigogine and S. A. Rice (John Wiley & Sons, Inc., 1971), Vol. 20, p. 39; K. M. Dyer, J. S. Perkyns, and B. M. Pettitt, Theor. Chem. Acc. 105, 244 (2001); T. Sun and A. S. Teja, J. Phys. Chem. 100, 17365 (1996). 18B. R. A. Nijboer and L. Van Hove, Phys. Rev. 85, 777 (1952). 19K. W. Kratky, J. Stat. Phys. 27, 533 (1982). 20M. Luban and A. Baram, J. Chem. Phys. 76, 3233 (1982); C. G. Joslin, ibid. 77, 2701 (1982); N. Clisby and B. M. McCoy, J. Stat. Phys. 114, 1343 (2004). 21A. J. Masters, J. Phys.: Condens. Matter 20, 283102 (2008). 22N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953); M. N. Rosenbluth and A. W. Rosenbluth, ibid. 22, 881 (1954); F. H. Ree and W. G. Hoover, ibid. 40, 939 (1964); 46, 4181 (1967); E. J. J. van Rensburg, J. Phys. A: Math. Gen. 26, 4805 (1993); M. Bishop, A. Masters, and J. H. R. Clarke, J. Chem. Phys. 110, 11449 (1999); A. Y. Vlasov, X. M. You, and A. J. Masters, Mol. Phys. 100, 3313 (2002); J. Kolafa, S. Labík, and A. Malijevský, Phys. Chem. Chem. Phys. 6, 2335 (2004); M. Bishop, A. Masters, and A. Y. Vlasov, J. Chem. Phys. 121, 6884 (2004); 122, 154502 (2005); S. Labík, J. Kolafa, and A. Malijevský, Phys. Rev. E 71, 021105 (2005). 23R. Finken, M. Schmidt, and H. Löwen, Phys. Rev. E 65, 016108 (2001); H. L. Frisch and J. K. Percus, ibid. 60, 2942 (1999); H. L. Frisch, N. Rivier, and D. Wyler, Phys. Rev. Lett. 55, 550 (1985); K. Jorge, P. Giorgio, and Z. Francesco, J. Stat. Mech.: Theory Exp. 2012, P10012; M. Robles, M. López de Haro, and A. Santos, J. Chem. Phys. 120, 9113 (2004); 126, 016101 (2007); R. D. Rohrmann, M. Robles, M. López de Haro, and A. Santos, ibid. 129, 014510 (2008); R. D. Rohrmann and A. Santos, Phys. Rev. E 76, 051202 (2007); M. Skoge, A. Donev, F. H. Stillinger, and S. Torquato, ibid. 74, 041127 (2006); S. Torquato, O. U. Uche, and F. H. Stillinger, ibid. 74, 061308 (2006); D. Wyler, N. Rivier, and H. L. Frisch, Phys. Rev. A 36, 2422 (1987). 24J. G. Loeser, Z. Zhen, S. Kais, and D. R. Herschbach, J. Chem. Phys. 95, 4525 (1991); G. Parisi and F. Slanina, Phys. Rev. E 62, 6554 (2000). 25J. K. Percus and G. J. Yevick, Phys. Rev. 110, 1 (1958). 26T. Morita, Prog. Theor. Phys. 20, 920 (1958); 21, 361 (1959). 27T. Morita, Prog. Theor. Phys. 23, 829 (1960). 28J. M. J. van Leeuwen, J. Groeneveld, and J. de Boer, Physica 25, 792 (1959); E. Meeron, J. Math. Phys. 1, 192 (1960); L. Verlet, Nuovo Cimento 18, 77 (1960); M. S. Green, J. Chem. Phys. 33, 1403 (1960); G. S. Rushbrooke, Physica 26, 259 (1960). 29S. Rast, P. H. Fries, and H. Krienke, Mol. Phys. 96, 1543 (1999); S. Labík, H. Gabrielová, J. Kolafa, and A. Malijevský, ibid. 101, 1139 (2003); D. J. Naresh and J. K. Singh, Fluid Phase Equilib. 285, 36 (2009); 279, 47

J. Chem. Phys. 142, 214110 (2015) (2009); T. B. Tan, A. J. Schultz, and D. A. Kofke, Mol. Phys. 109, 123 (2010); K. R. S. Shaul, A. J. Schultz, and D. A. Kofke, J. Chem. Phys. 137, 184101 (2012); J. H. Yang, A. J. Schultz, J. R. Errington, and D. A. Kofke, ibid. 138, 134706 (2013). 30S. K. Kwak and D. A. Kofke, J. Chem. Phys. 122, 104508 (2005). 31K. Hiroike, J. Phys. Soc. Jpn. 12, 326 (1957); L. L. Lee, J. Chem. Phys. 103, 9388 (1995). 32C. Hurst, Proc. Phys. Soc. 86, 193 (1965). 33J. S. Rowlinson, Mol. Phys. 9, 217 (1965). 34J. S. Rowlinson, Mol. Phys. 10, 533 (1966). 35F. Lado, J. Chem. Phys. 47, 4828 (1967); D. Henderson and R. O. Watts, Mol. Phys. 18, 429 (1970); G. Zerah and J. P. Hansen, J. Chem. Phys. 84, 2336 (1986). 36R. E. Caligaris, V. Kuz, and A. E. Rodriguez, J. Chem. Phys. 48, 4794 (1968); Y. Zhou and G. Stell, ibid. 92, 5533 (1990); P. Attard, ibid. 95, 4471 (1991); D. Henderson and S. Sokołowski, ibid. 104, 2971 (1996). 37T. Gaskell, Phys. Lett. A 27, 209 (1968). 38G. Stell, Mol. Phys. 16, 209 (1969); Y. Rosenfeld and N. W. Ashcroft, Phys. Rev. A 20, 1208 (1979). 39P. Hutchinson and W. R. Conkie, Mol. Phys. 21, 881 (1971). 40P. Hutchinson and W. R. Conkie, Mol. Phys. 24, 567 (1972). 41L. Verlet, Mol. Phys. 41, 183 (1980). 42D. S. Hall and W. R. Conkie, Mol. Phys. 40, 907 (1980). 43G. A. Martynov and G. N. Sarkisov, Mol. Phys. 49, 1495 (1983). 44F. J. Rogers and D. A. Young, Phys. Rev. A 30, 999 (1984). 45P. Ballone, G. Pastore, G. Galli, and D. Gazzillo, Mol. Phys. 59, 275 (1986). 46D. M. Duh and A. D. J. Haymet, J. Chem. Phys. 103, 2625 (1995); J. M. Bomont and J. L. Bretonnet, ibid. 114, 4141 (2001); N. Choudhury and S. K. Ghosh, ibid. 116, 8517 (2002). 47M. Marucho and B. Montgomery Pettitt, J. Chem. Phys. 126, 124107 (2007). 48R. J. Baxter, J. Chem. Phys. 47, 4855 (1967). 49G. S. Rushbrooke and H. I. Scoins, Proc. R. Soc. A 216, 203 (1953). 50J. Yvon, Actualites Scientifiques et Industrietles (Hermann & Campagnie, Paris, 1935); M. Born and H. S. Green, Proc. R. Soc. A 188, 10 (1946). 51J. G. Kirkwood, J. Chem. Phys. 3, 300 (1935); G. Stell, ibid. 36, 1817 (1962). 52W. G. Hoover and J. C. Poirier, J. Chem. Phys. 37, 1041 (1962). 53T. Morita and K. Hiroike, Prog. Theor. Phys. 23, 1003 (1960). 54S. J. Singer and D. Chandler, Mol. Phys. 55, 621 (1985). 55T. Morita and K. Hiroike, Prog. Theor. Phys. 25, 537 (1961). 56R. M. Adin, J. Phys. A 25, L433 (1992); A. Baram and M. Luban, ibid. 19, L585 (1986). 57A. Baram and J. S. Rowlinson, Mol. Phys. 74, 707 (1991). 58B. D. McKay, Congr. Numer. 30, 45 (1981); B. D. McKay and A. Piperno, J. Symbolic Comput. 60, 94 (2014). 59J. Perkyns and B. M. Pettitt, Theor. Chem. Acc. 96, 61 (1997); J. S. Perkyns, K. M. Dyer, and B. M. Pettitt, J. Chem. Phys. 116, 9404 (2002). 60Stampede, available at http://www.tacc.utexas.edu/stampede/. 61SuperMIC, available at http://www.hpc.lsu.edu/resources/hpc/system.php? system=SuperMIC. 62See supplementary material at http://dx.doi.org/10.1063/1.4921790 for a description of numerical and theoretical details. We include discussion of Fourier transforms in high dimensions, grid choices, series transformations, and details of some of the less common integral equations.

Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.