Benchmarking Orbital-Free Density-Potential Functional Theory of Electrified Metal-Solution Interfaces

TL;DR

This paper benchmarks a new computational model for studying electrified metal-solution interfaces using orbital-free density-functional theory.

Contribution

The study introduces a benchmark framework for evaluating DPFT models at interfaces and identifies optimal parameters for the TFvW functional.

Findings

The TFvW kinetic energy functional outperforms the Pauli-Gaussian functional in modeling the electrical double layer.

A smaller gradient coefficient in the TFvW functional improves accuracy for EDL modeling.

The benchmark framework enables better assessment of orbital-free DFT models for electrochemical interfaces.

Abstract

The electrical double layer (EDL) at the metal–solution
interface is a nanoscale region where quantum mechanical metal electrons
meet almost classical electrolyte species. Describing metal electrons
with orbital-free density-functional theory (DFT), the recently developed
density-potential functional theoretical (DPFT) model constitutes
a computationally efficient approach to modeling the EDL. However,
the performance of orbital-free DFT is less studied for interfaces
than for bulk phases. Herein, we develop a constant-potential Kohn–Sham–Poisson–Boltzmann
theory with exact kinetic energy as a benchmark for DPFT models. Solving
Kohn–Sham and Poisson–Boltzmann equations self-consistently,
we obtain electron density, electrostatic potential, and double-layer
capacitance of the EDL, which are then used to assess DPFT models
with Thomas–Fermi–von Weizsäcker (TFvW) or Pauli-Gaussian
kinetic energy functional. In general, TFvW outperforms the Pauli-Gaussian
kinetic energy functional for modeling EDL. In addition, a much smaller
gradient coefficient in the TFvW functional than the default value
of 5/3 is suggested for modeling the EDL. These findings are instrumental
to the future development of orbital-free DFT for electrochemical
interfaces.

Figures

Schematic diagram of 1D planar symmetrical geometry consisting of a metallic slab in the middle of two electrolyte solutions. L M and L S are the dimensionless thicknesses of the metal and electrolyte solution, respectively.

Flowchart for self-consistently solving the coupled Kohn–Sham–Poisson–Boltzmann equations.

Basic results of the Kohn–Sham–Poisson–Boltzmann model at μ̃ e = −3.6 eV. (a) Disbributions of the first five energy levels ε i of electrons, (b) distributions of the first five normalized wave functions φ̅ i , (c) distribution of electrostatic potential ϕ, (d) distributions of electron density and subband’s contributions, (e) distributions of dimensionless cation n̅ c and anion density n̅ a, (f) comparison of residuals among the fixed-point iteration, linear mixing method, and predictor–corrector method. Parameters used in calcualtion are as follows: n̅ cc = 0.01, n c,a bulk = 1 M (number density is concentration times Avogadro constant throughout this paper), L M = 20, L S = 30, d c,a = 4a 0, D l = 0.5/6 eV, and βc,a = 1.

Comparison between the DPFT model with the TFvW functional and Kohn–Sham–Poisson–Boltzmann model at μ̃ e = −3.57 eV in terms of (a) normalized electron density error Δn e,n = (n̅ e DPFT – n̅ e Kohn–Sham–Poisson–Boltzmann)/n̅ cc, (b) normalized electron density n e,n = n̅ e/n̅ cc, (c) electrostatic potential, (d) electrostatic potential at the region of the metal/solution interface. Parameters used in calculation are the same as those in Figure .

Comparison between the DPFT model with the Pauli-Gaussian functional and Kohn–Sham–Poisson–Boltzmann model at μ̃ e = −3.57 eV in terms of (a) normalized electron density error Δn e,n = (n̅ e DPFT – n̅ e Kohn–Sham–Poisson–Boltzmann)/n̅ cc, (b) normalized electron density n e,n = n̅ e/n̅ cc, (c) electrostatic potential, (d) electrostatic potential at the region of the metal/solution interface. Parameters used in calculation are the same as those in Figure .

Comparison between DPFT model using (a) TFvW functional, (b) Pauli–Gaussian functional and Kohn–Sham–Poisson–Boltzmann model in terms of double-layer capacitance at n̅ cc = 0.01. Parameters used in the calculation are the same as those in Figure .

Comparison between DPFT model with (a) TFvW functional and (b) Pauli–Gaussian functional and Kohn–Sham–Poisson–Boltzmann model results in terms of double-layer capacitance at n̅ cc = 0.005. Other parameters used in the calculation are the same as those in Figure .

Comparison between DPFT model with TFvW functional and Pauli–Gaussian functional and Kohn–Sham–Poisson–Boltzmann model, in terms of double-layer capacitance (a, b), μ̃e PZC (c, d), and mean error (e, f), respectively. n̅ cc = 0.01, other parameters are the same as those in Figure . The exchange-correlation potential is described using the PBE functional.

Comparison between the DPFT model with the TFvW functional and Kohn–Sham–Poisson–Boltzmann model at μ̃e = −2.45 eV in terms of (a) electrostatic potential ϕ, (b) normalized electron density n e,n, (c) normalized cation density n c,n, (d) normalized anion density n a,n, (e) normalized solvent density n s,n, (f) relative dielectric permittivity ϵ/ϵ0, (g) surface free charge σfree, (h) double-layer capacitance C dl, and (i) potential of zero charge μ̃e pzc. Parameters used in the calculation are as follows: n̅ cc = 0.01, r s = 2.75 Å, r a = 2.7 Å, r c = 3.6 Å, n c,a bulk = 0.1 M, n s bulk = 55.6 M, p s = 1.8 D, ϵop=ϵopM+(ϵopS−ϵopM)2⁡erfc(−βop(|X|−LM)) , ϵop M = 3ϵ0, ϵop S = 6ϵ0, D s = 0.5 eV, d s = 2a 0, βs = 1, L M = 20, L S = 120, other parameters are the same as those in Figure . The exchange-correlation potential is described using the PBE functional.

Tables

1. Model Parameters for the Base Case

symbol	item	values
Structural parameters
L M	Dimensionless metal thickness	20
L S	Dimensionless electrolyte solution thickness	30
Material properties
n cc	Cationic core density	0.01a 0 –3
n c,a bulk	Bulk ions density	1 M
a	Ion size	0 (neglecting ions size)
ϵop	Dielectric permittivity	ϵ0
Empirical parameters
d c,a	Equilibrium distance of cation/anion to the metal surface	4a 0
D l	Depth of the Morse potential well	0.5/6 eV
βc,a	Coefficient controlling the width of the Morse potential well	1
θ T	Empirical parameter in TFvW energy functional	Calibrated to Kohn–Sham–Poisson–Boltzmann models
μPG	Empirical parameter in Pauli–Gaussian energy functional

