L a n g m u i r 1995,11,4183-4184
The Method of Indeterminate Lagrange Multipliers in Nonlocal Density Functional Theory Alexander V. Neimark’ Department of Chemical Engineering, Yale University, 9 Hillhouse Avenue, New Haven, Connecticut 06520-8286 Received March 23, 1995. I n Final Form: July 31, 1995
4183
and references therein. It was shown that NLDFT gives density profiles almost identical to those generated by means of grand canonical Monte Carlo (GCMC) simulations,14J7~18J9~22 which are regarded as reference computer experiments. This supports the reliability of NLDFT results and its potential applicability to diverse experimental systems. In this note we propose an effective method for calculating the density profiles of adsorbates confined in pores of arbitrary geometry. The method is based on the implementation of indeterminate Lagrange multipliers (ILM) for minimizing the GPF. This method is not restricted to the adsorption systems under consideration; rather it can be employed in any heterogeneous system. In the following, the modern conventional formulation of NLDFT, similar to that employed in, e.g., refs 15 and 19, is used, and the foundations ofthe basis equations are not discussed. Consider a fluid confined in a volume M surrounded by solid walls, Le., an adsorbate in a pore, at a given chemical potential ,u and temperature T . The fluid-solid interactions are modeled by a particular spatially varying external potential Vext. The fluid is regarded as a hard sphere fluid with a given attractive potential Qatt. The GPF for this system is represented in a standard form, as a functional with respect to the spatially varying fluid density p ( r ) : 1 6
Density functional theory (DFT) is regarded as one of the most powerful methods of accounting for real intermolecular interactions and currently finds many important applications related to equilibrium, phase transitions, and dynamics in heterogeneous systems; for a recent review, see, e.g., ref 1. DFT, proposed by Hohenberg and Kohn,2 Kohn and S h a ~ and , ~ Mermin4 30 years ago, implies in general that the equilibrium spatial density distributions in inhomogeneous systems can be derived by means of minimization of the grand potential functional (GPF) a t given external and internal potentials with respect to the variable singe-particle densities. Even the simplest local version of DFT (LDFT),5which does not take into account short-ranged spatial correlations in the repulsive part of the GPF, gives quite reasonable results in many practical situations; see, for example, the studies of capillary condensation in pores by Evans et aL6and the calculations ofpore size distributions in porous solids from adsorption isotherms worked out by Seaton et al.’ The most advanced nonlocal version of DFT (NLDFT) was first developed by Tarazona,Ewho introduced the smoothed (or weighted) densities, and thereafter modified by Evans et al.,9 Rosenfeld,lo Denton and Ashcroft,ll Kierlik and Rosinberg,12and Patra and Ghosh.13 NLDFT has recently been successfully employed for analyses of a number of quite complicated problems in interfacial and Here the first term on the right-hand side represents the adsorption equilibrium, phase transitions, complex fluids, ideal free energy functional, the second term is the excess etc.; as the most recent research papers, see, e.g., those hard sphere Helmholtz free energy functional, where f e x of Kierlik et al.,14 Tang et al.,15 Lastoskie et al.,16J7 is the excess Helmholtz free energy per molecule. fex Sokolovski and Fisher,18Babluena and Gubbins,lg Patra and Ghosh,20,21 Olivier et a1.,22Chmiel et al.,23R o ~ e n f e l d , ~ ~ depends on the smoothed density, p(r), which is defined as
’ E-mail:
[email protected]. (1)Evans, R. InFundamentaZs oflnhomogeneous Fluids; Henderson, D., Ed.; Marcel Dekker: New York, 1992; Chapter 5. (2) Hohenberg, P.; Kohn, W. Phys. Rev. B 1964, 136, 864. (3) Kohn, W.; Shaw, L. J.Phys. Rev. A 1965, 140, 1133. (4) Mermin, N. D. Phys. Reu. A. 1965, 137, 1441. (5) Rowlinson, J. S.; Widom, B. Molecular Theory o f Capillarity; Oxford University: London, 1982. (6) Evans, R.; Marini Bettolo Marconi, U.; Tarazona, P. J. Chem. Soc., Faraday Trans. 2 1986, 82, 1763. (7) Seaton, N. A,; Walton, J. P. R. B.; Quirke, N. Carbon 1989,27, 853. ( 8 ) Tarazona, P. Phys. Rev. A 1986, 31, 2672. (9) Tarazona, P.; Marini Bettolo Marconi, U.; Evans, R. Mol. Phys. 1987, 60, 573. (10) Rosenfeld, Y. Phys. Rev. Lett. 1989, 63, 980. (11)Denton, A. R.; Ashcroft, N. W. Phys. Rev. A. 1991, 44, 8242. (12) Kierlik, E.; Rosinberg, M. L. Phys. Rev. A. 1991, 44, 5025. (13) Patra, C. N.; Ghosh, S. K. Phys. Rev. E 1993,47, 4088; 1993, 48. 1154. (14) Kierlik, E.; Rosinberg, M. L.; Finn, J. E.; Monson, P. A. Mol. Phys. 1992, 75, 1435. (15) Tang, 2.; Scriven, L. E.; Davis, T. H. J . Chem. Phys. 1992,97, 494, 9258. (16) Lastoskie, C.; Gubbins, K. E.; Quirke, N. J. Phys. Chem. 1993, 97, 4786. (17) Lastoskie, C.; Gubbins, K. E.; Quirke N. Langmuir 1993,9,2693. (18) Sokolovski, S.; Fisher, J. J. Chem. Soc., Faraday Trans. 1993, 85, 789. (19) Babluena, P. B.; Gubbins, K. E. Langmuir 1993, 9, 1801. (20) Patra, C. N.; Ghosh, S. K. Phys. Rev. E. 1994,49, 2826,1994, 50, 5123. (21)Patra, C. N.; Ghosh, S. K. J . Chem. Phys. 1994, 100, 5219. ~~
0743-7463/95/2411-4183$09.00/0
It is assumed that the function f&) can be explicitly calculated as well as its derivative fex’, for example, by using the Carnahan-Sterling equation of state. The different versions of NLDFT use different expressions for the weighting function w.7J0’11J2 In the most popular approach of Tarazona,Ewhich is in use in the majority of current publications, e.g., refs 16, 17, and 20, the secondorder expansion of the weighting function is employed, which leads to an algebraic equation with respect to the smoothed density. The problem of finding an equilibrium density distribution p(r)is the problem of minimizing the GPF (1)with respect to p(r)subject t o the integral relationship (2).The standard method employed first by TarazonaEis based on the Euler equation (22) Olivier, J. P.; Conklin, W. B.; v. Szombathely, M. I n Characterization of Porous Solids 111;Rouquerol, J., Rodriquez-Reinoso, F., Sing, K. S. W., Unger, K. K., Eds.; Elsevier: Amsterdam, 1994; p 81. (23) Chmiel, G.; Patrykiejew, A.; Rzysko, W.; Sokolowski, S. Mol. Phys. 1994, 83, 19. (24) Rosenfeld, Y. J . Phys. Chem. 1996, 99, 2857.
1995 American Chemical Society
Notes
4184 Langmuir, VoL. 11, No. 10, 1995
It is worth noting that in the case of the LDFT, wilr r’l,p(r))= d(r - r’) and p(r)= p(r), eq 8 degenerates into
i(r)= - g(r) fex’(p(r))/kT (9) and, correspondingly, eq 7 degenerates into the standard Euler equation for LDFT6 where the functional derivative of p(r)with respect to p(r) should be calculated from the relationship in eq 2. This leads to tedious calculations in the case when dp(r ) l d p ( r ) cannot be expressed explicitly. We introduce the Lagrange functional:
L(e,Q)= R
+ kT JM
d r ji(r){Q(r)-
IMdr’ g(r’) o A I r r’l, ,Ob)))
=R
+ kT IMd r Xr) Q(r)-
hT
JMjM d r dr‘ L(r)e(r’)w(lr - r’l,Q(z-1)
(4)
Here, I(r)is a function of the ILM. The implementation of the ILM function A(r)for minimizing the GPF allows one to avoid the functional derivative under the integral sign in the Euler equation. This is a key point in our approach. The Euler equations for the Lagrange functional (4) are
and
Here, the derivative au(lr - r’l ,p(r))/ap is not a functional derivative but rather a partial derivative of the weighting function w with respect to the smoothed density p which can be calculated explicitly as a function of /r- r‘l and p. These two equations (eqs 5 and 61, together with eq 2, constitute a closed system of equations for determination of three functions, the density p(r), the smoothed density p(r),and the ILM function k ) .Equations 5 and 6 can be rewritten in a form convenient for computations:
and
1 lnEe(r)le,l = - p ( V e x t ( r+) fex(e(r)) + efex’(e(r)) + JM
dr’ e(r’)Qatt(Ir - r’l)) (10)
In some versions of NLDFT, where the functional derivative dpir’)/dp(r) in the Euler equation (eq 3) can be expressed in an explicit form as a function of the density p(r), the ILM method gives equations equivalent to those used p r e v i o ~ s l y . ’ ~ J A ~ J substantial ~ advantage while implementing the ILM method is achieved when dp(r’)l d p ( r ) cannot be expressed explicitly. In this case, the only method for the solution of the Euler equation (eq 31, discussed in the literature earlier, is a method of successive approximations developed by Tarazona8 a while ago, who had himself criticized it as a nonefficient and even sometimes diverging method. In the general case, several iterative schemes can be proposed for the solution of the system of equations (eqs 7, 8, and 2) by means of successive approximations. Depending on the type of the weighting function w,one or another scheme may be preferable. Therewith, it is possible to vary the terms between the left-hand side and right-hand side of eq 6, which offers additional flexibility in choosing the iterative scheme. The typical iterative scheme may be chosen as follows. Let po(r) be a zero approximation, for example, the solution of the LDFT equation (eq 101, or the solution of NLDFT obtained for a neighbor value of the chemical potential. Then, the zeroth order approximation for the smoothed density po(r)is calculated from eq 2 with p(r) = po(r), and the zeroth order approximation for the ILM function &(r) is calculated from eq 8 as AO = iL[po,p0l. Substitution of po(r), po(r),and Ao(r)into the right-hand side of eq 7 gives the first-order approximation for the density pl(r). At later stages of the iteration procedure this series of calculations is repeated until the required accuracy is attained. The weighted iterations8 can be easily employed for better convergence, if necessary. For example
+
+
Q, = alQL* (1- al)gL-l,2, = a*X,X
(1 - a2)XL-1, e, = a3,0,*+ (1where PI*, i L vand-plx , are intermediate quantities cal-
culated from p L - l , i L - l , pl-l as described above without employing weights a, (j= 1-31. The weights a, can be easily optimized in order to reduce the time of computations. Thus, ILM grant an additional degree of freedom in constructing the iterative scheme compared with the iterative method of solving the Euler equation (eq 2) proposed by Tarazona.8 The ILM method has been recently employed by Neimark et aLZ5and Ravikovitch et a1.26for calculating adsorption-desorption isotherms in nanopores in the wide range ofpore sizes from 1to 10 nm and studies ofhysteresis phenomena in nanomaterials. Acknowledgment. Useful discussions with A. Rabinovich, P. Ravikovitch, V. Gusev, and J. O’Brien are gratefully acknowledged. This work has been supported in part by NSF Grant CTS-9215604. LA9502305
Here, the ideal gas of density pg, corresponding to the given fluid chemical potentialp = kT In A3p,, is introduced as a reference state.
(25) Neimark, A. V.; Ravikovitch, P. I.; ODomhnaill, S. C.; Schuth, F.; Unger, K. K. Fundam. Adsorp. May 1995 submitted for publication. ( 2 6 )Ravikovitch, P. I.; O’Domhnaill, S. C.; Neimark, A. V.; Schuth, F.;Unger, K. K. Langmuir, submitted for publication.
