The standard mean-field treatment of inter-particle attraction in classical DFT is better than one might expect

TL;DR

This paper demonstrates that the standard mean-field approach in classical DFT for attractive interactions provides more accurate inhomogeneous fluid predictions than expected, surpassing the simple RPA approximation in a one-dimensional model.

Contribution

The study reveals that mean-field DFT, when used with the test-particle method, outperforms RPA in predicting fluid structure, challenging assumptions about its accuracy based on bulk correlation functions.

Findings

Mean-field DFT yields superior inhomogeneous fluid structure results.

Results from a 1D model show better accuracy than RPA predictions.

The quality of DFT should not be judged solely by bulk correlation approximations.

Abstract

In classical density functional theory (DFT) the part of the Helmholtz free energy functional arising from attractive inter-particle interactions is often treated in a mean-field or van der Waals approximation. On the face of it, this is a somewhat crude treatment as the resulting functional generates the simple random phase approximation (RPA) for the bulk fluid pair direct correlation function. We explain why using standard mean-field DFT to describe inhomogeneous fluid structure and thermodynamics is more accurate than one might expect based on this observation. By considering the pair correlation function $g(x)$ and structure factor $S(k)$ of a one-dimensional model fluid, for which exact results are available, we show that the mean-field DFT, employed within the test-particle procedure, yields results much superior to those from the RPA closure of the bulk Ornstein-Zernike equation. We argue that one should not judge the quality of a DFT based solely on the approximation it generates for the bulk pair direct correlation function.

Figures

Equations

$$h(r)=c^{(2)}(r)+{\rho }_b\int \text{d}{\mathbf{r}}^{\prime }{c}^{(2)}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)h(r^{\prime })$$

$$\hat{h}(k)=\frac{\hat{c}(k)}{1-{\rho }_b\hat{c}(k)},$$

$$S(k)=1+{\rho }_b\hat{h}(k).$$

$$c^{(2)}(r)=h(r)-\ln (h(r)+1)-\beta v(r)+B(r)$$

$$v(r)=v_0(r)+v_1(r),$$

$$c_{\mathrm{RPA}}^{(2)}(r)\equiv c_0^{(2)}(r)-\beta v_1(r),$$

$$S_{\mathrm{RPA}}(k)=\frac{S_0(k)}{1+{\rho }_b\beta {\hat{v}}_1(k)S_0(k)}$$

$$F[\rho (\mathbf{r})]=$$

$$v_{\lambda }(r)=v_0(r)+\lambda v_1(r),0\le \lambda \le 1.$$

$${\rho }_{\lambda }^{(2)}(\mathbf{r},{\mathbf{r}}^{\prime })\approx \rho (\mathbf{r})\rho ({\mathbf{r}}^{\prime }),0\le \lambda \le 1.$$

$$F^{\mathrm{ex}}[\rho (\mathbf{r})]\approx F_0^{\mathrm{ex}}[\rho (\mathbf{r})]+\frac{1}{2}\int \text{d}\mathbf{r}\int \text{d}{\mathbf{r}}^{\prime }\rho (\mathbf{r})\rho ({\mathbf{r}}^{\prime })v_1(|\mathbf{r}-{\mathbf{r}}^{\prime }|).$$

$$c_{0}^{(2)}(|{\bf r}-{\bf r}^{\prime}|)=-\frac{\delta^{2}\beta F_{0}^{\mathrm{ex}}[\rho]}{\delta\rho({\bf r})\delta\rho({\bf r}^{\prime})}\bigg{|}_{\rho({\bf r})=\rho_{b}}.$$

$$\frac{\delta \Omega [\rho ]}{\delta \rho }=k_{B}T\ln [{\Lambda }^3\rho (r)]+\frac{\delta F_0^{\mathrm{ex}}[\rho ]}{\delta \rho }+\int \text{d}{\mathbf{r}}^{\prime }\rho (r^{\prime })v_1(|\mathbf{r}-{\mathbf{r}}^{\prime }|)+v(r)-\mu =0,$$

$$\mu=k_{B}T\ln[\Lambda^{3}\rho_{b}]+\frac{\delta F_{0}^{\mathrm{ex}}[\rho]}{\delta\rho}\bigg{|}_{\rho_{b}}+\rho_{b}\int{\textrm{d}}{\bf r}v_{1}(r).$$

$$\displaystyle\frac{\delta F_{0}^{\mathrm{ex}}[\rho]}{\delta\rho}=\frac{\delta F_{0}^{\mathrm{ex}}[\rho]}{\delta\rho}\bigg{|}_{\rho_{b}}+\int{\textrm{d}}{\bf r}^{\prime}(\rho({\bf r}^{\prime})-\rho_{b})\frac{\delta^{2}F_{0}^{\mathrm{ex}}[\rho]}{\delta\rho({\bf r})\delta\rho({\bf r}^{\prime})}\bigg{|}_{\rho_{b}}+H_{0}[\rho({\bf r})],$$