Full text

Introduction

In response to the potential difference applied at the interface of a metal with an electrolyte solution, metal electrons, electrolyte ions and solvent molecules arrange into an electrical double layer (EDL). Understanding the microstructure and dynamics in the EDL is key to modulating the capacitive response and the charge transfer kinetics at the metal-solution interface. Modeling and computation are now an integral part of efforts that strive to unravel the properties of the nanoscale EDL, complementing experimental investigations.

Modeling of the EDL has a long history and comes in different flavors, see recent reviews. ?−? ? The classical Gouy–Chapman–Stern model conceptualizes the EDL as a serial connection of a rigid Helmholtz plane and a diffuse Gouy–Chapman layer. ?−? ? Historically, improvements to the Gouy–Chapman–Stern model have been made along three directions. First, the original point-charge description of electrolyte ions has been improved by considering steric ion effects, ?−? ? solvent polarization, ?,? and nonlocal short-range correlations. ?−? ? Second, the Helmholtz plane has been divided into an inner Helmholtz plane for specific adsorption and chemisorption ?,? and an outer Helmholtz plane for nonspecifically adsorbed ions. The first water layer at the IHP exerts a great influence on the double-layer capacitance.? Third, the ideal-metal boundary in the original Gouy–Chapman–Stern model neglects the electron spillover and metal electronic effects. A jellium description of the metal has been devised to address this shortcoming. ?,?,?,? The metal is treated as an inhomogeneous electron gas, which is situated within a positive background charge contributed by the cationic cores of the metal. There exists a gap between the metal and the compartment of the electrolyte solution, which may shrink or expand as the electrode potential is changed. ?,?

In the last three decades, the rapid growth in computing power has enabled fast-paced progress in density functional theory (DFT)-based models of the EDL. ?−? ? ? ? In a typical setup, these models represent the metal as a thin slab composed of a few layers of metal atoms, the electrolyte solution either as an implicit dielectric medium or explicitly as a thin water film with a few ions in ab initio molecular dynamics simulations or using a hybrid implicit/explicit scheme. These DFT-based models have greatly enhanced the atomistic understanding of the EDL, especially in the roles of interfacial water molecules in determining the potential of zero charge, ?,?,? the Helmholtz capacitance, ?−? ? and negative capacitance response induced by chemisorption of water. ?,? However, the high computation cost involved in solving the Kohn–Sham equations has limited DFT-based models to a periodic structure with a few hundreds of atoms. Therefore, the statistical sampling of the electrolyte solution is usually poor, resulting in an unsatisfactory description of the local reaction environment under constant potential conditions.

Despite the progress in classical Gouy–Chapman–Stern-type models and DFT-based models, modeling the EDL still faces two key challenges. The first one relates to treating electronic, ionic and solvent effects on an equal footing under constant potential conditions. The second one is about nonperiodic, mesoscopic EDLs met in many applications render first-principles calculation difficult. These challenges have driven us to develop a density-potential functional theoretical (DPFT) method.

The DPFT method combines orbital-free DFT for metal electrons and statistical field theory for the adjacent electrolyte solution. The grand potential of the EDL is then expressed as a hybrid functional of the electron density (n e) and the electrostatic potential (ϕ). Resulted from a coarse-grained description for the electrolyte solution, ϕ becomes a primary variable with an equal status as n e, contrasting the Kohn–Sham DFT where they are dual variables. This is why we refer to it as a DPFT approach. The fundamental idea behind the DPFT method is that potential should be treated as a primary variable in modeling electrochemical interfaces, contrasting current DFT models where potential effects are often added afterward as a correction term. The gap between the metallic and electrolytic compartments, which has been an enigma issue in previous jellium models, ?,? is naturally formed out of short-range interactions between metal electrons and solution particles described using empirical relations such as Morse potentials or Lennard-Jones potentials. The DPFT approach has contributed to understanding the two peaks of double-layer capacitance curves at silver electrodes,? the solvent dependence of the potential of zero charge of Au (111)–organic solution interfaces,? the short-range disjoining pressure between two metal surfaces,? etc. The DPFT approach has demonstrated its effectiveness in describing the electro-ionic perturbations in the overlapping EDL at three-dimensional supported nanoparticles.?

A main shortcoming of the DPFT approach lies in the orbital-free approximation of the kinetic energy of metal electrons. Though orbital-free DFT has been carefully evaluated for atoms,? molecules? and bulk phases,? its performance of EDLs that involve a large density gradient at the metal-solution interface is poorly known. Our goal here is to develop a first-principles based benchmarking approach for DPFT. To focus on this specific goal, we will consider an ideal metal/solution interface that is uniform in two dimensions, and exhibits heterogeneity only in the perpendicular direction.? Therefore, the three-dimensional problem can be reduced to a one-dimensional problem. Metal electrons will be described using Kohn–Sham equations with exact kinetic energy for noninteracting electrons. Focusing on metal electronic effects, we will use a primitive model for the electrolyte solution, which is modeled as point charges dispersed in a dielectric continuum with a constant dielectric permittivity.? Taken together, the above descriptions result in a Kohn–Sham–Poisson–Boltzmann theory of metal/solution interfaces. The Kohn–Sham–Poisson–Boltzmann equations are solved using an in-house self-consistent numerical scheme, serving as benchmark for calibrating orbital-free kinetic energy functionals in DPFT models.

The rest of the paper is organized as follows. We first review previous efforts on similar problems and identify the specific technical gaps that this work contributes to filling out. Second, we introduce the setup of Kohn–Sham–Poisson–Boltzmann and DPFT models, including the geometric structure, the governing equations, and the boundary conditions. Third, we propose a self-consistent numerical scheme that solves the coupled Kohn–Sham–Poisson–Boltzmann equations. Finally, we compare Kohn–Sham–Poisson–Boltzmann and DPFT models.

Kohn–Sham–Poisson Models

Kohn–Sham–Poisson models have been developed for metal–vacuum systems ?−? ? and semiconductors ?,? to calculate the electronic structure. Regarded as an eigenvalue problem, the Kohn–Sham equation is solved using the eigenvalue method with the plane wave basis ?,? or atomic orbital basis, ?,? or the finite difference method.? With the assumption of homogeneous electron gas in the bulk metal, the Kohn–Sham equation can be formulated as an ordinary differential equation. The Poisson equation is a second-order partial differential equation, which can be solved using numerical methods such as finite difference method for regular geometry,? finite element method? and finite volume method? for irregular geometry. Solving coupled Kohn–Sham–Poisson equations is challenging. Charge distribution obtained from the Kohn–Sham equation may not conform to boundary conditions of the Poisson equation for the electrostatic potential. In the following, we review the progress in solving the Kohn–Sham-Poisson equations by formulating the Kohn–Sham equation as an ordinary differential equation and as an eigenvalue problem, respectively.

