Radial Distribution Function for a Hard Sphere Liquid: A Modified Percus-Yevick and Hypernetted-Chain Closure Relations

Abstract

Establishment of the radial distribution function by solving the Ornstein-Zernike equation is still an important problem, even more than a hundred years after the original paper publication. New strategies and approximations are common in the literature. A crucial step in this process consists in defining a closure relation which retrieves correlation functions in agreement with experiments or molecular simulations. In this paper, the functional Taylor expansion, as proposed by J. K. Percus, is applied to introduce two new closure relations: one that modifies the Percus Yevick closure relation and another one modifying the Hypernetted-Chain approximation. These new approximations will be applied to a hard sphere system. An improvement for the radial distribution function is observed in both cases. For some densities a greater accuracy, by a factor of five times compared to the original approximations, was obtained.

Keywords:
Ornstein-Zernike equation; correlation functions; functional Taylor expansion; hard sphere system

Introduction

The study of liquids is important since it involves systems with several applications in biological sciences, condensed matter systems and industry, for example. Predicting thermodynamic properties of these systems will demand some theoretical approach, such as molecular dynamics, Monte Carlo (MC) or the integral equation theory (IET). The IET has the advantage of being less time consuming. However, the accuracy of this approach will depend on some approximations such as the model potential of interaction in the system and the relation between correlation functions, which will be discussed.

One well known IET to obtain distribution functions, which can be used to calculate thermodynamic quantities, was proposed by Ornstein and Zernike.¹ Nevertheless, this approach needs an additional equation to be solved, called closure relation. Several approximations for the closure relation have been published, such as Percus-Yevick (PY),² Hypernetted-Chain (HNC)³,⁴ Verlet,⁵ Tsednee and Luchko⁶ and Carvalho and Braga.⁷

Another method consists in retrieving the distribution functions directly from experimental data, which has been carried out routinely by another methods, as in the works of Kaplow et al.,⁸ North et al.,⁹ Yarnell et al.¹⁰ and more recently the works by Amano et al.¹¹ and Carvalho and Braga.¹²-¹⁴

The PY closure relation is in good agreement with MC calculations for low densities, in which it is observed pressure consistency. However, both PY and HNC approximations diverge from the expected result at higher densities, consequently leading to the thermodynamic inconsistency. New closure relations to avoid this issue have been published recently.⁶,⁷

In this work, the functional Taylor expansion method, proposed by Percus,¹⁵ will be applied to obtain modified PY and HNC closure relations to correct the results at higher densities.

Methodology

Theoretical background

The radial distribution, g(r), and direct correlation, C(r), functions can be acquired by solving the Ornstein-Zernike equation, given by¹

(1)$$\mathrm{h}(\mathbf{r})=\mathrm{C}(\mathbf{r})+\mathrm{\rho }\int \mathrm{C}(\mathbf{r}-\mathbf{s})\mathrm{h}(\mathbf{s})\mathrm{d}\mathbf{s}$$

with ρ being the number density of the system. Since there are two unknown functions an additional equation, called closure relation, relating h(r) to C(r) is necessary. In 1962, J. K. Percus¹⁵ proposed the functional Taylor expansion procedure to obtain such relations. In this method one sets an external field, U_ex, acting on the system with this interaction generated by an extra particle placed at the origin, r₀. It is also necessary to define two functionals, A and B, called generating and independent functionals, respectively. The functional A is expanded with respect to B around U_ex = 0. This expansion is given as¹⁶

(2)$$\mathrm{A}({\mathbf{r}}_1\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})=\mathrm{A}\left({\mathbf{r}}_1\right)+\int [\mathrm{B}({\mathbf{r}}_2\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})-\mathrm{B}\left({\mathbf{r}}_2\right)]{\left(\frac{\mathrm{\delta }\mathrm{A}({\mathbf{r}}_1\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})}{\mathrm{\delta }\mathrm{B}({\mathbf{r}}_2\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})}\right)}_{{\mathrm{U}}_{\mathrm{c}\mathrm{x}}=0}{\mathrm{d}\mathrm{r}}_2+\dots$$

with

(3)$$\frac{\mathrm{\delta }\mathrm{A}({\mathbf{r}}_1\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})}{\mathrm{\delta }\mathrm{B}({\mathbf{r}}_2\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})}=\int \frac{\mathrm{\delta }\mathrm{A}({\mathbf{r}}_1\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})}{\mathrm{\delta }{\mathrm{U}}_{\mathrm{e}\mathrm{x}}\left({\mathbf{r}}_3\right)}\frac{\mathrm{\delta }{\mathrm{U}}_{\mathrm{e}\mathrm{x}}\left({\mathbf{r}}_3\right)}{\mathrm{\delta }\mathrm{B}({\mathbf{r}}_2\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})}$$