$$0=k_{B}T\ln \left(\frac{\rho (r)}{{\rho }_b}\right)+\int \text{d}{\mathbf{r}}^{\prime }(\rho (r^{\prime })-{\rho }_b)\left[-k_{B}T{c}_0^{(2)}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)+v_1(|\mathbf{r}-{\mathbf{r}}^{\prime }|)\right]+H_0[\rho (r)]+v(r),$$

$$\frac{(\rho (r)-{\rho }_b)}{{\rho }_b}=\frac{(\rho (r)-{\rho }_b)}{{\rho }_b}-\ln \left(\frac{\rho (r)}{{\rho }_b}\right)-\beta v(r)-\beta H_0[\rho (r)]+{\rho }_b\int \text{d}{\mathbf{r}}^{\prime }\frac{(\rho (r^{\prime })-{\rho }_b)}{{\rho }_b}\left[c_0^{(2)}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)-\beta v_1(|\mathbf{r}-{\mathbf{r}}^{\prime }|)\right].$$

$$h_0(r)=h_0(r)-\ln \left(h_0(r)+1\right)-\beta v_0(r)-\beta H_0[{\rho }_b{g}_0(r)]+{\rho }_b\int \text{d}{\mathbf{r}}^{\prime }{h}_0(r^{\prime })c_0^{(2)}(|\mathbf{r}-{\mathbf{r}}^{\prime }|).$$

$$S(k)=\frac{1-e^{-\beta [\mu (p+ik/\beta )-2\mu (p)+\mu (p-ik/\beta )]}}{(1-e^{-\beta [\mu (p+ik/\beta )-\mu (p)]})(1-e^{-\beta [\mu (p-ik/\beta )-\mu (p)]})}$$