Lang and Kohn treated the Kohn–Sham equation as an ordinary differential equation to study charge densities, surface energies, and work functions of metal–vacuum surfaces. ?,? In their approach, the electron density was updated iteratively using a perturbational Newton–Raphson-type functional, which is more robust than simple fixed-point iteration. The electrostatic potential was calculated from direct integration using the Columbic potential expression, which is prone to numerical instability. Monnier and Perdew generalized the method of Lang and Kohn improving solution of electrostatic potential with an integral equation and further by considering high-order discrete-lattice perturbation potential δv to study face-dependent surface energies, density profiles and work functions of nine simple metals.? Their results can be reduced to those of Lang and Kohn in the limit of weak δv. Subsequently, Yan? and Della Sala? employed the numerical routine of Lang and Kohn to solve the Kohn–Sham-Poisson equations and used it to calibrate the parameters of the orbital-free DFT. Posvyanskii and Shul’man improved the method of Lang and Kohn by introducing a predictor-corrector iteration scheme to update electron density.? They further extended this scheme to a quasi-two-dimensional electron gas in the accumulation layer at the surface of a semiconductor with a degenerate electron gas in the bulk phase.?

The idea of treating the Kohn–Sham equation as an eigenvalue problem can be traced back to the Kohn–Sham DFT.? To improve convergence of solving the Kohn–Sham–Poisson equations, many iteration schemes have been proposed. Stern systematically compared three different iteration methods of the electron density: including the fixed-convergence-factor method, the extrapolated-convergence-factor method, and the perturbational-iteration method to study semiconductor inversion layers.? The last method is similar to the perturbational approach of Lang and Kohn. ?,? The fixed-convergence-factor method is now also called the linear mixing or under-relaxation method. Hamann introduced the Anderson acceleration? into electronic structure calculation and studied semiconductor charge densities with hard-core and soft-core pseudopotentials.? Around the same time, Pulay independently proposed an acceleration method of the direct inversion in the iterative subspace and used it for molecular calculations.? Nowadays, the Anderson acceleration and Pulay’s method involve very similar algorithms. Dederichs and Zeller conducted a formal study of the linear mixing method in electronic structure calculations and compared its convergence with the Anderson acceleration for 3d impurities in Cu.? Using the predictor-corrector iteration method, Luscombe et al. studied the electron confinement in the GaAs quantum wire by solving the coupled Kohn–Sham–Poisson equations, and compared the results with those obtained from the Thomas–Fermi theory.? After that, Trellakis et al. systematically compared the convergence of the linear mixing method and the predictor–corrector method of solving the two-dimensional Kohn–Sham–Poisson equations in a GaAs-AlGaAs-based quantum wire.? They found that the predictor-corrector method converged much faster than the linear mixing method and the exchange-correlation potential impaired both convergence and stability. Recently, Gil-Corrales et al. studied the electronic properties of GaAs quantum well wires with various cross-sectional shapes using the predictor-corrector iteration method.? A more comprehensive summary and discussion on iteration methods of solving the coupled Kohn–Sham–Poisson equations can be found in recent review article by Novák et al.?

Although many studies about self-consistent solutions of the Kohn–Sham–Poisson equations exist, they have been focusing on the metal-vacuum interface or the semiconductor system under electroneutrality. The self-consistent solution of Kohn–Sham–Poisson–Boltzmann equations for metal-solution interfaces under constant potential conditions has not been reported to the best of our knowledge. In this work, we propose a method to solve Kohn–Sham–Poisson–Boltzmann equations under constant potential conditions, and use it to assess the accuracy of DPFT where the kinetic energy of electron densities is described using the Thomas–Fermi–von Weizsäcker ?−? ? and Pauli–Gaussian kinetic energy functionals.?

Model Setup

In this section, we specify the computational domain tailored to focus on metal electronic effects. Thereafter, we introduce the governing equations and their boundary conditions for both Kohn–Sham–Poisson–Boltzmann and DPFT models. Subsequently, we make a nondimensionalization for governing equations. In the following, we specify how to achieve constant-potential conditions in Kohn–Sham–Poisson–Boltzmann and DPFT models, and develop a numerical flowchart for self-consistently solving the Kohn–Sham–Poisson–Boltzmann model. Finally, we parametrize the models and calculate and analyze the solutions obtained.

1D Symmetrical Geometry

We consider a 1D symmetrical geometry consisting of a metallic slab in the middle of two electrolyte solutions, as shown in Figure. The symmetrical geometry simplifies the boundary conditions for electronic wave functions. The treatment of boundary conditions of wave functions for an asymmetrical geometry was discussed in earlier studies. ?,? Wave functions are assumed to be uniform in the other two dimensions, where translational symmetry is preserved. This symmetry allows the wave functions to be represented as plane waves parallel to the surface, while the variation along the x direction (normal to the interface) is treated explicitly. Therefore, the original 3D problem is reduced to a 1D problem exhibiting variations only in the coordinate normal to the metal surface (the detailed formulation is provided in supplementary note 1 in the Supporting Information).

Schematic diagram of 1D planar symmetrical geometry consisting of a metallic slab in the middle of two electrolyte solutions. L M and L S are the dimensionless thicknesses of the metal and electrolyte solution, respectively.

Governing Equations

Kohn–Sham–Poisson–Boltzmann Models

The governing equations for the Kohn–Sham–Poisson–Boltzmann model consist of the Kohn–Sham equation and Poisson–Boltzmann equation.

The 1D Kohn–Sham equation for a system of fictitious noninteracting electrons, subjected to an effective potential v eff can be expressed as?

$$[eqn]$$

where ℏ is the reduced Planck constant, m e is the mass of an electron, ε_ i _ is the energy, φ_ i _ is the wave function, and v eff is the effective external potential energy given by

$$[eqn]$$

where the first term is the electrostatic potential energy with ϕ­(x) being the local electrostatic potential and e 0 the elementary charge of an electron. Herein, −e _0_ϕ­(x) already includes the contributions of the Hartree potential energy and external potential energy, namely, −e _0_ϕ­(x) = v H + v ex. v XC(x) is the exchange-correlation (XC) potential energy, which is a function of local electron density n e(x) and its gradient ∇n e(x). For the base case, we use the primitive v XC(x), namely, the Dirac–Wigner functional? with the local density approximation:

$$[eqn]$$

where $[eqn]$ is the atomic energy with Bohr radius a 0 and vacuum permittivity ϵ_0_. The first term in the bracket, with $[eqn]$ and electron density n e, represents the exchange potential of a uniform electron gas; the second term, with c 4 = – 0.056, c 5 = 0.079, is Wigner’s expression for the correlation potential.?