and

(4)$$\frac{\mathrm{\delta }{\mathrm{U}}_{\mathrm{e}\mathrm{x}}\left({\mathbf{r}}_3\right)}{\mathrm{\delta }\mathrm{B}({\mathbf{r}}_2\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})}=-\mathrm{\beta }[\frac{\mathrm{\delta }({\mathbf{r}}_2-{\mathbf{r}}_3)}{{\mathrm{\rho }}^{(1)}({\mathbf{r}}_2\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})}-\mathrm{C}({\mathbf{r}}_2,{\mathbf{r}}_3\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})]$$

in which $$\mathrm{\beta }=\frac{1}{{\mathrm{k}}_{\mathrm{B}}\mathrm{T}}$$ (where T is the temperature and k_B is the Boltzmann constant), ρ(¹)(r|U_ex) is the one body density and C(r,r’|U_ex) the direct correlation function in the presence of an external field. The key step of this method is to define properly the functionals A and B. The first choice would be to expand the one-body density with respect to the external potential, which leads to the Yvon integral equation.¹⁷

To retrieve the PY closure relation, J. K. Percus set B(r|U_ex) = ρ(¹)(r|U_ex) and . This choice for the generating functional can be understood based on the following relation¹⁶

(5)$${\mathrm{\rho }}^{(1)}(\mathbf{r}\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})=\frac{{\mathrm{\rho }}^{(2)}({\mathbf{r}}_0,\mathbf{r})}{{\mathrm{\rho }}^{(1)}\left({\mathbf{r}}_0\right)}$$

In this equation, ρ(²)(r₀,r) and ρ(¹)(r₀) are the two- and one-body correlation functions for the N + 1 particle system. Since ρ(²)(r₀,r) = ρ(¹)(r₀)ρ(¹)(r)g(r₀,r),¹⁶ one is left with

(6)$${\mathrm{\rho }}^{(1)}(\mathbf{r}\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})={\mathrm{\rho }}^{(1)}(\mathbf{r})\mathrm{g}({\mathbf{r}}_0,\mathbf{r})$$

Therefore, the one body density for a system of N particles in an external field was written as the product of the radial distribution function for the N +1 particle system and the one body density for the N particle system in the absence of an external field. J. K. Percus defined ρ(¹)(r) as the generating functional considering the low-density limit for the radial distribution function, i.e.,, thus

(7)$$\mathrm{A}(\mathbf{r}\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})={\mathrm{\rho }}^{(1)}(\mathbf{r})=\frac{{\mathrm{\rho }}^{(1)}(\mathbf{r}\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})}{\mathrm{g}({\mathbf{r}}_0,\mathbf{r})}={\mathrm{\rho }}^{(1)}(\mathbf{r}\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}}){\mathrm{e}}^{\mathrm{\beta }\mathrm{U}\left({\mathbf{r}}_0\mathbf{r}\right)}$$

Applied in equation 2, this generating functional yields the Percus-Yevick (PY) closure relation,¹⁶ given by

(8)$$\mathrm{C}(\mathbf{r})=(\mathrm{h}(\mathbf{r})+1)(1-{\mathrm{e}}^{\mathrm{\beta }\mathrm{U}(\mathbf{r})})$$

The HNC closure relation can be obtained by setting the generating functional as

(9)$$\mathrm{A}(\mathbf{r}\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})=\ln \left({\mathrm{\rho }}^{(1)}(\mathbf{r}\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}}){\mathrm{e}}^{\mathrm{\beta }\mathrm{U}({\mathbf{r}}_0,\mathbf{r})}\right)$$

which is given by

(10)$$\mathrm{C}(\mathbf{r})=\mathrm{h}(\mathbf{r})-\ln (\mathrm{h}(\mathbf{r})+1)-\mathrm{\beta }\mathrm{U}(\mathbf{r})$$

As exemplified, each choice for the functional to be expanded along with the Ornstein-Zernike equation will result in a different integral equation to be solved. In this work two variations will be proposed as to modify the PY and HNC closure relations.

Results and Discussion

Theoretical results

A liquid of hard spheres particles, with potential given by

