pubs.acs.org/Langmuir © 2010 American Chemical Society
Free Energy Density Functional for Adsorption of Fluids in Nanopores Shiqi Zhou* State Key Laboratory of Powder Metallurgy, Central South University, Changsha, Hunan 410083, China, and School of Physics Science and Technology, Central South University, Changsha, Hunan 410083, China Received June 8, 2010. Revised Manuscript Received August 5, 2010 A classical free energy density functional, which is isomorphic to a usual effective hard sphere model þ mean field approximation for tail contribution, is proposed for treatment of real fluids in inhomogeneous states. In the framework of the classical density functional theory (DFT), the present functional is applied to two representative model fluids, namely, a Lennard-Jones fluid and a hard core attractive Yukawa fluid, subject to influence of various external fields. A comprehensive comparison with simulation results and a detailed analysis show that the present functional holds simultaneously all of the desirable properties inherent in an excellent functional, such as high accuracy, computational simplicity, consistency with a hard wall sum rule, nonrecourse to use of adjustable parameter(s) and weighted densities, reproduction of bulk second-order direct correlation function (DCF) in bulk limit, and applicability to subcritical fluid phenomena.
I. Introduction In the past several decades, classical density functional theory (DFT)1 has drawn the attention of many researchers working in fields of liquid physics, polymer statistical physics, colloid science, and others. In principle, only if an excess Helmholtz free energy density functional Fex[F] is known2 can all liquid problems be tackled in the framework of DFT. Unfortunately, the exact expression for a long-sought “divine” Fex[F]3 is never known; approximation therefore has to be introduced. As in the bulk state, the van der Waals picture4 also plays an important role in formulating an approximation for Fex[F] for real interparticle potentials usually comprising a highly steep repulsion at close separation and an attractive tail at large separation; thus, one two-step procedure is adopted in numerous efforts to approximate Fex[F]. First of all, a great deal of effort has been devoted to approximating Fex-hs[F], a hard sphere counterpart of Fex[F], and often employed as reference; and second, numerous activities are invested in making approximations for the tail functional Fex-tail[F]. Then, a simple sum of Fex-hs[F], for an effective hard sphere fluid, and Fex-tail[F] leads to the desired Fex[F]. Although several successful density functional approximations (DFA) have been advanced1,5-10 for Fex-hs[F], our knowledge about Fex-tail[F] is unfortunately still poor. Until now, a mean field approximation (MFA) exclusively used in the early years of development of DFT1,2,4 is still often employed in the literature. Although a lot of effort has been *Corresponding author. E-mail:
[email protected]. (1) Henderson, D. Fundamentals of Inhomogeneous Fluids; Marcel Dekker: New York, 1992. (2) Evans, R. Adv. Phys. 1979, 28, 143. (3) Mattsson, A. E. Science 2002, 298, 759. (4) Zhou, S.; Solana, J. R. Chem. Rev. 2009, 109, 2829. (5) Tarazona, P. Phys. Rev. A 1985, 31, 2672. (6) Curtin, W. A.; Ashcroft, N. W. Phys. Rev. A 1985, 32, 2909. (7) Rosenfeld, Y. Phys. Rev. Lett. 1989, 63, 980. (8) Kierlik, E.; Rosinberg, M. L. Phys. Rev. A 1990, 42, 3382. (9) Khein, A.; Ashcroft, N. W. Phys. Rev. Lett. 1997, 78, 3346. (10) Tarazona, P. Phys. Rev. Lett. 2000, 84, 694. (11) Walton, J. P. R. B.; Quirke, N. Chem. Phys. Lett. 1986, 129, 382. (12) Ustinov, E. A.; Do, D. D. Langmuir 2003, 19, 8349. (13) (a) Sweatman, M. B. Phys. Rev. E 2001, 63, 031102. (b) Tang, Y.; Wu, J. J. Chem. Phys. 2003, 119, 7388. (14) (a) Zhou, S. Commun. Theor. Phys. 2003, 40, 721. (b) Zhou, S.; Jamnik, A. J. Chem. Phys. 2005, 122, 064503. (c) Zhou, S.; Jamnik, A. J. Chem. Phys. 2005, 123, 124708. (d) Zhou, S.; Jamnik, A. Phys. Rev. E 2006, 73, 011202.
Langmuir 2010, 26(22), 17037–17047
invested in modifying11,12 the MFA, or proposing13 a new formalism type of WDA or other route14 for Fex-tail[F], there are still undesirable properties inherent in these promising approximations, as will be detailed next. The aim of this paper is to propose a novel formalism for Fex[F] to overcome these undesirable points. The organization of this paper is as follows. In section II, we set out the present an approximation for Fex[F]. In section III, we apply the present Fex[F] in the framework of DFT to two kinds of model potentials in the inhomogeneous state, i.e., a LennardJones (LJ) potential, usually used to model neutral atomic fluids and approximate spherical simple molecular fluids, and a hard core attractive Yukawa (HCAY) potential, frequently engaged in modeling charged colloidal particles immersed in electrolyte solution; an exhaustive comparison among the present free energy functional approximation, relevant simulation results, and other existing theories is also presented in section III. Then, the paper is ended with an in-depth discussion on the differences and advantages of the present approximation compared to existing approximations, and several concluding remarks are in section IV.
II. Methodology II. A. Free Energy Density Functional. The MFA specifies Fex-tail[F] as follows: Fex-tail ½F ¼
1 2
Z
dr0
Z
dr Fðr0 ÞFðrÞuatt ðjr0 - rjÞ
ð1Þ
Here, F(r) is particle density in space position r, and uatt(r) is the tail contribution to a full pair potential u(r) under consideration. Then, a simple sum of Fex-hs[F; d] for the effective hard sphere fluid with an effective hard sphere diameter d, differing from a size parameter σ of u(r) under consideration, and Fex-tail[F] leads to an expression for the perused functional, henceforth also referred to as MFA and denoted by Fex-mfa[F] Z Z 1 dr0 dr Fðr0 ÞFðrÞuatt ðjr0 - rjÞ Fex-mfa ½F ¼ Fex-hs ½F; d þ 2 ð2Þ Drawbacks inherent in the MFA are that (i) a simple sum of Fex-tail[F] and Fex-hs[F; d] ignores correlation effects originating from a coupling between a hard core repulsion and the tail contribution, which is reflected in the fact that a bulk second-order direct correlation
Published on Web 10/18/2010
DOI: 10.1021/la102341a
17037
Article
Zhou
function (DCF) C(2) 0-mfa(r, Fb) (Fb is density of a coexistence bulk fluid with which the system under consideration is in equilibrium, henceforth referred to as coexistence bulk density), obtained by performing a functional derivative on the sum and then taking a bulk limit, is not in agreement with that solved from an OrnsteinZernike (OZ) integral equation theory (IET); (ii) there exists no unified approach to specify consistently uatt(r) and d entering Fex-hs[F; d]; (iii) the bulk equation of state (EOS) resulting from the sum in the bulk limit is not very accurate. The bulk second-order DCF C(2) 0-mfa(r, Fb) corresponding to the MFA Fex-mfa[F] is obtainable by a basic relationship ð2Þ
F f Fb
ð3Þ
Here, β = 1/kBT with kB the Boltzmann constant and T the absolute temperature, C(2) 0-eff-hs(r, Fb) is the bulk second-order DCF for the effective hard sphere fluid whose concrete expression depends on the expression used for Fex-hs[F; d]. As the plain MFA is only a rough approximation, unavoidably the relevant C(2) 0-mfa(r, Fb) is unfavorably influenced by the uncontrollable approximations and is not expected to be accurate. Throughout the text, we denote a bulk second-order DCF obtained from other more accurate (2) routes by C(2) 0 (r, Fb); then, C0-mfa(r, Fb) is in all probability not (2) (2) equal to C0 (r, Fb). After substituting C(2) 0-mfa(r, Fb) by C0 (r, Fb) in eq 3, one new tail contribution, denoted by u0att(r), is obtained as follows: ð2Þ
ð2Þ
βu0 att ðrÞ ¼ C 0-eff-hs ðr, Fb Þ - C 0 ðr, Fb Þ
ð4Þ
Replacing uatt(r) in eq 2 by u0att(r), one acquires the present approximation for Fex[F] as follows:
1 þ 2
Z
βF ex ½F ¼ βF ex-hs ½F; d Z h ð2Þ dr0 dr Fðr0 ÞFðrÞ C0-eff-hs jr0 - rj, Fb ð2Þ
- C0
jr - r0 j, Fb
i
ð5Þ
Possibly, one may want to know what we obtain from the above substitution. A preliminary examination immediately draws forth two apparent advantages of the present eq 5 over the traditional MFA. First, the present functional reproduces the bulk second(2) order DCF C(2) 0 (r, Fb) in the bulk limit. Since C0 (r, Fb) is more (2) accurate than the MFA C0-mfa(r, Fb), and it is known that a free energy functional approximation improves by an in-line more accurate bulk second-order DCF, this provides corroborating evidence of the potentially improved performance of the present functional compared to the MFA. Second, in the new approximation (eq 5), one is no longer perplexed by the self-consistency with which one specifies both parameters d and uatt(r) in the MFA; what one has to designate is simply d, and the self-consistency problem no longer exists. Our procedure, used to determine d in the present eq 5, is elaborated below. In the bulk limit, the present βFex[F] gives an expression for a bulk excess Helmholtz free energy per particle fex(Fb) for the full pair potential u(r) βf ex ðFb Þ ¼ βf ex-hs ðFb ; dÞ þ
1 F 2 b
Z
h i ð2Þ ð2Þ dr C0-eff-hs r, Fb - C0 r, Fb
ð6Þ Here, fex-hs(Fb; d) is an effective hard sphere counterpart of fex(Fb). A nonuniform first-order DCF C(1)(r, [F]) for the full pair potential u(r) is connected with βFex[F] by a functional relationship -δβFex ½F C ð1Þ r, ½F ¼ δFðrÞ 17038 DOI: 10.1021/la102341a
i ð2Þ - C 0-eff-hs jr - r0 j, Fb Fðr0 Þ
ð7Þ
ð8Þ
where C(1) hs (r, [F]; d), an effective hard sphere counterpart of eff-hs C(1)(r, [F]), must be consistent with the engaged βFex [F], i.e. ð1Þ
C hs ðr, ½F; dÞ ¼
2
- δ βF ex-mfa ½F δFðrÞ δFðr0 Þ ð2Þ ¼ C 0-eff-hs jr - r0 j, Fb - βuatt jr - r0 j
C0 - mfa ðjr - r0 j, Fb Þ ¼ lim
then a combination of eqs 5 and 7 immediately leads to Z h ð1Þ ð2Þ C ð1Þ r, ½F ¼ Chs r, ½F; d þ dr0 C 0 jr - r0 j, Fb
-δβFex-hs ½F; d δFðrÞ
ð9Þ
Performing a bulk limit on eq 8 generates an expression for (1) C(1) 0 (Fb), a bulk counterpart of C (r, [F]), Z h i ð1Þ ð1Þ ð2Þ ð2Þ C0 ðFb Þ ¼ C0hs ðFb ; dÞ þ Fb dr C0 r, Fb - C0-eff-hs r, Fb ð10Þ (1) where C(1) 0hs(Fb; d), a bulk counterpart of Chs (r, [F]; d), is mathe(r, [F]; d). The expressions for both matically a bulk limit of C(1) hs C(1) 0 (Fb) and βfex(Fb) directly lead to a bulk compressibility factor Zf inherent in the present functional scheme Z h 1 ð1Þ ð2Þ Zf ¼ 1 - C0hs ðFb ; dÞ - βfex-hs ðFb ; dÞ þ Fb dr C0-eff-hs r, Fb 2 i ð2Þ ð11Þ - C0 r, Fb
By forcing Zf relevant to the present functional equal to an externally supplied bulk compressibility factor, d can be uniquely identified.
II. B. Local Self-Consistent OZ Integral Equation Theory.
Another important input parameter, C(2) 0 (r, Fb), is obtainable by solving a bulk OZ IET,15a which is described briefly below. The OZ IET is defined by an OZ IE, and a closure equation, the former, is expressed as Z ð2Þ ð2Þ hðrÞ ¼ C0 ðr, Fb Þ þ Fb dr0 hðr0 ÞC 0 jr - r0 j, Fb ð12Þ and the latter is given as hðrÞ þ 1 ¼ expf- βuðrÞ þ γ þ BðrÞg
ð13Þ
Fb) are total and where h(r) = g(r) - 1 and γ(r) = h(r) indirect correlation functions, respectively, and g(r) is a radial distribution function (rdf); B(r) is a so-called bridge function which refers to special multidimensional integrals and has to be approximated. Different OZ IET versions are distinguished by different approximations used for B(r), and in this paper, we use a local version of a recently proposed global self-consistent OZ IET16a to acquire C(2) 0 (r, Fb). Relevant details are summarized below. In ref 16a, a new hard sphere bridge function, denoted by BHS(r/σhs, Fhsσhs3) (σHS and FHS being hard sphere diameter and hard sphere number density, respectively), is proposed; B(r) of the present interest is then approximated by BHS with σHS being substituted by an effective hard sphere diameter deff, which is different from the size parameter σ and the effective hard sphere diameter d entering the present functional. One thus has C(2) 0 (r;
B r=σ, Fb σ 3 ¼ BHS r=d eff , Fb d eff 3
ð14Þ
(15) (a) Hansen, J. P.; McDonald, I. R. Theory of simple liquids, 2nd ed.; Academic Press: London, 1986; (b) Wertheim, M. S. J. Math. Phys. 1964, 5, 643. (c) Thiele, E. J. Chem. Phys. 1963, 39, 474. (16) (a) Zhou, S. Phys. Rev. E 2009, 79, 011126. (b) Zhou, S. Theor. Chim. Acta 2007, 117, 555. (17) Henderson, J. R. In Fundamentals of Inhomogeneous Fluids; Henderson, D., Ed.; M. Dekker: New York, 1992.
Langmuir 2010, 26(22), 17037–17047
Zhou
Article
To determine the adjustable parameter deff entering the approximation for B(r), a global self-consistent condition is used in ref 16a, which imposes an equality of pressure between the virial route and the fluctuation route. As calculation of the pressure from the fluctuation route is involved with a thermodynamic integration from zero density to the density of interest, the global procedure is thus involved with calculations for many state points even if one is only interested in one state point; moreover, when a subcritical region is being dealt with, the global procedure will incur complexity of numerical solution, as a bulk vapor-liquid coexistence region then has to be avoided to reach a state point situated at the liquid side. As we will discuss, the present functional only needs as an input the information of the coexistence bulk state, and desirably is not involved with information of many bulk state points with varying coexistent bulk densities; thus, one perhaps tends to use a “local” procedure16b to determine deff since the “local” procedure is free of the two disadvantages relevant to the “global” procedure. In the “local” procedure,16b an equality of a virial reduced isothermal compressibility χ0/χTv and a reduced isothermal compressibility χ0/χTc throughout the fluctuation route is used to specify deff, namely, the following equation is imposed χ0 χ ¼ c0 χvT χT
ð15Þ
The virial reduced isothermal compressibility χ0/χTv is calculated by a finite difference method as below χ0 ∂P ¼ β χvT ∂Fb βPv Fb þ dFb , T , d eff - βPv Fb , T , d eff ¼ dFb
II. C. Application in Framework of Density Functional Theory. As for βFex-hs[F; d], we will use a fundamental measure
functional (FMF).7,8 The Kierlik and Rosinberg FMF8 is found18a to actually be equivalent to a Rosenfeld FMF version7 and simultaneously has no recourse for vectorial weighted densities. However, in this paper, we will use a Mansoori-Carnahan-StarlingLeland (MCSL)-augmented18b FMF version. An analytical expression consistent with the hard sphere MCSL equation of (2) state should be used for C(1) 0hs(Fb; d); as for C0-eff-hs(r, Fb), an analytical expression stemming from a Percus-Yevick approximation15b,c for the bulk hard sphere fluid is used for simplicity. Grand potential Ω[F] of the inhomogeneous fluid is connected with Fex[F] via a Legendre transform Z ð20Þ Ω½F ¼ F½F þ dr FðrÞ jext ðrÞ - μ where μ is chemical potential of the coexistence bulk fluid and jext(r) is an external potential responsible for the generation of F(r). One notes that F [F] is decoupled into an idea part Fid[F] and an excess part Fex[F] F½F ¼ Fid ½F þ Fex ½F The ideal part is exactly given by Z Fid ½F ¼ kT dr FðrÞ ln λ3 FðrÞ -1
ð21Þ
ð22Þ
v
ð16Þ
Throughout this paper, the size parameter σ and an energy parameter ε entering the potential function of interest are used as length and energy units, respectively; thus, the bulk density Fb is normalized to Fb* = Fbσ3, and the temperature T is normalized to T* = kT/ε. |dF*| b is a conveniently small density increment (typically 10-4), and dFb* = dFbσ3. The virial pressure Pv in eq 16 is calculated from the virial route Z βPv βF duðrÞ gðrÞ ð17Þ ¼ 1- b dr 4πr3 dr Fb 6 0 whereas the reduced isothermal compressibility χ0/χTc throughout the fluctuation route reads Z χ0 ð2Þ dr C0 r, Fb ¼ 1 F ð18Þ b c χT Just as is clearly shown in eqs 15-18, numerical implementation of the “local” procedure only refers to information of the state point of interest and the most adjacent point; no thermodynamic integration is needed. Below the critical temperature, there is an instability region in the phase space where the OZ IET is no longer solvable physically or numerically; to prevent the most adjacent point from entering into the instability region and make the numerical code run safely, dFb* should be chosen to be positive when the state point is adjacent to the liquid side of the vapor-liquid coexistence curve and negative when the state point is adjacent to the vapor side of the vapor-liquid coexistence curve, whereas for the supercritical state, dFb* can be positive or negative, since there is no instability region available for phase space above the critical temperature. Upon numerical solution of the OZ IET, a virial compressibility factor Zv is obtained by integrating the rdf g(r), as shown in eq 17 Zv ¼
βPv Fb
ð19Þ
Throughout the text, as a general rule we will determine d by equating Zv to Zf if not particularly highlighted. Langmuir 2010, 26(22), 17037–17047
Here, λ is a thermal de Broglie wavelength. Fex[F], representative of many-body correlation effects originating from the interparticle interactions and never known exactly as discussed, is approximately given by eq 5 in this paper. The local structural property, represented by the density profile F(r), is immediately calculable by a density profile equation which, denoted by eq 23 and stemming from minimizing Ω[F], is given as n o ð1Þ FðrÞ ¼ Fb exp - βjext ðrÞ þ C ð1Þ ðr; ½FÞ - C0 ðFb Þ ð23Þ Given that C(1)(r; [F]) and C(1) 0 (Fb) are approximated by eqs 8 and 10, respectively, F(r) is amenable to numerical calculation for any external field. Once the density profile is available, the grand potential of the inhomogeneous system is obtainable by eq 20, in which the chemical potential μ is decoupled into ideal contribution μid and excess contribution μex, and μid is analytically known as μid = ln(λ3Fb) and μex = -kTC(1) 0 (Fb). Once the two basic quantities, F(r) and Ω[F], are obtained, other derivatives are readily computable, as will be detailed next.
III. Results First, the present functional will be applied to the HCAY fluid under the influence of external fields, which are, respectively, due to a large hard sphere particle and two planar hard walls. The HCAY potential function reads uðrÞ ¼ ¥ r σ
ð24Þ
Here, κ* is a reduced screening parameter measuring the range of the Yukawa potential; as pointed out above, throughout this paper, σ and ε are, respectively, the size parameter and energy parameter entering the potential function of interest. jext(r) derived from the large hard sphere particle reads jext ðrÞ ¼ ¥ r < R ¼ 0 r>R
ð25Þ
(18) (a) Phan, S.; Kierlik, E.; Rosinberg, M. L.; Bildstein, B.; Kahl, G. Phys. Rev. E 1993, 48, 618. (b) Zhou, S. Commun. Theor. Phys. (accepted).
DOI: 10.1021/la102341a
17039
Article
Zhou
whereas the two planar hard walls generate jext(r) given as 0 0:5σ < z < H - 0:5σ ð26Þ φext ðrÞ ¼ otherwise In this paper, the present functional is used to calculate the density profiles of the HCAY fluid in the two external field cases, which had been simulated by a grand canonical Monte Carlo (GCMC) simulation, and these simulation results14b are used to test the present functional approximation. One notes that the density profiles due to the large hard sphere particle were also calculated by a previous nonuniform fifth-order-PTS-TPT;19 in this paper, we additionally calculate the HCAY fluid density profiles due to the two planar hard walls by using the nonuniform fifth-order-PTS-TPT. All of these results are displayed in Supporting Information Figures 1S-4S for the case of the large hard sphere particle and in Figures 1-4 for the case of the two planar hard walls. Before discussing the HCAY fluid density profiles, we first test whether the present functional satisfies a single hard wall sum rule,20 which declares that reduced contact density F(0.5σ)/Fb of the fluid near a single planar hard wall is equal to the bulk compressibility factor used as input; in the present paper, it is Zv. One notes that the single planar hard wall external field is a limiting case of the external field eq 26 on the condition that H is so large that the density profiles near each of the two hard walls do not interfere with each other. Table 1 presents the calculated results for the reduced contact density respectively from the present scheme and the sum rule for the case of T* = 1.3, κ* = 1.8, and the coexistence bulk density F* varying over a range. From Table 1, one finds that the sum rule is met with an accuracy of percent relative error (PRE) lower than 0.23; the presence of small residual error, we think, stems from an inaccuracy of the bulk second-order DCF used as input, and we explain this point below. The sum rule is a deep statistical mechanical property that need not be satisfied in any approximate theories;20 by construction, it holds exactly in weighted-density approximations (WDA),1 but in other DFAs, it is not necessarily met. However, it is always desirable that the sum rule is met in any DFAs, since the contact reduced density will be predicted with an accuracy with which the input compressibility factor approaches the real bulk compressibility factor once the sum rule is met. In WDA, C(2) 0 (r; Fb) is used to construct weighting functions, and the weighting functions are always made to meet a normalization condition; as a result, any inaccuracy in C(2) 0 (r; Fb) does not undermine the normalization condition, and therefore the sum rule is always met. In the present functional, C(2) 0 (r; Fb) enters into the functional itself in a more fundamental way, and consequently, the properties of C(2) 0 (r; Fb) truly influence whether the sum rule is met. In the present paper, the sum rule is able to be met as much as possible by guaranteeing self-consistency between C(2) 0 (r; Fb) and the compressibility factor; this is why we use the local self-consistency condition16b to close the OZ IE. It is found that, if the bridge function does not enable the local selfconsistency condition,16b the sum rule will only be met poorly. It should be pointed out that the “local” procedure implicit in eqs 15-18 ignores the dependence of deff on the density argument, and the influence of overlooking this on the final performance of the OZ IET is a complex issue. Although it is generally considered that the “local” procedure is somewhat vulnerable, in most cases the “local” OZ IET in fact predicts structure functions in very (19) Zhou, S. Theor. Chim. Acta 2009, 124, 279. (20) Henderson, J. R. In Fundamentals of Inhomogeneous Fluids; Henderson, D., Ed.; M. Dekker: New York, 1992.
17040 DOI: 10.1021/la102341a
Figure 1. Density profiles of the HCAY fluid of κ* = 1.8 at supercritical temperature of T*=1.25 and, with varying coexistent bulk densities as shown, the responsible external field is due to two planar hard walls situated at z=0 and z=4σ, respectively. Symbols are for simulation data from ref 14b; dashed lines are calculated in this paper using the nonuniform fifth-order-PTS-TPT,19 whereas solid lines denote the results based on the present functional.
good agreement with those predicted by the “global” OZ IET; as a result, the influence of overlooking is most often negligible. However, such overlooking may make the self-consistency condition only partially met in the “local” OZ IET; particularly in the high-density domain, the resulting partial violation of the selfconsistency condition may tend to become somewhat evident. This is why the PRE documented in Table 1 tends to rise when one enters the high-density domain. In practice, the following situation often appears: the empirical equation of state (EOS) is available Langmuir 2010, 26(22), 17037–17047
Zhou
Figure 2. Same as Figure 1 but for κ* = 3, T* = 0.76, and the coexistent bulk densities as shown.
for the potential function of interest, and one wishes to incorporate the accurate EOS information into the present functional approximation, particularly using the EOS to specify d. For this purpose, we propose a viable procedure as below. First, the virial compressibility factor is substituted by the new one, and then C(2) 0 (r; Fb) for r