The 1D Poisson–Boltzmann equation obtained from the DPFT functional can be expressed as (the detailed derivation is provided in supplementary note 2 in the Supporting Information):?

$$[eqn]$$

where ϵ_op_ is the dielectric permittivity, n cc is the metal catonic core density, n c/n a is the cation/anion density in the solution, and q c/q a is the absolute value of cation/anion valence.

n c/n a is described using the modified Boltzmann equation that incorporates ion size effect via a lattice-gas description: ?,?

$$[eqn]$$

where n c,a ^bulk^ is the cation/anion density in the solution bulk, k B is the Boltzmann constant, T is the temperature, χ_c,a_ = n c,a ^bulk^ a ^3^ is the bulk volume fraction of ions, with a being the ion size, and w c,a is the Morse potential that describes short-range interactions between metal and electrolyte solution. We use the repulsive part of w c,a to prevent solution cations/anions from penetrating into the metal, expressed as:?

$$[eqn]$$

where D l is the strength of the Morse potential, β_c,a_ is the coefficient controlling the decay of the Morse potential, d and d c,a denote the distances from the position and equilibrium position of cation/anion to the metal surface, respectively.

The Kohn–Sham and Poisson–Boltzmann equations are coupled via ϕ and n e. The Kohn–Sham equation solves for φ_ i _ and ε_ i _ rather than n e, therefore, we need to calculate n e from φ_ i _ and ε_ i _, ?,?

$$[eqn]$$

where μ̃_e_ is the electrochemical potential of electrons, equivalent to the Fermi level. μ̃_e_ is usually defined as μ̃_e_ = μ_e_ – e 0_ϕ, where μ_e is the chemical potential composed of the kinetic energy and exchange-correlation potential energy, namely, $[eqn]$ . The detailed derivation of n e is provided in supplementary note 5 in the Supporting Information.

DPFT Models

The governing equations for the DPFT model consist of an equation describing the electron density and the Poisson–Boltzmann equation. The latter is the same as that in the Kohn–Sham–Poisson–Boltzmann model.

For the DPFT model using the Thomas–Fermi–von Weizsäcker (TFvW) kinetic energy functional, a variational analysis of the grand potential g with respect to n e using the Euler–Lagrange equation gives:

$$[eqn]$$

where g is the volumetric grand potential of the electric double layer shown in eq (S16) in the Supporting Information.

After algebraic manipulations provided in supplementary note 3 in the Supporting Information, we obtain the controlling equation of the electron density:?

$$[eqn]$$

where θ_T_ is the coefficient of gradient terms in kinetic energy, $[eqn]$ is the Thomas–Fermi kinetic energy functional, μ̃ e is the electrochemical potential of electrons, ∇̅ = a 0∇ is the dimensionless Laplace operator, and n̅ e = n e a 0 ^3^ is the dimensionless electron density.

For the DPFT model using the Pauli–Gaussian kinetic energy functional, we obtain a similar equation for the electron density:

$$[eqn]$$

where $[eqn]$ is the reduced gradient term, and μ_PG_ is a coefficient in the Pauli–Gaussian functional. A detailed derivation of eq is provided in the supplementary note 4 in the Supporting Information. The Pauli–Gaussian functional is reduced to the TFvW kinetic energy functional when μ_PG_ = 0, and eq is reduced to eq.

Boundary Conditions

Boundary conditions to close the Kohn–Sham and Poisson–Boltzmann equations are shown in Figure,

$$[eqn]$$ $$[eqn]$$

which means that the wave function cannot propagate into the bulk solution and the electrostatic potential in the bulk solution is taken as a reference, see more discussion in supplementary note 6 in the Supporting Information.

For the governing equation of the electron density in the DPFT, the boundary condition is as follows,

$$[eqn]$$

which means that metal electrons are absent in the bulk solution.

Nondimensionalization

To facilitate numerical solutions, we define dimensionless variables, marked with an overbar, as follows:

$$[eqn]$$

We obtain the following dimensionless Kohn–Sham and Poisson–Boltzmann equations:

$$[eqn]$$ $$[eqn]$$

The dimensionless boundary conditions in the bulk solution, x̅ = ±(L S + L M), naturally become φ = 0 and ϕ̅ = 0.

The dimensionless electron density is written as:

$$[eqn]$$

where $[eqn]$ is the normalized wave function.?

Constant-Potential Conditions

Here we explain how constant potential conditions are preserved in DPFT and Kohn–Sham–Poisson–Boltzmann models. For DPFT models, controlling the constant potential is realized by controlling μ̃ e, namely, the electrochemical potential of electrons, in governing equations in eqs and ?. For Kohn–Sham–Poisson–Boltzmann models, adjusting μ̃ e in eq can also achieve constant-potential calculations.

Numerical Flowchart of Self-Consistent Solution

In this part, we first introduce the numerical method of solving the Kohn–Sham equation and Poisson–Boltzmann equation, separately, then present the flowchart of self-consistent solutions.

Finite Difference Method for Solving the Kohn–Sham Equation

Following Halpern et al.,? we use the finite difference method with a second-order center difference scheme to discretize the Kohn–Sham equation. For a 1D uniform mesh grid, the discretization of eq leads to:

$$[eqn]$$

where the index i denotes the grid point and d x̅ is the uniform step width on the 1D mesh. ** H ** _ ** ij ** _ is the Hamiltonian matrix:

$$[eqn]$$

We notice ** H ** _ ** ij ** _ is a symmetric tridiagonal matrix with the uniform mesh, which can be easily solved with numerical method. While the nonuniform mesh step is preferable for certain problems, this would destroy the symmetry of the matrix ** H ** _ ** ij ** _ and increase numerical complexities. Tan et al. used a proper matrix transformation that preserves the symmetry of the discretized Kohn–Sham equation even with a nonuniform mesh step, thereby reducing the computation time.?

Finite Element Method for Solving the Poisson–Boltzmann

Equation

The Poisson–Boltzmann equation shown in eq can be solved using the finite difference method, finite element method, and finite volume method. Considering the stability and capability of extension to more complex cases, we use the finite element method to solve it, with details provided in supplementary note 6 in the Supporting Information.

Since the Poisson–Boltzmann equation is nonlinear, an iteration strategy is needed. We define a residual vector R(ϕ̅ _ i _ ^ k ^), which is expected to be zero,

$$[eqn]$$

where the superscript k denotes the k ^th^ iteration, K _ ij _(ϕ̅ _ i _ ^ k ^) is the stiffness matrix, F(ϕ̅ _ i _ ^ k ^) is the source vector, and c(ϕ̅ _ i _ ^ k ^) are coefficients to be determined.

