Formally `exact' first-order Taylor series expansion for density functional theory

The functional expansion of the non-uniform first-order direct correlation function (DCF) around the bulk density was truncated at the first order, then the functional counterpart of the Lagrangian theorem of differential calculus was employed to make the truncation formally exact. The actual procedure is as follows. According to the Lagrangian theorem and the definition of the DCF, the original expansion coefficient, i.e. the uniform second-order DCF, was replaced by the non-uniform second-order DCF whose argument is the appropriate mixture of the density distribution and the bulk density with an adjustable parameter determined by a hard-wall sum rule. With reference to an earlier paper (Khein A and Ashcroft N W 1999 Phys. Rev. E 59 1803), the non-uniform second-order DCF was then approximated by its uniform counterpart with a weighted density as its density argument. The truncated expansion was incorporated into the density functional theory formalism to predict the non-uniform hard-sphere fluid density distribution - in very good agreement with simulation data for three confining geometries: a single hard wall, a spherical cavity and a bulk hard-sphere particle whose resulting external potential leads to a radial distribution function of the bulk hard-sphere fluid whose prediction by the present theory was also in good agreement with the corresponding simulation data.


36.2
fundamental measure theory (FMT); and a recently proposed DFT based on the universality principle of the free-energy density functional and the test particle trick. A review of the above types of DFT can be found in [2,3] and the references therein. The present paper concentrates on the FPEA at the level of the first-order direct correlation function (DCF), but we did not strive to approximate the expansion coefficients of higher orders, as was done previously in the literature.
A formal 'Taylor series' expansion of the non-uniform first-order DCF C (1) (r; [ρ]) around the bulk density ρ b of a uniform system can always be written down: (1) Here each functional derivative, i.e. the expansion coefficient C (n) 0 (r, r 1 , . . . , r n−1 ; ρ b ), n ≥ 2, is evaluated at the initial ρ b . When such a series is terminated, it can be made an accurate representation by having the last functional derivative evaluated not at the initial ρ b but at some ρ b + λ(ρ − ρ b ), with λ between 0 and 1 [4]; this is actually a functional counterpart of the well-known Lagrangian theorem of differential calculus which states that equation (2), below, is exact if the value of λ is correctly chosen. According to the above procedure, we truncated the series at the first order; then equation (1) reduces to In the above equation, subscript 0 stands for the uniform case; its absence corresponds to the non-uniform case. It should be noted that the expansion coefficient C (2) 0 in equation (1) was replaced by C (2) in equation (2) because of the replacement of the bulk density ρ b by the non-uniform density field ρ b + λ(ρ − ρ b ); this replacement has to be done to be in agreement with the definition of the DCF. Equation (2) is exact and it does not include any approximation.
In the formalism of DFT, the density profile equation of a non-uniform single-component fluid reads where ϕ ext (r) is the external potential responsible for generation of the density distribution ρ(r), and β = 1/kT with k the Boltzmann constant and T the absolute temperature. Substituting equation (2) into (3) leads to The key is to approximate the non-uniform second-order DCF C (2) (r, r 1 ; [ρ b + λ(ρ − ρ b )]) and choose the appropriate mixing parameter λ. The former is approximated by extending the SWDA [5] for the non-uniform first-order DCF to its second-order counterpart; i.e., we made the approximation C (2) (r, r ; [ρ]) = C (2) 0 (|r − r |;ρ(r )).
From [6], one knows that one can approximate the non-uniform second-order DCF reasonably well through the uniform second-order DCF (see the bottom of the second page of [6]). So our extension of the SWDA for the non-uniform first-order DCF to its second-order 36.3 counterpart is reasonable. Considering the symmetry condition, we made r = (r + r )/2. Here,ρ is the weighted density and is as usually defined bỹ ρ(r) = dr ρ(r )w(|r − r |;ρ(r)). (6) In the spirit of the SWDA, which employed the normalized second-order DCF as a weighting function and made the bulk density ρ b the density argument of the weighting function, we can specify the weighting function w(|r − r |;ρ(r)) in equation (6) by After taking the SWDA for the DCF, equation (5) with equations (6) and (7) was introduced; equation (4) is now of the following form: whereρ The mixing parameter λ was specified by a hard-wall sum rule which specifies the bulk pressure P in terms of a hard-wall contact density ρ w : ρ w can be obtained from ρ(0.5σ) in equation (4 ) when the external potential has the following form: where σ stands for the fluid particle diameter, while P can be obtained from the equation of state for the corresponding bulk fluid. Although the hard-wall sum rule is valid only for a fluid in contact with a hard wall, the specified value of λ for this special case can also be used for other cases of the external potential, such as a spherical cavity and a soft wall, at the same bulk thermodynamic state. The reasoning is as follows: λ is included in expression (2), is the functional derivative of the excess Helmholtz free-energy density functional with respective to the density distribution, the excess Helmholtz free-energy density functional is independent of the external potential; thus C (1) (r; [ρ]) is also independent of the external potential. This means that a certain value of the associated λ suitable for some external potential (for example, a single hard wall) is also suitable for any other external potential (for example, a spherical cavity) at the same bulk thermodynamic state point [2,3], because equation (2) expressed C (1) (r; [ρ]) as a functional of the density distribution using the weighted density. To exemplify this point, we will calculate the density distribution profile for a hard-sphere fluid subjected to two different external fields from the value of λ determined for the same bulk hard-sphere fluid near a hard wall.
The numerical solution of equation (4 ) can be obtained by iteration; i.e., at the beginning, one assumes the density distribution ρ(r) in the region in which the fluid particle is confined to be ρ b , and the density distribution ρ(r) of other regions to be zero; then equation (4 ) combined with equation (8) can be used to produce a new ρ(r) in the left-hand side of equation (4 )  putting the assumed ρ(r) in the right-hand side of the same expression, equation (4 ). After the first calculation, the new density distribution is mixed with the old density distribution according to the following relation: where δ has to be a very small number, for example 0.1, or even 0.01 when the bulk density is high, to avoid divergence. Then one puts ρ(r) mix in the right-hand side of equation (4 ), and iterates until convergence is achieved. We test the performance of the present method with a hard-sphere fluid near a single hard wall and confined in a spherical cavity. For the former, the external potential is of the same form as in equation (10); for the latter, the following form: The function needed, i.e. the uniform second-order DCF, is obtained from the Percus-Yevick approximation [7], while the equation of state for the bulk hard-sphere fluid is the Carnahan-Starling empirical one [8]. The calculated density distribution is plotted against the corresponding simulation data [9,10] in figures 1 and 2 for the two cases of external potentials. The agreement with the simulation data is very good.
To further test the universality of λ-i.e. whether the value of λ determined for the bulk hard-sphere fluid near a single hard wall for certain state point can be used for the calculation for any other cases of external potentials (e.g. different geometry and/or different external potential parameters) for the same state point-we applied the present approach to calculate the radial distribution function of the bulk hard-sphere fluid.  Figure 2. Density profiles of a hard-sphere fluid confined in a spherical cavity with a hard-wall R = 4.5σ for two cases of bulk densities. The lines correspond to the predictions of the present method; the symbols stand for the corresponding computer simulation data [10]. The ordinate for the upper curve should be shifted down by one unit.
Following the test particle trick of Percus [11], the radial distribution function g(r) of the bulk fluid can be related to the non-uniform density distribution ρ(r) by the following relation: where the external potential responsible for the generation of ρ(r) results from a so-called test particle chosen from the bulk and situated at the origin. The interaction potential between the test particle and a any other particle in the bulk is exactly the external potential, i.e.
In figure 3, the calculated radial distribution functions derived from the present DFT approach are displayed together with the corresponding simulation data [12] for the hard-sphere fluid at several bulk densities. It was shown that the agreement with the simulation data is also very good, even for very high bulk density such as ρ b σ 3 = 0.8. Since the external potential generated from a spherical particle is completely different from that of a single hard wall, the good accuracy shown in figure 3 indicates the universality of the mixing parameter λ, i.e. its independence of the external potential.
In the case of an impenetrable solid, the one-particle cavity correlation function t(z) satisfies the following relation: where µ is configurational chemical potential, while the present approach gives  Figure 3. Comparison of the radial distribution function of the bulk hard-sphere fluid from the present DFT with the Monte Carlo simulation data [12] for several bulk densities. The ordinate for the curve for ρ b = 0.8 should be shifted down by two units, while that for ρ b = 0.7 should be shifted down by one unit.
Although the mathematical forms of equations (15) and (16) differ greatly from each other, one can evaluate them numerically. By employing the analytical expressions for βµ and C (2) 0 from the PY approximation and the specified value of λ as shown in figures 1-3, which is derived from the combination of the hard-wall sum rule and the CS equation of state, we found that the maximal calculated percentage relative error (PRE) between βµ and − dr ρ b C (2) is 17.337% for the highest calculated bulk density ρ b = 0.813 (λ = 0.47), while the minimum calculated PRE is 11.282% for the lowest calculated bulk density ρ b = 0.575 (λ = 0.51). The largeness of the PRE given by the present approach arises from two sources. First, to calculate the above PRE, we used βµ and C (2) 0 from the PY approximation, while the value of λ employed is determined from the CS equation of state; this is not self-consistent. When we replaced the CS equation of state by the PY compressibility equation of state to determine the value of λ, we found that the PRE was reduced greatly. For example: for state point ρ b = 0.575 the PRE is now only 7.77% (λ = 0.482), for state point ρ b = 0.75 the PRE is now 10.28% (λ = 0.45), while the previous PRE is 15.42% (λ = 0.48). The other source of largeness of the PRE is obviously the approximation for the non-uniform second-order DCF; this is the only approximation included in the present DFT. It is well known that not all of the existing DFTs satisfy the exact relation (15). In fact, whether the exact equation (15) is satisfied or not is not important for the DFT calculation for density profiles, as shown by the present calculation. Although we employed the CS equation of state for the determination of λ and this value of λ resulting from the CS equation of state leads to a larger PRE than the value of λ resulting from the PY compressibility equation of state, the former predictions about the density profile are far more accurate than the latter, at least near the hard-wall contact, because the PY compressibility equation of state predicts higher bulk pressure than the CS equation of state; thus it has to predict a higher single-hard-wall contact density than 36.7 the CS equation of state, while the single-hard-wall contact density prediction resulting from the CS equation of state is very accurate as we see from figure 1. The other example is the previously proposed SWDA [5], which satisfies equation (15) exactly, but its predictions for the density profile are not good. From the above calculation for λ, one finds that the value of λ resulting from the PY compressibility equation of state is lower than the value of λ resulting from the CS equation of state for the same bulk state point. Obviously a lower value of λ leads to prediction of a higher hard-wall (not necessarily single or parallel) contact density for all confining geometries. This conclusion can be arrived at from the fact that the usual first-order perturbation DFT always predicts higher hard contact density for various confining geometries and the present approach reduces to the usual first-order perturbation DFT when the value of λ is reduced to zero. From figure 2 one can see that even the value of λ resulting from the CS equation of state led to a prediction of a slightly higher hard contact density than the simulation data for the bulk state, ρ b = 0.75. It can be envisaged that the value of λ resulting from the PY compressibility equation of state will lead to far too high a predicted hard contact density for the same state point.
To summarize, the present procedure applies the functional counterpart of the Lagrangian theorem of the differential calculus to DFT; this is the first time that this has been done since the introduction of the DFT in classical statistical thermodynamics several decades ago. In fact, perturbation expansion is a basic and general tool in the field of theoretical physics; from the viewpoint of mathematics the functional counterpart of the Lagrangian theorem of differential calculus can be applied to any theoretical problems concerned with the functional Taylor series expansion. So the present procedure introduces a general methodology to the field of theoretical physics; its range of application will be very wide.
The present method can be applied to non-uniform fluids of other interaction potentials by two routes. First, regarding the non-hard-sphere interaction potential as a whole, the C (2) 0 (r; ρ b ) needed is obtained by solving the OZ equation numerically. This route has been pursued numerically for the non-uniform Lennard-Jones fluid based on the WDA [13]. Secondly, one divides the interaction potential into a hard-sphere-like part and a perturbation part, the former treated using the previously available accurate and simple DFT for non-uniform hard-sphere fluids and the latter treated by the present method. Pursuing the present idea at the level of the free energy will open a new route for development of the DFT. These studies will be reported in separate papers.