(11)$$\mathrm{U}(\mathbf{r})=\{\begin{array}{l}\mathrm{\infty },\text{if}|\mathrm{r}|<\mathrm{\sigma }\\ 0,\text{if}|\mathrm{r}|\ge \mathrm{\sigma }\end{array}$$

in which σ is the diameter of the sphere, will be discussed. The PY closure relation gives good results for this system at relatively low densities, σ³ρ < 0.4, as observed by Trokhymchuk et al.¹⁸ For higher densities, PY approximation underestimates the radial distribution function first peak intensity. In this sense, one may consider the approximation for g(r) in equation 7 sufficient only for small densities.

It is proposed a better approximation for slightly higher densities, for which PY closure relation does not give good results. One has¹⁹

(12)$$\mathrm{g}({\mathbf{r}}_1,{\mathbf{r}}_2)={\mathrm{e}}^{-\mathrm{\beta }\mathrm{U}({\mathbf{r}}_1,{\mathbf{r}}_2)}[1+\mathrm{\rho }\int ({\mathrm{e}}^{-\mathrm{\beta }\mathrm{U}({\mathbf{r}}_1,{\mathbf{r}}_3)}-1)({\mathrm{e}}^{-\mathrm{\beta }\mathrm{U}({\mathbf{r}}_2,{\mathbf{r}}_3)}-1)\mathrm{d}{\mathbf{r}}_3+\dots ]$$

Following the same reasoning, it can be proposed two generating functionals given by

(13)$$\mathrm{A}(\mathbf{r}\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})={\mathrm{\rho }}^{(1)}(\mathbf{r}\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})\mathrm{f}\left(\mathrm{\beta }\mathrm{U}({\mathbf{r}}_0,\mathbf{r})\right)$$

or

(14)$$\mathrm{A}(\mathbf{r}\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})=\ln \left({\mathrm{\rho }}^{(1)}(\mathbf{r}\mid {\mathrm{U}}_{\mathrm{e}\mathrm{x}})\mathrm{f}\left(\mathrm{\beta }\mathrm{U}({\mathbf{r}}_0,\mathbf{r})\right)\right)$$

in which $$\mathrm{f}\left(\mathrm{\beta }\mathrm{U}({\mathbf{r}}_0,\mathbf{r})\right)=\frac{1}{\mathrm{g}({\mathbf{r}}_0,\mathbf{r})}$$.

The generating functional given in equation 13 was recently applied in another work by the authors,⁷ considering f as the Mittag-Leffler function. In the present study, the closure relation proposed consists in using the explicit terms of equation 12 in equation 13, leading to the result

(15)$$\mathrm{C}(\mathbf{r})=(\mathrm{h}(\mathbf{r})+1)(1-\frac{{\mathrm{e}}^{\mathrm{\beta }\mathrm{U}(\mathbf{r})}}{1+\mathrm{\rho }\int ({\mathrm{e}}^{-\mathrm{\beta }\mathrm{U}({\mathbf{r}}_1,{\mathbf{r}}_3)}-1)({\mathrm{e}}^{-\mathrm{\beta }\mathrm{U}({\mathbf{r}}_2,{\mathbf{r}}_3)}-1)\mathrm{d}{\mathbf{r}}_3})$$

referred hereafter as modified Percus-Yevick (mPY) closure relation. The integral was considered constant with respect to variations in external field.

The second closure relation proposed is obtained if in equation 14 one defines

(16)$$\mathrm{f}\left(\mathrm{\beta }\mathrm{U}({\mathbf{r}}_0,{\mathbf{r}}_1)\right)=\mathrm{\Gamma }\left({\mathrm{\beta }}^{\mathrm{\prime }}\right){\mathrm{E}}_{{\mathrm{\alpha }}^{\mathrm{\prime }},{\mathrm{\beta }}^{\mathrm{\prime }}}\left(\mathrm{\beta }\mathrm{U}({\mathbf{r}}_0,{\mathbf{r}}_1)\right)=\mathrm{\Gamma }\left({\mathrm{\beta }}^{\mathrm{\prime }}\right)\sum _{\mathrm{k}=0}^{\mathrm{\infty }}\frac{{\left(\mathrm{\beta }\mathrm{U}({\mathbf{r}}_0,{\mathbf{r}}_1)\right)}^{\mathrm{k}}}{\mathrm{\Gamma }(\mathrm{k}{\mathrm{\alpha }}^{\mathrm{\prime }}+{\mathrm{\beta }}^{\mathrm{\prime }})}$$

with $${\mathrm{E}}_{{\mathrm{\alpha }}^{\mathrm{\prime }},{\mathrm{\beta }}^{\mathrm{\prime }}}\left(\mathrm{\beta }\mathrm{U}({\mathbf{r}}_0,{\mathbf{r}}_1)\right)$$ the two parameters Mittag-Leffler function,²⁰ which generalizes the exponential function and Γ(x) the gamma function. The multiplying factor guarantees $$\underset{{\mathrm{U}}_{\mathrm{e}\mathrm{x}}\to 0}{lim}\mathrm{f}\left(\mathrm{\beta }\mathrm{U}({\mathbf{r}}_0,{\mathbf{r}}_1)\right)=1$$. The free parameters, α’ and β’, will give flexibilization to the closure relation, referred as Mittag-Leffler Hypernetted-Chain (MLHNC) closure relation, represented as