Using the first-order Taylor expansion, we obtain

$$[eqn]$$

with δϕ̅_ i _ ^ k ^ being the changing value between k ^th+1^ and k ^th^ iteration.

Therefore, we obtain

$$[eqn]$$

Equation is rewritten as

$$[eqn]$$

where $[eqn]$ is the so-called Jacobian matrix.?

In the k ^th^+1 iteration, the updated ϕ̅_ i _ ^ k+1^ = ϕ̅_ i _ ^ k ^ + δϕ̅_ i _ ^ k ^, which is referred to as the Newton–Raphson iteration.? Repeating iteration until ϕ̅_ i _ ^ k ^ is unchanged, we can obtain convergent solutions of ϕ̅_ i _.

Notably, the boundary conditions will modify the stiffness matrix and source term, see more details in supplementary note 7 in the Supporting Information.

Flowchart for Self-Consistently Solving the Kohn–Sham–Poisson–Boltzmann

Equations

A self-consistent method is developed to solve the coupled Kohn–Sham–Poisson–Boltzmann equations,? as shown in Figure. Initial guesses of ϕ̅_0_ and n̅ e ^0^ used are detailed in supplementary note 8 in the Supporting Information. Using the finite difference method with a dimensionless grid step of 0.04 to solve the Kohn–Sham equation, we obtain the eigenvalues and eigenvectors, namely, the distributions of φ_ i _ and ε̅, respectively. Thereafter, we calculate an updated electron density n̅ e ^ m ^ using eq, with the superscript m denoting the m ^th^ iteration. Subsequently, the updated electron density is used to solve the (modified) Poisson–Boltzmann equation to obtain a new electrostatic potential ϕ̅_ m . This step is carried out using the finite element method combined with the Newton–Raphson iteration, which ensures high numerical accuracy and stability when dealing with nonlinear terms. The newly computed potential ϕ̅ m+1_ and electron density n̅ e ^ m+1^ are compared to their previous values to evaluate the convergence tolerance, defined as |ϕ̅_ m+1_ – ϕ̅_ m |&|n̅ e ^ m+1^ – n̅ e ^ m ^|, which is required to be smaller than 10^–6^. If both quantities are smaller than this threshold, the self-consistent field cycle is considered to be converged, and the results are reported. Otherwise, we repeat the iteration with updated ϕ̅ m _ and n̅ e ^ m ^ until the errors are smaller than the predefined values. To evaluate the influence of mixing parameters, Figure(f) further compares convergence performance obtained with linear mixing parameters 0.05 and 0.5. For the DPFT models described using the TFvW and Pauli–Gaussian functionals, a denser mesh with a step of 0.02 at the metal–solution interface is required to accurately capture the large gradients.

Flowchart for self-consistently solving the coupled Kohn–Sham–Poisson–Boltzmann equations.

*Basic results of the Kohn–Sham–Poisson–Boltzmann model at μ̃ e = −3.6 eV. (a) Disbributions of the first five energy levels ε i of electrons, (b) distributions of the first five normalized wave functions φ̅ i , (c) distribution of electrostatic potential ϕ, (d) distributions of electron density and subband’s contributions, (e) distributions of dimensionless cation n̅ c and anion density n̅ a, (f) comparison of residuals among the fixed-point iteration, linear mixing method, and predictor–corrector method. Parameters used in calcualtion are as follows: n̅ cc = 0.01, n c,a bulk = 1 M (number density is concentration times Avogadro constant throughout this paper), L M = 20, L S = 30, d c,a = 4a 0, D

l = 0.5/6 eV, and βc,a = 1.*

Model Parameters

We divide the parameters into three groups: structural parameters including L M and L S, materials properties including n cc, n c,a ^bulk^, a and ϵ_op_, and empirical parameters d c,a, D _ l , β_c,a, θ_T_ and μ_PG_. Table gives the model parameters for the base case.

1: Model Parameters for the Base Case

For the structural parameters, L M and L S denote dimensionless thicknesses of the metal and electrolyte solution ranging from 1 nm to tens of nanometers in typical DFT calculations.? Here we set L M = 20 and L S = 30 referenced to a 0 (0.529 Å).

For materials properties, n cc denotes the density of cationic cores in the metal bulk which can be estimated from the effective Wigner-Seitz radius or the cubic cell structure.? The dimensionless n̅ cc with respect to a 0 ^–3^ ranges from 0.001 to 0.01 for common metals (see supplementary note 9 in the Supporting Information.) n c,a ^bulk^ is the cation/anion density in the solution bulk that ranges from 1 mM to 1 M.? The ion size a ranges from 1 Å to 5 Å.? ϵ_op_ denotes the permittivity which is different between metal and electrolyte solution. For typical metals, ϵ_op_ is in the range from 2 ϵ_0_ to 9 ϵ_0_.? While for the electrolyte solution, ϵ_op_ ranges from 5 ϵ_0_ to 20 ϵ_0_ (organic solvent)? to 78.5 ϵ_0_ (aqueous solvent). ?,? Here we set n _ c,a _ ^bulk^ = 1 M, a = 0 and ϵ_op_ = ϵ_0_, which means neglecting the ion size and assuming ϵ takes the vacuum’s value without spatial distributions.

For the empirical parameters, d c,a is the equilibrium distance of cation/anion to the metal surface, which ranges from 2 Å to 10 Å. ?,?

D _ l _ is the depth of the Morse potential well that ranges from 0.01 to 0.1 eV. ?,? β_c,a_ is the coefficient that controls the width of the potential well to be adjusted for a good agreement with experiments. θ_T_ and μ_PG_ are empirical parameters to be determined by first-principle calculations or experiments.? We set d c,a = 4a 0, D _ l _ = 0.5/6 eV, β_c,a_ = 1, while θ_T_ and μ_PG_ are calibrated by the solutions of the Kohn–Sham–Poisson–Boltzmann model.

Results and Discussion

First, we present results of the Kohn–Sham-Poisson–Boltzmann model with basic setups, including the distributions of energy levels, wave functions, electrostatic potential, electron density, and comparison of convergence plots using different iteration schemes. Second, we compare the electrostatic potential and electron density calculated from the Kohn–Sham–Poisson–Boltzmann and DPFT models using the TFvW and Pauli–Gaussian functionals, respectively. Third, we compare the C dl values calculated from the Kohn–Sham–Poisson–Boltzmann and DPFT models using the TFvW and Pauli–Gaussian functionals, respectively. Fourthly, we compare the DPFT model and Kohn–Sham–Poisson–Boltzmann model using Perdew–Burke–Ernzerhof (PBE) functional for the exchange-correlation functional. Lastly, we improve the electrolyte solution model by accounting for finite ion size via a lattice-gas model and solvent polarization through a Langevin function and compare the results calculated from the DPFT model and the Kohn–Sham-Poisson–Boltzmann model.