$$v_{0}(x)=\left\{\begin{array}[]{l l}\infty&\quad{|x|\leq\sigma}\\ 0&\quad{|x|>\sigma}\end{array}\right.$$

$$\beta v_{1}(x)=\left\{\begin{array}[]{l l}0&\quad{|x|\leq\sigma}\\ -z_{p}(\sigma+\sigma_{p}-|x|)&\quad{\sigma<|x|\leq\sigma+\sigma_{p}}\\ 0&\quad{|x|>\sigma+\sigma_{p}}.\end{array}\right.$$

$$F_0^{\mathrm{ex}}[\rho ]\equiv F_{\mathrm{HR}}^{\mathrm{ex}}[\rho ]=-\int \rho (x)\ln \left[1-t(x)\right]\mathrm{d}x,$$

$$t(x)={\int }_{x-\sigma }^x\rho (x^{\prime })\mathrm{d}{x}^{\prime }.$$

$$\frac{\delta F_0[\rho ]}{\delta \rho }=\ln \frac{\Lambda \rho (x)}{1-t(x)}+{\int }_x^{x+\sigma }\frac{\rho (x^{\prime })}{1-t(x^{\prime })}\mathrm{d}{x}^{\prime }.$$

$$c_{\mathrm{HR}}^{(2)}(|x-x^{\prime }|)=-\Theta (\sigma -|x-x^{\prime }|)\frac{1-{\rho }_b|x-x^{\prime }|}{(1-{\rho }_b\sigma {)}^2},$$

$$F^{\mathrm{ex}}[\rho (\mathbf{r})]\approx F^{\mathrm{ex}}[{\rho }_b]+{\mu }_{\mathrm{ex}}\int \text{d}\mathbf{r}(\rho (\mathbf{r})-{\rho }_b)$$

$$-\frac{1}{2\beta }\int \text{d}\mathbf{r}\int \text{d}{\mathbf{r}}^{\prime }(\rho (\mathbf{r})-{\rho }_b)(\rho ({\mathbf{r}}^{\prime })-{\rho }_b)c(|\mathbf{r}-{\mathbf{r}}^{\prime }|).$$

Full text

I Introduction

Classical density functional theory (DFT)Evans (1979, 1992); Löwen (2002); Lutsko (2010); Hansen and McDonald (2013) is based on the idea that the thermodynamic grand potential of an inhomogeneous fluid can be expressed as a functional of the average one-body density profile $\rho({\bf r})$. Minimizing an (approximate) functional with respect to $\rho({\bf r})$ provides an estimate of the equilibrium density profile and the grand potential. DFT has proved to be a versatile tool for determining the thermodynamic quantities arising in the physics of adsorption and fluid interfaces. These include, for example, interfacial tensionsEvans (1992) and the solvation force (excess pressure) for confined fluids.Evans and Marconi (1987); Stewart and Evans (2014) Since DFT provides directly the grand potential it is particularly well-suited to investigations of surface phase behaviour and perhaps it is here where DFT has had most success, revealing new phenomena and new physical insight.Evans (1992); Hansen and McDonald (2013) DFT also provides a direct measure of density fluctuations in the inhomogeneous fluid via the local compressibility $(\partial\rho({\bf r})/\partial\mu)$, where $\mu$ is the chemical potential and recent papers have investigated this quantity for liquids at solvophobic planar substratesEvans and Stewart (2015) and confined between a variety of nanostructures.Chacko et al. (2017) In addition to the one-body density, higher order correlation functions can be obtained by taking further functional derivatives of the free energy functional. In particular, the two-body pair direct correlation function $c^{(2)}({\bf r},{\bf r}^{\prime})$ is proportional to the second functional derivative of the excess Helmholtz free energy with respect to the density.Evans (1979, 1992); Lutsko (2010); Hansen and McDonald (2013) It is tempting to assess the quality of an approximate DFT for a given model fluid by comparing the resulting $c^{(2)}(r)$ for a homogeneous fluid with that obtained from computer simulation of the radial distribution function $g(r)$ or structure factor $S(k)$ or from integral equation theories; see Refs. Evans, 1992; Lutsko, 2010; Hansen and McDonald, 2013 and references therein. In this paper we explain why this might not be the best means of testing the accuracy of a free energy functional.

We focus on the much-used excess Helmholtz free energy functional [Eq. (11) below] that treats the attractive inter-particle interactions in a mean-field (MF) or van der Waals approximation. This standard MF DFT is considered the ‘work-horse’ for applications of DFT to interfaces, adsorption and many other inhomogeneous situations, see e.g. Refs. Evans, 1992; Löwen, 2002; Lutsko, 2010; Hansen and McDonald, 2013; Evans and Marconi, 1987; Stewart and Evans, 2014; Evans and Stewart, 2015; Chacko et al., 2017. Taking two derivatives of this functional and evaluating for a uniform density $\rho_{b}$, corresponding to the (bulk) fluid, leads to the pair direct correlation function $c^{(2)}_{\mathrm{RPA}}(r)$ that forms the basis for the well-known random phase approximation (RPA).Evans (1979); Hansen and McDonald (2013); Likos (2001); Barker and Henderson (1976) Recognizing this connection between the MF DFT and the RPA, one might assume that the fluid structure, thermodynamics and phase behaviour predicted by the MF DFT is of similar quality to that resulting from the RPA closure to the bulk Ornstein-Zernike (OZ) equation. Here, we argue: (i) Results from the MF DFT are more accurate than one might expect from examining bulk pair correlation functions extracted from the RPA closure to the OZ equation and (ii) one should not judge the quality of the MF DFT, or any other approximate DFT, solely upon what the approximate functional generates by taking two functional derivatives. We make these arguments explicit by considering the test-particle limit of Percus,Percus (1962) i.e. we calculate the inhomogeneous fluid density profile $\rho(r)$ around a fixed particle of the same type as the surrounding fluid. This enables us to calculate the fluid radial distribution function $g(r)=\rho(r)/\rho_{b}$. By rewriting the Euler-Lagrange equation obtained from minimizing the MF DFT in this test particle limit and comparing with the OZ equation, we identify additional terms somewhat akin to a hybrid closure relation and certainly different from the RPA. The additional terms also suggest that the MF DFT, treated in the test particle limit, is much superior to standard OZ with the RPA closure. We argue that this observation should carry through more generally for arbitrary external potentials, not just in the test-particle limit.

Our paper proceeds as follows: In Sec. II we describe briefly the RPA in the context of the OZ equation. In Sec. III we describe the MF DFT that generates the RPA and indicate why one might expect this DFT to yield much better results for $g(r)$, within the test-particle procedure, than in the standard RPA treatment. In Sec. IV we illustrate and confirm our arguments by presenting results for $g(x)$ and $S(k)$ for a model one-dimensional (1D) fluid where the pair correlation functions are known exactly and where the MF DFT can be implemented with the exact reference free energy functional, i.e. that pertaining to hard-rods. Finally, in Sec. V we discuss the general implications of our results. We also explain that for the same choice of reference system (hard rods) the MF DFT results are identical to those one would obtain from the local molecular field theory (LMF) of Weeks and co-workersWeeks et al. (1997, 1998); Weeks (2002) for this particular model fluid.

II Integral equation approach: OZ equation and RPA

The two-body pair direct correlation function, $c^{(2)}(r)$, is usually defined via the OZ equation,Hansen and McDonald (2013) which for a uniform and isotropic fluid is

[TABLE]

where $h(r)$ is the total correlation function and $\rho_{b}$ is the density of the (bulk) fluid. The OZ approach to calculating $h(r)$ is to split the correlations present in $h(r)$ into a direct part, which describes the ‘direct’ correlations acting over a range of order that of the interaction pair potential, and an ‘indirect’ part, i.e. the remainder described by the convolution integral. Note that the total correlation function $h(r)=g(r)-1$, where $g(r)$ is the radial distribution function. Fourier transforming Eq. (1) yields an algebraic relation:

[TABLE]

where $\hat{h}(k)$ and $\hat{c}(k)$ are the Fourier transforms of $h(r)$ and $c^{(2)}(r)$, respectively. $\hat{h}(k)$ is closely related to the static structure factorHansen and McDonald (2013)

[TABLE]

Consider a fluid composed of particles interacting via the pair potential $v(r)$. To calculate $h(r)$ one must supplement the OZ Eq. (1) with a further, closure relation between $c^{(2)}(r)$ and $h(r)$. The exact closure relation is usually expressed asHansen and McDonald (2013)

[TABLE]

where $\beta=(k_{B}T)^{-1}$, $k_{B}$ is Boltzmann’s constant, $T$ is the temperature and $B(r)$ is termed the bridge function. $B(r)$ is not known exactly for any 3D fluid. In order to make progress approximations must be made. For example, the hypernetted-chain approximation (HNC) corresponds to $B(r)=0$.Hansen and McDonald (2013); Barker and Henderson (1976) Formally $B(r)$ is the sum of the bridge or elemental diagrams missing in HNC. Suppose $v(r)$ can be split as follows:

[TABLE]

where $v_{0}(r)$ is a suitably chosen reference potential, usually the purely repulsive part of $v(r)$. Then the remainder $v_{1}(r)$ usually incorporates the attractive part of the interaction between particles. The simple closure relation

[TABLE]

where $c_{0}^{(2)}(r)$ is the pair direct correlation function for the (purely repulsive) reference system with the same density $\rho_{b}$, defines the RPA. Note that Eq. (6) enforces the correct asymptotic behaviour: $c^{(2)}(r)\sim-\beta v(r)$, $r\to\infty$ for a fluid away from its critical point. Inserting the Fourier transform of (6) into (2) yields the standard RPA result for the structure factor:Hansen and McDonald (2013)

[TABLE]

where $S_{0}(k)$ is the structure factor of the reference system. The Fourier transform $\hat{v}_{1}(k)$ is assumed to exist. Often the further approximation $c_{0}^{(2)}(r)\approx c_{\mathrm{HS}}^{(2)}(r)$ is made, where $c_{\mathrm{HS}}^{(2)}(r)$ is the pair direct correlation function for a hard-sphere (HS) fluid at the same density with suitably chosen effective particle diameter $\sigma$.Hansen and McDonald (2013); Barker and Henderson (1976) Recall that accurate expressions for the reference $c_{\mathrm{HS}}^{(2)}(r)$ exist and for a fluid of 1D hard-rods $c_{\mathrm{HS}}^{(2)}(r)$ is known exactly (see below). The RPA closure relation in Eq. (6) has been used extensively in the theory of simple and complex liquids.Hansen and McDonald (2013); Evans (1992, 1979); Likos (2001); Löwen (2002) The reliability of the corresponding OZ result for the structure factor (7) depends on the particular model system and the choice of reference potential.

III The mean-field DFT and the RPA

III.1 The non-uniform fluid

Consider now the fluid composed of particles interacting via the pair potential $v(r)$, split as in Eq. (5) and subject to an external potential $V({\bf r})$. The corresponding one-body density is $\rho({\bf r})$. Suppose too that we have an accurate DFT for the reference system, indicated by the subscript “0”, i.e. with particles interacting via the potential $v_{0}(r)$. The intrinsic Helmholtz free energy functional can be written as:Evans (1992)

[TABLE]

where $F_{0}[\rho]$ is the corresponding functional for the reference system. This exact expression is obtained from a thermodynamic integration ‘turning on’ the potential $v_{1}(r)$ between the particles via the integration parameter $\lambda$. The two-body density distribution function $\rho_{\lambda}^{(2)}({\bf r},{\bf r}^{\prime})$ is that for the system with interaction potential

[TABLE]

In deriving Eq. (III.1) one must impose an external potential, varying with $\lambda$, that ensures the equilibrium one-body density remains $\rho({\bf r})$ at each value of $\lambda$.Evans (1992) The standard MF DFT approximation is obtained by assuming that (i) $\rho_{\lambda}^{(2)}({\bf r},{\bf r}^{\prime})$ does not change much as $\lambda$ is varied from 0 to 1 and, more drastically, (ii)

[TABLE]

Eq. (10) clearly constitutes a mean-field treatment of the ‘perturbation’ $\lambda v_{1}(r)$. It follows that the MF DFT approximation for the excess (over ideal) free energy functional111The correlations neglected in (11) are incorporated into the functional $F_{\mathrm{corr}}[\rho]$ defined in Eq. (3.4.11) of Ref. Hansen and McDonald, 2013. is:

[TABLE]

Within DFTEvans (1979, 1992); Hansen and McDonald (2013) two functional derivatives of $-\beta F^{\mathrm{ex}}[\rho]$ with respect to the density yields the pair direct correlation function. From Eq. (11), and evaluating for a uniform (bulk) fluid $\rho({\bf r})=\rho_{b}$, we obtain the RPA approximation in Eq. (6), since

[TABLE]

The approximations inherent in Eq. (10) imply correlations are omitted and so one must be sceptical about the accuracy of the MF functional (11). Recognising that taking two functional derivatives of (11) yields the RPA (6, 7) for bulk correlation functions provides some useful insight into the status of the MF DFT. It is tempting then to argue that employing the MF DFT (11) should lead to results with similar accuracy to those obtained from the RPA for bulk liquids. However, this argument is at best misleading. In practical applications of (11), or any other DFT approximation, one works at the one-body level which requires only a single functional derivative. We explain and illustrate this below within the context of the test particle procedure for calculating $g(r)$.

III.2 The Percus test particle procedure

Percus provedPercus (1962) that one can determine the radial distribution function $g(r)$ by calculating the density profile $\rho({\bf r})=\rho(r)$ around a fixed particle that exerts an external potential $V({\bf r})\equiv v(r)$ on the fluid. Then the radial distribution function $g(r)=\rho(r)/\rho_{b}$. Within DFT, $\rho(r)$ is obtained by minimising the grand potential functional $\Omega[\rho]=F[\rho]-\int{\textrm{d}}{\bf r}[\mu-V({\bf r})]\rho({\bf r})$, where $\mu$ is the chemical potential. Using (11), the resulting Euler-Lagrange equation is

[TABLE]

where $\Lambda$ is the (irrelevant) thermal de-Broglie wavelength. For $r\to\infty$, away from the fixed test-particle, the density $\rho(r)\to\rho_{b}$, so within the approximation Eq. (11) we obtain the following relation between the chemical potential $\mu$ and the bulk density $\rho_{b}$:

[TABLE]

We make a functional Taylor expansion about the bulk density:

[TABLE]

where $H_{0}[\rho({\bf r})]$ denotes all higher order terms; these are $\sim{\cal O}([\rho-\rho_{b}]^{2})$ and higher. From Eqs. (12)–(15), we obtain:

[TABLE]

which eliminates $\mu$. Multiplying through by $-\beta$ and adding $(\rho(r)-\rho_{b})/\rho_{b}$ to both sides we obtain:

[TABLE]

This is the equation for the density profile, equivalent to $g(r)$, in the test particle treatment of Percus, as determined by the MF DFT in Eq. (11). If we set $v_{1}(r)=0$, then we obtain the following equation for the total correlation function $h_{0}(r)=g_{0}(r)-1$ of the reference system:

[TABLE]

Suppose we know the exact functional $F_{0}[\rho]$, and therefore $c_{0}^{(2)}(r)$, then comparison with Eq. (4) and use of the OZ equation (1) allows us to identify $-\beta H_{0}[\rho_{b}g_{0}(r)]$ as the exact bridge-function $B_{0}(r)$ of the reference system.

Returning to the full system, we see that Eq. (17) is an Ornstein-Zernike-like equation [see Eq. (1)] with the RPA closure (6) for the pair direct correlation function $c^{(2)}(r)$ inside the convolution integral, but with a different closure approximation for $c^{(2)}(r)$ [see Eq. (4)] appearing outside. Specifically, the sum of the first four terms on the right-hand side of (17) correspond formally to the exact expression for $c^{(2)}(r)$ but with the bridge function $B(r)$ replaced by $B_{0}(r)$.

The form of Eq. (17), which follows from the standard MF DFT functional (11), suggests that calculating the radial distribution function $g(r)$ via the test-particle route might yield results better than those given by $g_{\mathrm{RPA}}(r)$, obtained by solving the OZ equation together with the RPA closure (6), i.e. by Fourier inverting the RPA structure factor (7). This is evident for a fluid in which the pair potential $v(r)$ has a hard-core of diameter $\sigma$ since solving the Euler-Lagrange equation (13) guarantees the exact core condition $\rho(r)=0$, $r<\sigma$ is satisfied which is, of course, not the case for $g_{\mathrm{RPA}}(r)$, given by (7). In the latter $g_{\mathrm{RPA}}(r)$ is not identically zero inside the hard-core. We note that the application of any reasonable non-local DFT in the test particle procedure enforces the core condition. In the footnote222There are, of course, other theories that enforce the hard-core condition on $g(r)$. The well-known Mean-Spherical Approximation (MSA)Hansen and McDonald (2013); Lebowitz and Percus (1966) enforces this and sets $c^{(2)}(r)=-\beta v(r)$ outside the core. Perhaps more pertinent to our present discussion is the Optimized RPA (ORPA).Hansen and McDonald (2013); Andersen et al. (1972) The ORPA invokes the RPA closure (6) but, in addition, seeks to vary the perturbation potential $v_{1}(r)$ inside the hard core, where this is not uniquely defined, so that $g(r)=0$. This constraint is equivalent to requiring the functional derivative of the RPA free energy with respect to $v_{1}(r)$ to be zero inside the hard core. In our present approach, there is no attempt to minimize the free energy with respect to the perturbation potential, so for $r>\sigma$, $g(r)$ depends weakly on the choice of $v_{1}(r)$ inside the core. There is no reason to expect the MF DFT to perform as well as the ORPA – at least for a Lennard-Jones type fluid at high densities where the ORPA is known to be very accurate.Hansen and McDonald (2013); Andersen et al. (1972) For completeness, we also mention the Reference HNC (RHNC), see e.g. Ref. Hansen and McDonald, 2013, which improves upon the HNC integral equation by approximating the bridge function by that of a reference system, usually hard-spheres. Although the integral equation (17) that emerges from our MF DFT shares a feature of the RHNC, in that $B(r)$ is replaced by $B_{0}(r)$, this is implemented only in one part of the right hand side. There is no reason to expect the MF DFT to be as accurate as the RHNC. we mention briefly relationships to other theories of liquids.

This observation concerning the test-particle procedure has repercussions for more general external potentials $V({\bf r})$. Solving the corresponding Euler-Lagrange equation, based on the seemingly crude MF functional (11), yields equilibrium density profiles $\rho({\bf r})$ that are often very accurate – see Refs. Evans, 1992; Hansen and McDonald, 2013 and references therein. Investigation of hard/impenetrable potentials is once again illuminating. For a planar hard wall, with $V(z)=\infty$ for $z<0$, the density profile satisfies $\rho(z)=0$, $z<0$. Moreover, for a sensible choice of a (non-local) DFT for the reference system the profile will satisfyEvans (1992) the wall contact sum-rule: $\beta\rho(0^{+})=p(\rho_{b})$. The right hand side is the pressure of the bulk fluid, far from the hard wall, obtained from the bulk free energy $F[\rho_{b}]$ with $\rho_{b}=\rho(\infty)$.

In the following Sec. IV we compare the results for $g_{\mathrm{RPA}}(x)$ and $S_{\mathrm{RPA}}(k)$ (7) with those obtained from standard MF DFT (11), implemented within the Percus test particle prescription for a 1D fluid. The reference system is the hard-rod fluid for which the free energy functional is known exactly, and therefore its correlation functions and thermodynamics. Moreover, the pair correlation functions and the thermodynamics of the full system for uniform (bulk) densities are also known exactly. By considering this model fluid we can make a stringent examination of some of the basic approximations employed in classical DFT.

IV Illustrative results for a 1D fluid

In order to illustrate and support the observations made above, we consider a system of 1D hard-particles (rods on a line) with an additional attractive interaction between pairs of neighbouring rods. As mentioned above, this 1D system is chosen because we know the ingredients of the theory exactly and can therefore test carefully the accuracy of the various approximations. For a fluid of 1D rods with only nearest neighbour interactions $v(x)$, one finds the following exact expression for the structure factorPercus (1982)

[TABLE]

where $p$ is the pressure and $\mu(p)$ is the chemical potential which are known exactly.Takahashi (1942); Brader and Evans (2002) The above equation along with Eqs. (1) – (3) can be used to obtain the distribution functions of the uniform 1D system.

We consider 1D rods on a line interacting via the pair potential $v(x)=v_{0}(x)+v_{1}(x)$, where $x$ is the distance between the centres of a pair of particles – c.f. Eq. (5). $v_{0}(x)$ is the hard-rod potential

[TABLE]

and the attractive tail potential (also considered in Ref. Archer and Evans, 2013) is:

[TABLE]

The above potential is the 1D analogue of the Asakura-Oosawa potential for the effective colloid-colloid potential of hard-sphere colloids of diameter (length) $\sigma$ mixed with ideal polymers: $\sigma_{p}$ is the ‘length’ of polymer coils and $z_{p}$ is the fugacity of ideal polymers.Dijkstra et al. (1999); Archer and Evans (2013); Brader and Evans (2002)

For the hard-rods reference system “0” PercusPercus (1976, 1982) derived an exact expression for the excess Helmholtz free energy functional

[TABLE]

where the weighted density is

[TABLE]

Taking the functional derivative of (22) and adding the contribution from the ideal-gas free energy $F^{\mathrm{id}}[\rho]$, we obtain

[TABLE]

By differentiating again and using (12) the exact direct correlation function of a uniform fluid of hard-rods of length $\sigma$ and density $\rho_{b}$ is derived:

[TABLE]

where $\Theta$ is the Heaviside step function. This can be used to construct the RPA approximation (6) for the 1D system.

Comparison of the exact solution (solid line) with DFT results [calculated using the test particle equation (13) with (24) and $v(x)$ (dotted)] and results for $g_{\mathrm{RPA}}(x)$ and $S_{\mathrm{RPA}}(k)$ from the RPA approximation to the OZ equation (7) with $v_{1}(x)$ (dashed), is shown in Figs. 1, 2 and 3.333Note that for $\lambda=1$, $S_{\mathrm{MF}\lambda}(k)$, introduced in Eq. (41) of Ref. Archer and Evans, 2013, is identical to $S_{\mathrm{RPA}}(k)$ defined here. In Fig. 1 we fix the bulk density $\rho_{b}\sigma=0.5$ and the attraction range $\sigma_{p}=0.9\sigma$ while varying the attraction strength parameter $z_{p}$. In Fig. 2 we fix the attraction strength $\sigma z_{p}=1$ and range $\sigma_{p}=0.9\sigma$, while varying the bulk density $\rho_{b}\sigma$. Finally, in Fig. 3 we fix the attraction strength $\sigma z_{p}=2$ and the bulk density $\rho_{b}\sigma=0.5$ while varying the attraction range $\sigma_{p}$.

As expected, $g_{\mathrm{RPA}}(x)$ fails to satisfy the core condition $g(x)=0$ when $x<\sigma$, whereas the DFT enforces this. Outside the core of the hard-rod, both the DFT and RPA differ from the exact solution. However, the DFT results for $g(x)$ are much closer to the exact solution. As we increase the attraction strength or the density, both the DFT and RPA results deviate increasingly from the exact solution. $g_{\mathrm{RPA}}(x)$ displays only weakly damped oscillations. The same type of pattern is also observed in the sequence of structure factors displayed in Figs. 1(b) and 2(b). The RPA greatly overestimates the height of the first peak in $S(k)$ as $z_{p}\sigma$ (attractive well-depth) or $\rho_{b}\sigma$ are increased. It is important to note that when these parameters are increased beyond the values considered here, $S_{\mathrm{RPA}}(k)$ diverges at the first peak.Archer and Evans (2013) Note that for a given choice of pair potential the number density at which the divergence occurs is identical within both the RPA and MF DFT since the linearized form of both theories is identical; only the terms that are non-linear in $h(x)$ are different. However, the height of the first peak in the structure factor obtained from the MF DFT is always lower and closer to the exact value. It is also noteworthy that both the RPA and MF DFT capture well the behaviour in $S(k)$ at larger $k\sigma$.

The results in Fig. 3 show that for fixed attraction strength $z_{p}$ both the MF DFT and the RPA become less accurate as the range of the potential $\sigma_{p}$ increases. In particular the height of the first peak in $S(k)$ is overestimated. It appears that the integrated strength of the attraction, measured by the product $z_{p}\sigma_{p}$, is important in determining under what conditions both approximations are accurate.

V Discussion

In order to assess the physical content, and hence the usefulness, of an approximation for the excess free energy functional, it is not sufficient to take two functional derivatives, compute the bulk fluid pair direct correlation function $c^{(2)}(r)$, and then enquire how well this quantity performs when used to calculate the pair correlation function $h(r)$ via the OZ route, Eqs. (1)-(3). In most practical applications of DFT one is concerned with solving the Euler-Lagrange equation for the one-body density profiles and calculating the associated thermodynamic quantities (free energies) that result from minimizing a given (approximate) grand potential functional. This requires taking only a single functional derivative–not two. We have illustrated this point of view by focusing on the standard MF DFT, defined by (11), which treats the attractive part of the inter-particle potential via a simple MF approximation. Working in 1D and computing the density profile induced by a test particle exerting the potential (20, 21), we determined $g(x)$ and $S(k)$ from the MF DFT. These were compared with the exact results and with those obtained using (6) with the OZ equation. The latter corresponds to the usual RPA, given by (7). That the test particle route yields more accurate structure, within a DFT treatment, is not surprising. What is significant is that the MF DFT applied in the test particle situation performs much better for strong inter-particle attraction and for higher densities than the standard RPA. This is evident in Figs. 1–3.

That it is generally more appropriate to assess the performance of an approximate functional at the one-body rather than at the two-body level is known to the DFT community. The former requires only a single functional derivative with respect to density whereas the latter requires two. Naturally errors build up as further derivatives are taken. Here we are concentrating upon the efficacy of the particular functional (11). Why might the results, from what appears to be a crude approximation, be much better at the one-body density and free energy level than one might expect? We give three separate arguments: (i) as indicated in Sec. III, the Euler-Lagrange equation (17) for $\rho(r)$, when viewed as an integral equation for $g(r)$, implies a more sophisticated closure approximation than the RPA (6) implemented directly in the OZ equation (2), which leads to the RPA expression (7). (ii) Oettel Oettel (2005) discusses (11) in the context of a powerful and rather general reference functional approach for constructing approximate free energy functionals. By invoking the assumption that the bridge functional for the full system is well-approximated by that of the reference system and by considering expansions about the bulk density he argues that the reference functional approach predicts roughly MF behaviour for the density deviations (from bulk). He concludes that for adsorption problems, such as wetting and drying, a description based on the MF DFT (11), with an accurate reference functional $F_{0}[\rho]$, should capture all the essential physical features. OettelOettel (2005) also emphasizes that (11) has the advantage, over the more sophisticated reference functional approach, of satisfying identically the Gibbs adsorption equation and the wall-contact sum rule. (iii) Weeks and co-workersWeeks et al. (1997, 1998); Weeks (2002) introduced a local molecular field theory (LMF) that has proved to be highly successful in describing the structure and thermodynamics of a variety of non-uniform liquids. The derivationWeeks et al. (1998); Rodgers and Weeks (2008) of the LMF equation for the effective reference field $\phi_{R}({\bf r})$ starts with the Yvon-Born-Green equation and uses insightful arguments about the form of the conditional singlet densities $\rho({\bf r}|{\bf r}^{\prime})$ in the full and reference (mimic) systems. It does not employ concepts from DFT. Although LMF operates at the one-body level, like MF DFT, at first sight there does not appear to be a direct connection between the two approaches. This is not the case. Archer and EvansArcher and Evans (2013) showed that the LMF equation follows directly from the standard mean-field treatment of attractive interactions as embodied in MF DFT (11) and if one has access to the exact functional $F_{0}[\rho]$ for the same reference system the two theories are equivalent. We note that the derivation of the LMF equation and the relation to DFT is also discussed in the Supporting Information in a recent paper on solvation free energies.Remsing et al. (2016) In the 1D system described in Sec. IV the free energy functional of the reference (hard-rod) fluid is known exactly; it is given by the Percus result (22). It follows that our present results for $g(x)$ obtained using MF DFT and the test particle route are identical to those that would emerge from LMF using hard-rods as the reference system. More generally, in three dimensions, LMF with a hard-sphere reference system would lead to the same $g(r)$ as MF DFT using the test particle route and a very accurate hard-sphere functional for $F_{0}[\rho]$. Given the success of LMF for a wide variety of fluids, one might argue, albeit circuitously, that the physical arguments and plausible approximations that lead to LMFWeeks et al. (1998); Rodgers and Weeks (2008); Remsing et al. (2016) provide an alternative justification as to why MF DFT might perform better than one might expect. Of course, there is a caveat. The justification for LMF relies upon the judicious choice of reference fluid, described by $v_{0}(r)$, so that the LMF equation used to treat $v_{1}(r)$, the longer ranged part of the pair potential, captures the essential physics for a given model fluid.Rodgers and Weeks (2008); Remsing et al. (2016) Although the MF DFT treatment of $v_{1}(r)$ is formally equivalent to that of the LMF, the limitation and drawback of the former is finding an accurate free energy functional for the reference fluid.Archer and Evans (2013)

As a final note of caution on assessing the quality of a DFT on the basis of what two derivatives of the free energy functional yields for the (bulk) pair direct correlation function, one should also recall the following functional:

[TABLE]

This is the well-known Ramakrishnan-Yousouf functional,Ramakrishnan and Yussouff (1979) constructed to yield the ‘exact’ pair direct correlation function, $c(r)=c^{(2)}_{\mathrm{exact}}(r)$, with $c^{(2)}_{\mathrm{exact}}(r)$ calculated at the relevant bulk density $\rho_{b}$.444Note that minimizing the grand potential corresponding to (26) within the test particle procedure leads to the HNC closure for the uniform fluid, provided the direct correlation function is determined self-consistently using the OZ equationEvans (1992); Hansen and McDonald (2013) However, the functional (26) has significant weaknesses. For example, it is unable to describe wetting or drying phenomenaEvans (1992); Evans et al. (1983) at substrates nor critical adsorptionEvans and Marconi (1986) owing to the fact that it is only quadratic in the density deviation.

Acknowledgements

B.C. is supported by an EPSRC studentship and R.E. by a Leverhulme Emeritus Fellowship: EM-2016-031.