(17)$$\mathrm{C}(\mathbf{r})=\frac{[\mathrm{h}(\mathbf{r})-\ln (\mathrm{h}(\mathbf{r})+1)-\ln (\mathrm{\Gamma }\left({\mathrm{\beta }}^{\mathrm{\prime }}\right){\mathrm{E}}_{{\mathrm{\alpha }}^{\mathrm{\prime }}{\mathrm{\beta }}^{\mathrm{\prime }}}(\mathrm{\beta }\mathrm{U}(\mathbf{r})))]\mathrm{\Gamma }({\mathrm{\alpha }}^{\mathrm{\prime }}+{\mathrm{\beta }}^{\mathrm{\prime }})}{\mathrm{\Gamma }\left({\mathrm{\beta }}^{\mathrm{\prime }}\right)}$$

In both equations 15 and 17 r = r₀ - r₁. These closure relations will be used to solve Ornstein-Zernike equation for the hard sphere model. Since the structure of the fluid will only depend on the distance from the reference particle, one can represent the correlation functions as functions of r = |r|.

Relation with the bridge function

A general closure relation can be acquired through the diagrammatic expansion method¹⁷ and is given by

(18)$$\mathrm{C}(\mathbf{r})=\mathrm{h}(\mathbf{r})-\ln (\mathrm{g}(\mathbf{r}))-\mathrm{\beta }\mathrm{U}(\mathbf{r})+\mathrm{B}(\mathbf{r})$$

in which B(r) is the Bridge function. For a specific approximation of B(r) it is possible to retrieve closure relations as, for example, the PY or HNC.

A direct connection between equations 15 and 17 and the diagrammatic expansions may be very difficult to achieve. However, from equation 18 it is possible to determine the Bridge functions which leads to the closure relations obtained previously. For equation 15 one has

(19)$$\mathrm{B}(\mathbf{r})=\ln (\mathrm{\tau }(\mathbf{r}))+\ln (\mathrm{h}(\mathbf{r})+1-\mathrm{C}(\mathbf{r}))-\mathrm{h}(\mathbf{r})+\mathrm{C}(\mathbf{r})$$

with $$\mathrm{\tau }(\mathbf{r})=1+\mathrm{\rho }\int ({\mathrm{e}}^{-\mathrm{\beta }\mathrm{U}({\mathbf{r}}_1,{\mathbf{r}}_3)}-1)({\mathrm{e}}^{-\mathrm{\beta }\mathrm{U}({\mathbf{r}}_2,{\mathbf{r}}_3)}-1)\mathrm{d}{\mathbf{r}}_3$$. Equation 17 is acquired by setting

(20)$$\mathrm{B}(\mathbf{r})=\mathrm{\beta }\mathrm{U}(\mathbf{r})-\ln (\mathrm{\Gamma }\left({\mathrm{\beta }}^{\mathrm{\prime }}\right){\mathrm{E}}_{{\mathrm{\alpha }}^{\mathrm{\prime }}{\mathrm{\beta }}^{\mathrm{\prime }}}(\mathrm{\beta }\mathrm{U}(\mathbf{r})))+(1-\frac{\mathrm{\Gamma }({\mathrm{\alpha }}^{\mathrm{\prime }}+{\mathrm{\beta }}^{\mathrm{\prime }})}{\mathrm{\Gamma }\left({\mathrm{\beta }}^{\mathrm{\prime }}\right)})\mathrm{C}(\mathbf{r})$$

In the first case, if ρ → 0 it is obtained the usual Bridge function for PY closure relation. For the last case, if α’ = β’ = 1 one will be left with B(r) = 0, which is the Bridge function which retrieves the HNC approximation. Therefore, a connection between the closure relations proposed in this work and the Bridge function was possible.

Numerical details to solve Ornstein-Zernike equation

The numerical solution of Ornstein-Zernike equation can be simplified if the indirect correlation function, defined as γ(r) = h(r) - C(r), is used. As one can observe, in equation 15 the closure relation involves the term e^βU(r) on numerator, which can cause numerical instability while computing solutions.

In this sense equations 15 and 17 can be modified to

(21)$$\mathrm{C}(\mathrm{r})=(\mathrm{\gamma }(\mathrm{r})+1)[\frac{1+\mathrm{\rho }\mathrm{I}(\mathrm{r})}{{\mathrm{e}}^{\mathrm{\beta }\mathrm{U}(\mathrm{r})}}-1]$$

which has the exponential on denominator, and