A Glance at the Kohn–Sham–Poisson–Boltzmann

Model

Following the procedure shown in Figure, we obtain numerical solutions at μ̃_e_ = – 3.6 eV, as shown in Figure. Only one metal–solution interface is shown here. Figure(a) displays the first five energy levels ε_ i _ of electrons. The energy gap Δε_ i _ increases at higher ε_ i , which can be understood using the standard “particle in a box” model because the metal bulk approximates the scenario of an infinite potential well. In the “particle in a box” model, Δε i _ is proportional to 2i + 1.?

The first five normalized wave functions φ̅_ i _ are shown in Figure(b). As ε_ i _ increases, the period of φ̅_ i _ in metal decreases and the distributions become denser. Figure(c) displays the distribution of electrostatic potential ϕ. It is almost constant in the metal bulk and decreases to zero within a region of a few Å in the solution bulk. Figure(d) shows the distributions of dimensionless electron density n̅ e and subband’s contribution n̅ e (ε_ i _). n̅ e exhibits Friedel oscillations at the metal–solution interface, which is absent in DPFT or orbital-free DFT. Figure(e) displays the distributions of dimensionless cation density, n̅ c, and anion density, n̅ a, where minor difference is observed at the metal-solution interface. This means the system is nearly at the electroneutrality case with a controlled μ̃ e = −3.6 eV, namely, under the potential of zero charge condition.

Figure(f) compares the residuals among the fixed-point iteration, linear mixing method and predictor-corrector method. The convergence speeds are almost the same between the fixed-point iteration and predictor-corrector method, whereas the linear mixing method is slower. Increasing the mixing factor α accelerates convergence as the next-step solution will be closer the convergent value. Additionally, the linear mixing method exhibits some hysteresis, which arises from the next-step solution retaining information from the previous-step solution. Therefore, hysteresis is suppressed with increasing α from 0.05 to 0.5.

Comparison between Kohn–Sham–Poisson–Boltzmann

and DPFT Models in Terms of the Electron Density and Electrostatic Potential

We compare the electron density and electrostatic potential calculated from the Kohn–Sham–Poisson–Boltzmann and DPFT models using the TFvW and Pauli–Gaussian kinetic energy functionals. The comparison is conducted at μ̃ e = −3.57 eV, which is the potential of zero charge for the Kohn–Sham–Poisson–Boltzmann model at n̅ cc = 0.01, see supplementary note 10 in the Supporting Information.

DPFT with TFvW Kinetic Energy Functional

To compare the differences between the two approaches, we define two normalized variables, the normalized electron density n e,n = n̅ e/n̅ cc and the normalized electron density error Δn e,n = (n̅ e ^DPFT^ – n̅ e ^Kohn–Sham–Poisson–Boltzmann^)/n̅ cc.

Figure(a) shows the distributions of Δn e,n with increasing θ_T_ at E F = −3.57 eV. Three different peaks (peak 1: ∼−2.5 Å, peak 2: ∼−1 Å, peak 3: 0–1 Å) for Δn e,n near the metal/solution interface are observed. We notice that all three peaks change with increasing θ_ T . Specifically, the first two peaks are more pronounced while the third peak exhibits no clear trend with increasing θ T . Figure(b) shows the distributions of n e,n with increasing θ T . We notice that the electron spillover is more pronounced with increasing θ T _ due to the stronger gradient correction at the metal/solution interface (larger electron density gradient region). The TFvW kinetic functional with θ_ T _ = 0.7 gives a better approximation in terms of n e,n. Figure(c) shows the comparison of the electrostatic potential calculated from the Kohn–Sham–Poisson–Boltzmann and DPFT models with increasing θ_ T . The surface potential drop increases as θ T _ increases, which is caused by enhanced electron spillover as shown in Figure(b). Figure(d) shows the electrostatic potential calculated from the Kohn–Sham–Poisson–Boltzmann and the DPFT models in the region of the metal/solution interface, where θ_ T _ = 0.6 gives a more accurate approximation. Comparison of Figure(b) and (d) shows that the optimal value of θ_ T _ is different for n e,n and ϕ distributions. Moreover, we observe oscillations within the metal in the electrostatic potential calculated using the Kohn–Sham–Poisson–Boltzmann model, which are absent in the DPFT model. These oscillations arise from the oscillatory distributions of the electron density shown in Figure(a). Since $[eqn]$ , the oscillation of electron density will cause oscillatory distributions of ϕ­(x) within the metal.

Comparison between the DPFT model with the TFvW functional and Kohn–Sham–Poisson–Boltzmann model at μ̃ e = −3.57 eV in terms of (a) normalized electron density error Δn e,n = (n̅ e DPFT – n̅ e Kohn–Sham–Poisson–Boltzmann)/n̅ cc, (b) normalized electron density n e,n = n̅ e/n̅ cc, (c) electrostatic potential, (d) electrostatic potential at the region of the metal/solution interface. Parameters used in calculation are the same as those in Figure .

DPFT with Pauli–Gaussian Kinetic Energy Functional

Figure(a) shows the distributions of Δn e,n with increasing μ_PG_ at a fixed θ_ T _ = 5/3, a default value in the Pauli–Gaussian functional.? Three different peaks (peak 1: ∼−2.5 Å, peak 2: ∼−1 Å, peak 3: 0–1 Å) for Δn e,n near the metal/solution interface are also observed. The difference from the Figure (a) is that only the second and third peaks change with increasing μ_PG_ while the first peak keeps unchanged. This indicates that a large θ_ T _ will overestimate the electron density which cannot be corrected by the e ^–μ_PG_ s ^2^ ^ term. For the Pauli–Gaussian functional, $[eqn]$ overestimates the electron density (peak 1), compared with that from the Kohn–Sham–Poisson–Boltzmann model. Figure(b) shows the distributions of n e,n with increasing μ_PG_. We notice that the Pauli–Gaussian functional with μ_PG_ = 0.35 gives a more accurate approximation in terms of n e,n. Figure(c) shows the comparison of the electrostatic potential calculated from the Kohn–Sham–Poisson–Boltzmann and DPFT models with increasing μ_PG_. The surface potential drop decreases as μ_PG_ increases, which is consistent with the understanding that a larger μ_PG_ suppresses the electron spillover shown in Figure(b). Figure(d) shows the electrostatic potential calculated from the Kohn–Sham–Poisson–Boltzmann and DPFT models in the region of the metal–solution interface, where μ_PG_ = 0.25 gives a more accurate approximation in terms of the distribution of the electrostatic potential. Figure(b) and (d) also suggest different optimal values of μ_PG_ for n e,n and ϕ. The Kohn–Sham–Poisson–Boltzmann model exhibits oscillations in the electrostatic potential within the metal calculated, which are absent in the DPFT model with the Pauli–Gaussian kinetic energy functional. These oscillations also arise from the oscillatory distributions of the electron density shown in Figure(a).

