Restoring the consistency with the contact density theorem of a classical density functional theory of ions at a planar electrical double layer

Classical density functional theory (DFT) of fluids is a fast and efficient theory to compute the structure of the electrical double layer in the primitive model of ions where ions are modeled as charged, hard spheres in a background dielectric. While the hard-core repulsive component of this ion-ion interaction can be accurately computed using well-established DFTs, the electrostatic component is less accurate. Moreover, many electrostatic functionals fail to satisfy a basic theorem, the contact density theorem, that relates the bulk pressure, surface charge, and ion densities at their distances of closest approach for ions in equilibrium at a smooth, hard, planar wall. One popular electrostatic functional that fails to satisfy the contact density theorem is a perturbation approach developed by Kierlik and Rosinberg [Phys. Rev. A 44, 5025 (1991)] and Rosenfeld [J. Chem. Phys. 98, 8126 (1993)], where the full free-energy functional is Taylor-expanded around a bulk (homogeneous) reference fluid. Here, it is shown that this functional fails to satisfy the contact density theorem because it also fails to satisfy the known low-density limit. When the functional is corrected to satisfy this limit, a corrected bulk pressure is derived and it is shown that with this pressure both the contact density theorem and the Gibbs adsorption theorem are satisfied.

Figure 1

Verification that the corrected pressure ${\tilde{P}}^{\text{ES}}$ in Eq. (21) satisfies the contact density theorem in Eq. (6). On the $x$ axis are plotted the sums of the contact densities [the left-hand side of Eq. (6)] for 180 different cases (described below). The corresponding right-hand sides of Eq. (6) are plotted on the $y$ axis (symbols). Each data point is connected to the $x$ axis with a gray line to aid the eye in seeing clustered points. The line depicts when they are equal; if the symbols lie on the line, then the contact density theorem is satisfied. Because the wide range of different bulk densities produces different orders of magnitude in pressure, both axes use logorithmic scales. In all cases the system consists of one smooth, hard, planar wall with 0 or negative surface charge $\sigma$ $(log\left|\sigma \right|=-0.5,-1,-1.5$, and $-2$ with $\sigma$ in units of C/${\mathrm{m}}^2)$ and two ion species. The anion always has valence $-1$ and diameter $0.4$ nm. The cation has valences $+1,+2$, and $+3$ and diameters $0.2,0.4,0.6$, and $0.8$ nm. The bulk concentrations of the cation ${\rho }_{\text{c}}^{\text{bulk}}$ in Molar take values $log\left({\rho }_{\text{c}}^{\text{bulk}}\right)=0$ to $-5$ in steps of $-1$. Calculating the contact densities accurately enough to verify the contact density theorem (to 4 or more significant figures) is inherently difficult because of two major sources of error: the bulk boundary conditions must be placed far enough away from the wall and the grid points must be close enough together $(\Delta x$ is small enough) that the contact densities have sufficiently converged to their limiting values as $\Delta x\to 0$. (Note that fast Fourier transforms are generally used to solve the DFT equations, and so $\Delta x$ is uniform throughout the system.) The former condition requires a large number of points far from the wall when the density is low and the latter condition requires a large number of points near the wall when $\sigma$ is large. It was emperically found that ${10}^3$ nm of distance from the wall was sufficient for $log\left({\rho }_{\text{c}}^{\text{bulk}}\right)=-5$ and ${10}^{0.5}$ nm for $log\left({\rho }_{\text{c}}^{\text{bulk}}\right)=0$; for intermediate values of $log\left({\rho }_{\text{c}}^{\text{bulk}}\right)$, for every increase of 1 in $log\left({\rho }_{\text{c}}^{\text{bulk}}\right)$, the $log$ of the distance was increased by $0.5$ to scale with the change in the Debye length. For surface charges, when $log\left|\sigma \right|=-0.5,log\left(\Delta x\right)=-6$ was found to be sufficient $(\Delta x$ in nm) and $log\left(\Delta x\right)=-3$ for $log\left|\sigma \right|=-2$; for intermediate values of $log\left|\sigma \right|$, for every increase of $0.5$ in $log\left|\sigma \right|,log\left(\Delta x\right)$ was decreased by $-1$. For $\sigma =0,\Delta x={10}^{-3}$ nm was used. In the worst case, this produces $>{10}^9$ grid points, which is far too large for the 16 GB memories on the computers used. Therefore, all cases with $>2^{23}\approx 8\times {10}^6$ grid points (FFT requires the number of points to be a power of 2) were not considered. In all, 180 of the possible 288 cases were computed and are shown in the figure. For the ${\mu }_i^{\text{ref,ES}}$, the standard MSA was used [10, 11] and for $c_{i,j}^{\left(2\right),\text{ES}}(\left\{{\rho }_k^{\text{ref}}\right\};\mathbf{x},{\mathbf{x}}^{\prime })$, the approximate version given by Blum and Rosenfeld was used [20].