(22)$$\mathrm{C}(\mathrm{r})={\mathrm{e}}^{-\ln (\mathrm{\Gamma }\left({\mathrm{\beta }}^{\mathrm{\prime }}\right){\mathrm{E}}_{\mathrm{\alpha },\mathrm{\beta }}(\mathrm{\beta }\mathrm{B}(\mathrm{r}))+\mathrm{\gamma }(\mathrm{r})+\left(1\frac{\mathrm{\Gamma }\left({\mathrm{\beta }}^{\mathrm{\prime }}\right)}{\mathrm{\Gamma }({\mathrm{\alpha }}^{\mathrm{\prime }}+{\mathrm{\beta }}^{\mathrm{\prime }})}\right)\mathrm{c}(\mathrm{r})}-\mathrm{\gamma }(\mathrm{r})-1$$

with I(r) representing the integral in equation 15 and can be calculated using bipolar coordinates.²¹ The usual PY and HNC closure relations are retrieved for ρ → 0 in equation 18 and α’ = β’ = 1 in equation 22. The closure relation 22 is a transcendental equation and must be solved at each iteration, demanding more computational time.

Applying the Fourier transform in equation 1 and using $$\hat{\mathrm{\gamma }}(\mathrm{k})=\hat{\mathrm{h}}(\mathrm{k})-\hat{\mathrm{C}}(\mathrm{k})$$, for $$\hat{\mathrm{f}}(\mathrm{k})$$ representing the Fourier transform of function f(r), one has

(23)$$\hat{\mathrm{\gamma }}(\mathrm{k})=\frac{\mathrm{\rho }\hat{\mathrm{C}}(\mathrm{k}{)}^2}{1-\mathrm{\rho }\hat{\mathrm{C}}(\mathrm{k})}$$

An initial guess for γ(r) is made and equation 21 or 22 is used to obtain the direct correlation function. The Fourier transform of this function is calculated and the result is used in equation 23 for acquisition of $$\hat{\mathrm{\gamma }}(\mathrm{k})$$. Its inverse Fourier transform is carried out and a new indirect correlation function is retrieved. If this algorithm become unstable for some density or temperature one has to mix the solutions as γ_i+1(r) = ωγ_i(r) + (1 - ω)γ_i+1(r), with 0 ≤ ω ≤ 1.²²

The calculations were carried out in reduced units, i.e., $${\mathrm{r}}^{\ast }=\frac{\mathrm{r}}{\mathrm{\sigma }}$$ and ρ* = σ³ρ. The Fourier transform and its inverse were calculated using r* and k* from 0 to 15 with step of 0.01.

The integral I(r),

(24)$$\mathrm{I}(\mathrm{r})=\frac{2\mathrm{\pi }}{\mathrm{r}}{\int }_0^{\mathrm{\infty }}\mathrm{s}({\mathrm{e}}^{-\mathrm{\beta }\mathrm{U}(\mathrm{s})}-1){\int }_{|\mathrm{r}-\mathrm{s}|}^{\mathrm{r}+\mathrm{s}}\mathrm{t}({\mathrm{e}}^{-\mathrm{\beta }\mathrm{U}(\mathrm{t})}-1)\mathrm{d}\mathrm{t}\mathrm{d}\mathrm{s}$$

was solved in the interval 0 ≤ r ≤ 3 using bipolar coordinates, since this function vanishes for larger values of r. It was set 0 ≤ s ≤ 1.2 and ds = dt = 5 × 10^-4. The Mittag-Leffler function was calculated using an algorithm by Podlubny.²³

Results for g(r) using mPY closure relation

The mPY closure relation proposed has a density dependence since the second term in expansion (equation 12) is considered, unlike the J. K. Percus approximation. Since the PY closure relation, which is derived under the low-density limit approximation, underestimates the first peak intensity, adding a second term to the generating functional will gives a closure relation which corrects the radial distribution function for higher densities, as can be seen in Figure 1. The results are plotted along those acquire by Yu and Wu,²⁴ which has a great accuracy if compared with Monte Carlo (MC) calculations.

The absolute difference between the predicted g(r) and that found in literature is presented in Figure 2. Since the behavior of this difference is similar in all cases, only data for ρ* = 0.8 is given. It is possible to see an improvement with respect to the first peak intensity while for larger values of r the absolute error oscillates between higher and smaller values, although these differences are about 0.1. However, the norm of the differences is smaller for the new closure relation proposed at all densities, as can be seen in Table 1.

It was observed an improvement of the radial distribution function, but at cost of the direct correlation function for this closure relation. This can be explained since PY closure relation already gives accurate results for C(r).

Figure 1
Radial distribution function results for hard spheres obtained by the Percus-Yevick approximation (dashed line), the modified Percus-Yevick (full line) and results reported by Yu and Wu²⁴ (×) at ρ* = 0.7, 0.8 and 0.9. For clarity the curves are shifted in the y-axis by a factor of 0, 0.5 and 1, respectively.

Figure 2
Difference between results reported by Yu and Wu²⁴ and those acquired by the Percus-Yevick closure relation (o) and the one proposed in this work (×) for ρ* = 0.8.

Results for g(r) using MLHNC closure relation

The second closure relation proposed using the two parameters Mittag-Leffler function, equation 16, does not contain any density dependence. However, the parameters give flexibilization to the closure relation, since for each set of parameters it will be necessary to solve a different transcendental equation, represented by equation 19. This feature enables one to obtain results in good agreement with those reported by Yu and Wu,²⁴ as represented in Figure 3 for α’ = 1 and β’ = 0.7.

The MLHNC closure relation represents an overall improvement if compared to the usual HNC closure relation, which superestimates the first peak intensity of g(r) for the hard sphere model. In Figure 4 it is represented the difference between the results obtained by HNC and MLHNC closure relations at ρ* = 0.8 with those reported in literature.²⁴

Table 1
Error norms obtained using the Percus-Yevick (PY), Hypernetted-Chain (HNC), Mittag-Leffler-Hypernetted-Chain (MLHC) and the modified Percus-Yevick (mPY) closure relations

r*	||g(r) - g(r)PY||	||g(r) - g(r)mPY||	||g(r) - g(r)HNC||	||g(r) - g(r)MLHNC||
0.7	0.5043	0.3853	0.7367	0.2403
0.8	0.7342	0.2564	0.8884	0.1663
0.9	1.3781	0.4049	1.5778	0.2708

³

Figure 3
Same as in Figure 1, changing PY and mPY by Hypernetted-Chain (dashed line) and Mittag-Leffler-Hypernetted-Chain (full line) closure relations.

Figure 4
Same as in Figure 2, changing PY and mPY by Hypernetted-Chain (o) and Mittag-Leffler-Hypernetted-Chain (×) closure relations.

As already mentioned, only data for ρ* = 0.8 is given due the simillarity with other densities. From Figure 4 it is observed a similar absolute difference variation for both closure relations at higher values of r, although MLHNC closure relation give the smaller values.

In this approximation C(r) is acquired with a smaller error if compared to the HNC closure relation, but it is not accurate if compared with those obtained by MC calculations.²⁵

The results presented for the radial distribution error norm given in Table 1 demonstrates an overall improvement for both closure relations proposed in this work. The smaller values are observed for the MLHNC approximation, followed in crescent order by mPY, PY and HNC.

Thermodynamic properties

In this section some thermodynamic properties, such as the compressibility factor, isothermal compressibility and chemical potential, will be calculated using the proposed closure relations and compared with those acquired by the PY and HNC approximations.

For a hard sphere system the compressibility can be calculated from virial pressure, P_v, as²⁶

(25)$$\frac{\mathrm{\beta }{\mathrm{P}}_{\mathrm{v}}}{\mathrm{\rho }}=1+4\mathrm{\eta }\mathrm{g}\left({\mathrm{\sigma }}^{+}\right)$$

in which $$\mathrm{\eta }=\frac{\mathrm{\pi }{\mathrm{\rho }}^{\ast }}{6}$$. As for the isothermal compressibility, K_T, one has

(26)$${\mathrm{K}}_{\mathrm{T}}=1-4\mathrm{\pi }\mathrm{\rho }\int {\mathrm{r}}^2\mathrm{C}(\mathrm{r})\mathrm{d}\mathrm{r}$$

Since the thermodynamic consistency was not imposed in this work, the compressibility pressure, P_C, can also be calculated from the isothermal compressibility as

(27)$$\frac{\mathrm{\beta }{\mathrm{P}}_{\mathrm{C}}}{\mathrm{\rho }}=\frac{1}{\mathrm{\rho }}{\int }_0^{\mathrm{\rho }}{\mathrm{K}}_{\mathrm{T}}\mathrm{d}{\mathrm{\rho }}^{\mathrm{\prime }}=\frac{1}{\mathrm{\rho }}{\int }_0^{\mathrm{\rho }}[1-4\mathrm{\pi }{\mathrm{\rho }}^{\mathrm{\prime }}\int {\mathrm{r}}^2\mathrm{C}(\mathrm{r})\mathrm{d}\mathrm{r}]\mathrm{d}{\mathrm{\rho }}^{\mathrm{\prime }}$$

The accuracy can be evaluated by comparing with results obtained from the Carnahan-Starling (CS) equation.²⁷ In this case one has

(28)$$\frac{\mathrm{\beta }{\mathrm{P}}_{\mathrm{C}\mathrm{S}}}{\mathrm{\rho }}=\frac{1+\mathrm{\eta }+{\mathrm{\eta }}^2-{\mathrm{\eta }}^3}{(1-\mathrm{\eta }{)}^3}$$

and

(29)$${\mathrm{K}}_{\mathrm{T}}=\frac{8\mathrm{\eta }-2{\mathrm{\eta }}^2}{(1-\mathrm{\eta }{)}^4}+1$$

Table 2
Results for compressibility acquired by the different approximations and predicted by CS equation

ρ*	$$\frac{\mathrm{\beta }{\mathrm{P}}_{\mathrm{v},\mathrm{P}\mathrm{Y}}}{\mathrm{\rho }}$$	$$\frac{\mathrm{\beta }{\mathrm{P}}_{\mathrm{v},\mathrm{m}\mathrm{P}\mathrm{Y}}}{\mathrm{\rho }}$$	$$\frac{\mathrm{\beta }{\mathrm{P}}_{\mathrm{v},\mathrm{H}\mathrm{N}\mathrm{C}}}{\mathrm{\rho }}$$	$$\frac{\mathrm{\beta }{\mathrm{P}}_{\mathrm{v},\mathrm{M}\mathrm{L}\mathrm{H}\mathrm{N}\mathrm{C}}}{\mathrm{\rho }}$$	$$\frac{\mathrm{\beta }{\mathrm{P}}_{\mathrm{C}\mathrm{S}}}{\mathrm{\rho }}$$
0.7	5.1751	5.8870	6.3692	5.8988	5.7102
0.8	6.7459	7.7668	8.8868	7.9541	7.7497
0.9	8.8942	10.3554	12.6036	10.8502	10.7461

³β=1kBT_Bv

Table 3
Results for isothermal compressibility for each approximation and predicted by CS equation

ρ*	KT,PY	KT,mPY	KT,HNC	KT,MLHNC	KT,CS
0.7	17.7880	8.2688	13.0871	14.7158	17.5392
0.8	27.9770	13.4057	18.6924	21.8747	27.3070
0.9	45.0320	22.2209	27.0617	33.2809	43.5457

³_T

The results calculated by the PY, mPY, HNC and MLHNC closure relations for virial pressure and isothermal compressibility are given in Tables 2 and 3.

As can be inferred from Figures 1 and 3, the virial pressure calculated from PY approximation understemimates the results calculated from CS equation, while those calculated using the HNC closure relation superestimates these values. It is observed good accuracy for the results acquired by both closure relations proposed in this work.

One can conclude that the mPY and MLHNC closure relations do not satisfy the pressure consistency condition from Table 3. Since the values obtained using mPY are lower than those acquired by CS equation, the compressibility pressure will not coincide with P_v. Thus, integration of the isothermal compressibility for calculating P_C is not necessary. For MLHNC closure relation it is possible to infer that the results for P_C will be better than those acquired by HNC, although not pressure consistent.

Therefore, the MLHNC closure relation represents an overall improvement with respect to the pressure and isothermal compressibility calculations. For mPY and PY it is observed a change in accuracy depending on the property of interest, i.e., mPY will be more accurate for properties calculated through the g(r) and PY will lead to more accurate results which depends on C(r).

From equation 19 the chemical potential can be easily calculated, since for hard spheres one has ln y(0) = βµ, in which y(r) is called cavity distribution function. By equation 18 it is possible to write

(30)$$\ln \mathrm{y}(\mathrm{r})=\mathrm{\gamma }(\mathrm{r})+\mathrm{B}(\mathrm{r})$$

The results for PY and mPY approximatons are given in Figure 5. It should be noted for MLHNC closure relation that the Bridge function given in equation 20 is numerically divergent for the range of r < σ. Thus, the procedure adopted for mPY can not be carried out here.

Figure 5
Results for ln y(r) in the PY (dashed lines) and mPY (full lines). From bottom to top, the densities are 0.7, 0.8 and 0.9.

It is observed that both PY and mPY underestimates the results for ln(y(r))²⁶ and, consequently, the chemical potential predicted by the CS equation, calculated as

(31)$$\mathrm{\beta }{\mathrm{\mu }}_{\mathrm{C}\mathrm{S}}=\frac{8\mathrm{\eta }-9{\mathrm{\eta }}^2+3{\mathrm{\eta }}^3}{(1-\mathrm{\eta }{)}^3}$$

for which the results at ρ* = 0.7, 0.8 and 0.9 are, respectively, βµ_CS = 7.36, 10.15 and 14.16.

Conclusions

Two new closure relations which modifies the usual PY and HNC approximations were obtained by means of functional Taylor expansion method. The modified generating functional, which takes into account the second term in the expansion for g(r) to obtain the mPY approximation, is based on physical grounds, correcting the radial distribution function at higher densities. The change proposed to the HNC closure relation enables one to modify the exponential function correcting the closure relation behavior to guarantee better results for g(r).

Therefore, the method proposed by J. K. Percus can be used to derive new closure relations as to increase the results accuracy for important model systems. Since the HNC closure relation is appropriate for charged systems, the MLHNC proposed in this work could also be used to refine results which deviates from simulated or experimental data.

The Ornstein-Zernike equation solutions for the hard sphere system demonstrated an improvement of the results for the radial distribution function considering both new closure relations proposed at higher densities. Since closure relations that describes correlation functions in liquids are very important to understand such systems, the new approximations given in this work are important to complement the already existing set of closure relations and contribute to demonstrate the usefulness of J. K. Percus method.

Acknowledgments

We would like to thank CNPq for the financial support.

References

• ¹
  Ornstein, L. S.; Zernike, F.; Proc.-Acad. Sci. Amsterdam 1914, 17, 793.
• ²
  Percus, J. K.; Yevick, G. J.; Phys. Rev. 1958, 110, 1.
• ³
  van Leeuwen, J. M. J.; Groeneveld, J.; de Boer, J.; Physica 1959, 25, 792.
• ⁴
  Morita, T.; Prog. Theor. Phys. 1958, 20, 920.
• ⁵
  Verlet, L.; Mol. Phys. 1980, 41, 183.
• ⁶
  Tsednee, T.; Luchko, T.; Phys. Rev. E 2019, 99, 032130.
• ⁷
  Carvalho, F. S.; Braga, J. P.; Phys. A 2021, 576, 126065.
• ⁸
  Kaplow, R.; Strong, S. L.; Averbach, B. L.; Phys. Rev. 1965, 138, A1336.
• ⁹
  North, D. M.; Enderby, J. E.; Egelstaff, P. A.; J. Phys. C: Solid State Phys. 1968, 1, 1075.
• ¹⁰
  Yarnell, J. L.; Katz, M. J.; Wenzel, R. G.; Koenig, S. H.; Phys. Rev. A 1973, 7, 2130.
• ¹¹
  Amano, K.-I.; Sawazumi, R.; Imamura, H.; Sumi, T.; Hashimoto, K.; Fukami, K.; Kitaoka, H.; Nishi, N.; Sakka, T.; Chem. Lett. 2020, 49, 1017.
• ¹²
  Carvalho, F. S.; Braga, J. P.; Alves, M. O.; Gonçalves, C. E. M.; Theor. Chem. Acc. 2020, 139, 29.
• ¹³
  Carvalho, F. S.; Braga, J. P.; J. Mol. Model. 2020, 26, 193.
• ¹⁴
  Carvalho, F. S.; Braga, J. P.; Braz. J. Phys. 2020, 50, 489.
• ¹⁵
  Percus, J. K.; Phys. Rev. Lett. 1962, 8, 462.
• ¹⁶
  Balescu, R.; Equilibrium and Nonequilibrium Statistical Mechanics; Wiley: New York, USA, 1975.
• ¹⁷
  Hansen, J.-P.; McDonald, I. R.; Theory of Simple Liquids: with Applications to Soft Matter; Academic Press: Oxford, UK, 2013.
• ¹⁸
  Trokhymchuk, A.; Nezbeda, I.; Jirsák, J.; Henderson, D.; J. Chem. Phys. 2005, 123, 024501.
• ¹⁹
  Carley, D. D.; Lado, F.; Phys. Rev. 1965, 137, A42.
• ²⁰
  Podlubny, I.; Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of their Solution and Some of their Applications; Elsevier: New York, USA, 1998.
• ²¹
  Lee, L. L.; Molecular Thermodynamics of Nonideal Fluids; Butterworth-Heinemann: London, UK, 1988.
• ²²
  Broyles, A. A.; J. Chem. Phys. 1960, 33, 456.
• ²³
  Podlubny, I.; Mittag-Leffler Function, MATLAB Central File Exchange, 2018, available at https://www.mathworks.com/matlabcentral/fileexchange/8738-mittagleffler-function accessed in August 2021.
  » https://www.mathworks.com/matlabcentral/fileexchange/8738-mittagleffler-function
• ²⁴
  Yu, Y. X.; Wu, J.; J. Chem. Phys. 2002, 117, 10156.
• ²⁵
  Groot, R. D.; Van der Eerden, J. P.; Faber, N. M.; J. Chem. Phys. 1987, 87, 2263.
• ²⁶
  Lee, L. L.; J. Chem. Phys. 1995, 103, 9388.
• ²⁷
  Carnahan, N. F.; Starling, K. E.; J. Chem. Phys 1969, 51, 635.

Publication Dates

• Publication in this collection
  26 Nov 2021
• Date of issue
  Dec 2021

History

• Received
  15 July 2021
• Accepted
  18 Aug 2021