Comparison between the DPFT model with the Pauli-Gaussian functional and Kohn–Sham–Poisson–Boltzmann model at μ̃ e = −3.57 eV in terms of (a) normalized electron density error Δn e,n = (n̅ e DPFT – n̅ e Kohn–Sham–Poisson–Boltzmann)/n̅ cc, (b) normalized electron density n e,n = n̅ e/n̅ cc, (c) electrostatic potential, (d) electrostatic potential at the region of the metal/solution interface. Parameters used in calculation are the same as those in Figure .

Comparison between Kohn–Sham–Poisson–Boltzmann

and DPFT Models in Terms of Double-Layer Capacitance Curves

In this part, we explain the calculation of the surface free charge and double-layer capacitance, and study the effect of θ_ T _ and μ_PG_ and the charge density of metal cationic cores on the double-layer capacitance.

Definitions of the Surface Free Charge and Double-Layer Capacitance

The surface free charge density for the EDL is defined as:?

$$[eqn]$$

where the second equal sign is due to the dimensionless definition.

Double layer capacitance (C dl) curve can be obtained by differentiating σ_free_ with electrochemical potential of electrons μ̃_e_, namely,

$$[eqn]$$

C dl exhibits a minimum extremum at the potential of zero charge condition in dilute solutions,? while a maximum extremum in highly concentrated solutions. ?,?

Effect of θ

T and μPG

Figure(a) shows the comparison of C dl calculated from Kohn–Sham–Poisson–Boltzmann and DPFT models with the TFvW functional. From the C dl curve calculated by the Kohn–Sham-Poisson–Boltzmann model, we determine the potential of zero charge to be approximately −3.6 eV, which is close to −3.57 eV calculated from the surface free charge (see supplementary note 5 in the Supporting Information). For the DPFT model described using the TFvW functional, the potential of zero charge becomes more negative with increasing θ_ T . We calibrate θ T _ = 0.7 in the DPFT model to capture the potential of zero charge calculated by the Kohn–Sham–Poisson–Boltzmann model. The DPFT model with θ_ T _ = 0.7 also captures well the electron spillover in Figure(b). However, C dl is overrated at less negative μ̃_e_, shown in Figure(a), which results from the overestimation of the electron spillover.

Comparison between DPFT model using (a) TFvW functional, (b) Pauli–Gaussian functional and Kohn–Sham–Poisson–Boltzmann model in terms of double-layer capacitance at n̅ cc = 0.01. Parameters used in the calculation are the same as those in Figure .

The DPFT model described using the Pauli–Gaussian functional exhibits poorer performance in predicting C dl shown in Figure(b), regardless of whether the system is positively or negatively charged. Especially at less negative μ̃_e_, a large C dl arises from the larger electron density error, though it predicts the potential of zero charge better at μ_PG_ = 0.25 in accordance with Figure(d).

To quantify the differences between the DPFT and Kohn–Sham–Poisson–Boltzmann models, we perform a comprehensive error analysis (detailed results are provided in supplementary note 11 in the Supporting Information). We notice that θ_ T _ = 0.7 and μ_PG_ = 0.25 are optimized parameters in terms of μ̃_e_ ^pzc^, whereas the DPFT models at θ_ T _ = 0.5 and μ_PG_ = 0.35 yield the minimum mean error in terms of C dl.

Effect of the Charge Density of Metal Cationic Cores

Metal cationic core density n̅ cc is a characteristic parameter of metals in our model. As shown in Figure and Figure, the C dl curve changes at a smaller n̅ cc in both Kohn–Sham–Poisson–Boltzmann and DPFT models. Figure(a) and (b) show that the overrated C dl calculated by the DPFT models is suppressed at less negative μ̃_e_ as n̅ cc decreases. This is because that the overestimation of the electron spillover is suppressed at a smaller n̅ cc. The potential of zero charge determined by the Kohn–Sham–Poisson–Boltzmann model shifts from ∼−3.6 eV to −3.3 eV as n̅ cc decreases from 0.01 to 0.005. Though the DPFT model described using the TFvW functional with θ_ T _ = 0.7 still fits the potential of zero charge prediction, θ_ T _ = 0.6 fits the C dl best. The Pauli–Gaussian functional still exhibits a poorer performance than the TFvW functional and predicts the potential of zero charge value best at μ_PG_ = 0.3. A quantitative error analysis is detailed in supplementary note 12 in the Supporting Information. We notice that θ_ T _ = 0.7 and 0.8, as well as μ_PG_ = 0.2 are optimized parameters in terms of μ̃_e_ ^pzc^, whereas the DPFT models at θ_ T _ = 0.5 and μ_PG_ = 0.35 have the minimum mean error in terms of C dl. From the comparison of C dl, we conclude that the TFvW functional provides a more accurate description of metal-solution interfaces than the Pauli–Gaussian functional, though the latter is reported to perform better in bulk metals and semiconductors.? This highlights the importance of examining orbital-free DFT in the context of metal–solution interfaces.

Comparison between DPFT model with (a) TFvW functional and (b) Pauli–Gaussian functional and Kohn–Sham–Poisson–Boltzmann model results in terms of double-layer capacitance at n̅ cc = 0.005. Other parameters used in the calculation are the same as those in Figure .

Effect of Exchange-Correlation Functional

We further compare the results calculated from the DPFT and Kohn–Sham–Poisson–Boltzmann models using an improved description of XC potential based on a general gradient approximation, namely, the Perdew–Burke–Ernzerhof (PBE) functional.

The exchange-correlation potential energy described by the PBE functional is defined as

$$[eqn]$$

where v XC ^PBE^ is the volumetric exchange-correlation energy including an exchange part and a correlation part:?

$$[eqn]$$

where v X is expanded as

$$[eqn]$$

with the volumetric exchange energy of a uniform electron gas: $[eqn]$ . θ_X_ is a gradient coefficient tuning the contribution of the gradient term in the exchange energy. Similarly, v C is expanded as,

$$[eqn]$$

where $[eqn]$ is a reduced density gradient and v C ^0^ is the volumetric correlation energy of a uniform electron gas, for which we use the interpolation relation of Perdew et al.,?

$$[eqn]$$ $$[eqn]$$

with α_1_ = 0.0310907, α_2_ = 0.21370, α_3_ = 7.5957, α_4_ = 3.5876, α_5_ = 1.6382, and α_6_ = 0.49294. θ_ C _ is a gradient coefficient in the correlation energy.

Replacing the Dirac–Wigner exchange-correlation functional with the PBE functional, we calculated the C dl from the DPFT model and Kohn–Sham–Poisson–Boltzmann model, respectively.

Figure(a) and (b) compares C dl calculated from Kohn–Sham–Poisson–Boltzmann and DPFT models with the TFvW functional and Pauli–Gaussian functional at n̅ cc = 0.01, respectively. For the DPFT model with the TFvW functional, θ_ T _ = 0.6 is a good approximation. The mean error of C dl defined as $[eqn]$ with N tot being total number of C dl ^DPFT^ data points, shown in Figure(e), exhibits a minimum value at θ_ T _ = 0.6, in agreement with the trend observed in Figure(a). The DPFT model with the Pauli–Gaussian functional exhibits poorer performance in predicting C dl, as shown in Figure(b), though it predicts the PZC better at μ_PG_ = 0.25, as shown in Figure(d). The mean error of C dl shown in Figure(f) decreases with increasing μ_PG_, in accordance with Figure(b).

Comparison between DPFT model with TFvW functional and Pauli–Gaussian functional and Kohn–Sham–Poisson–Boltzmann model, in terms of double-layer capacitance (a, b), μ̃e PZC (c, d), and mean error (e, f), respectively. n̅ cc = 0.01, other parameters are the same as those in Figure . The exchange-correlation potential is described using the PBE functional.

We note that the calibrated parameters depend on the choice of the XC potential, with the minimum error at θ_T_ = 0.5 for the Dirac-Wigner functional and at θ_T_ = 0.6 for the PBE functional, whereas the theoretical framework and benchmarking workflow are general.

Effect of Finite Ion Size

The finite ion size and solvent polarization have been known to importantly influence the EDL, which have been neglected in the foregoing comparison. Next, we describe the electrolyte solution more accurately by accounting for finite ion size via a lattice-gas model and solvent polarization through a Langevin function, as shown in eq S19, and perform a comprehensive comparison between the Kohn–Sham–Poisson–Boltzmann and DPFT models using the TFvW kinetic energy functional. In both models, the exchange-correlation potential is described using the PBE functional.

Figure(a)–(f) show the comparison of the electrostatic potential, normalized electron density, cation density, anion density, solvent density, and relative dielectric permittivity calculated from the Kohn–Sham–Poisson–Boltzmann, and DPFT models with increasing θ_ T , respectively. Relative to the earlier results in Figure, where ion size and solvent effects were neglected, the present results show pronounced maxima in the cation and anion densities, and the spatial variation of the dielectric permittivity closely follows that of the solvent density. For electrostatic potential and electron density, θ T _ = 0.3 gives a more accurate approximation inside the metal, whereas θ_ T _ = 0.5 performs better outside the metal. For cation and anion densities, θ_ T _ = 0.4 gives more accurate approximations. A similar trend is observed for the solvent density and relative permittivity, for which θ_ T _ = 0.4 also performs better. Figure(g) shows the surface charge with increasing θ_ T . For negatively charged surfaces, θ T _ = 0.3 and θ_ T _ = 0.4 give more accurate approximations, whereas θ_ T _ = 0.5 performs better for positively charged surfaces. Figure(h) presents the comparison of C dl with increasing θ_ T , where the DPFT results are overall slightly larger than those obtained from the Kohn–Sham–Poisson–Boltzmann model. The C dl exhibits a camel shape compared with earlier results in Figure and Figure, which is caused by the effect of finite ion size.? Figure(i) compares μ̃_e ^pzc^, determined from DPFT and the Kohn–Sham–Poisson–Boltzmann models, showing that θ_ T _ = 0.4 gives a more accurate approximation.

*Comparison between the DPFT model with the TFvW functional and Kohn–Sham–Poisson–Boltzmann model at μ̃e = −2.45 eV in terms of (a) electrostatic potential ϕ, (b) normalized electron density n e,n, (c) normalized cation density n c,n, (d) normalized anion density n a,n, (e) normalized solvent density n s,n, (f) relative dielectric permittivity ϵ/ϵ0, (g) surface free charge σfree, (h) double-layer capacitance C dl, and (i) potential of zero charge μ̃e pzc. Parameters used in the calculation are as follows: n̅ cc = 0.01, r s = 2.75 Å, r a = 2.7 Å, r c = 3.6 Å, n c,a bulk = 0.1 M, n s bulk = 55.6 M, p s = 1.8 D, ϵop=ϵopM+(ϵopS−ϵopM)2⁡erfc(−βop(|X|−LM)) , ϵop M = 3ϵ0, ϵop

S = 6ϵ0, D s = 0.5 eV, d s = 2a 0, βs = 1, L M = 20, L S = 120, other parameters are the same as those in Figure . The exchange-correlation potential is described using the PBE functional.*

Conclusions

In this article, we have presented a Kohn–Sham–Poisson–Boltzmann theory to simulate one-dimensional metal–solution interfaces under constant-potential conditions and use it as a benchmark for orbital-free density-potential functional theory (DPFT). The Kohn–Sham–Poisson–Boltzmann model exhibits Friedel oscillations of electron density which are absent in the current DPFT model. Empirical parameters used in the DPFT models are calibrated with the solutions of the Kohn–Sham–Poisson–Boltzmann model. The comparison of electron density, electrostatic potential and double-layer capacitance shows that the Thomas–Fermi–von Weizsäcker kinetic functional provides a more accurate description of metal-solution interfaces than the Pauli–Gaussian functional, in contrast with their performance for bulk materials. The double-layer capacitance is overrated at less negative electrochemical potential of electrons for both functionals, which results from the overestimation of the electron spillover. Decreasing the density of cationic cores can suppress the electron spillover at the metal surface. The calibrated empirical parameters in DPFT models depend on the choice of the XC potential, whereas the overall trends are similar. Considering the effects of finite ion size and solvent polarization changes the ion-density distributions, dielectric permittivity, the shape of C dl, and the calibrated empirical parameter in the DPFT models. It should be noted that our assumption of a laterally uniform metal/solution interface reduces the problem to one dimension, which is a standard and widely used setup for evaluating kinetic energy functionals in planar systems. Under this ideal planar geometry, the density varies only along the surface normal, allowing us to isolate and compare the intrinsic behavior of the TFvW and Pauli–Gaussian functionals without additional geometric effects. Extending such benchmarking from one-dimensional planar to two-dimensional curved surfaces or nanoscale structures would require accounting for the additional geometric complexity, including curvature-induced anisotropy in the density gradients and spatial variations parallel to the interface.

Supplementary Material