ESM method¶
The effective screening medium (ESM) method treats the electrostatics of charged surfaces and interfaces, and of surfaces under an applied bias, in plane-wave first-principles calculations. It removes the periodic boundary condition perpendicular to the slab (the \(z\) direction) and solves the Poisson equation with a Green's function under the boundary conditions set by a continuum medium (vacuum, ideal metal, etc.) placed outside the slab. The Kohn-Sham equation is still solved under the usual periodic boundary conditions, so only the electrostatic part of an existing plane-wave pseudopotential code needs to be changed. Calculations at a fixed electrode potential are covered in Constant-μ method.
Charged surfaces and periodic boundary conditions¶
From the mid-1990s, driven mainly by fuel-cell research, there was a growing need for first-principles calculations of electrode/electrolyte interfaces under an applied voltage (electrochemical interfaces). An electric double layer forms at such an interface, and an electric field arises between the excess charge on the electrode surface and the counter charge that compensates it.
Plane-wave calculations impose periodic boundary conditions that repeat the cell in three dimensions. If one tries to represent an electric double layer under these conditions, the electrostatic potential becomes discontinuous at the cell boundary. If the cell carries a net charge, the electrostatic energy diverges; adding a uniform background charge to prevent this removes the physical meaning of the absolute electron energy levels. Conventional approaches based on a sawtooth potential or a uniform background charge were therefore inadequate as models of electrochemical interfaces.
The ESM method proceeds as follows.
- Keep the periodic boundary conditions in the surface-parallel (\(xy\)) directions and remove them in the \(z\) direction.
- Place a screening medium characterized by a dielectric permittivity \(\varepsilon(\boldsymbol{r})\) outside the slab and impose boundary conditions appropriate to the model (isolated slab, biased surface, capacitor, etc.).
- The screening medium carries the charge that compensates the excess charge of the slab, so the whole system, including the medium, is always neutral.
Total-energy functional and the generalized Poisson equation¶
Let \(\rho_\mathrm{e}\) be the electron charge density, \(\rho_\mathrm{I}\) the ionic charge density, and \(\rho_\mathrm{tot}=\rho_\mathrm{e}+\rho_\mathrm{I}\) the total charge density. Following the convention of plane-wave pseudopotential codes, we take the electron charge as positive. Treating the electrostatic potential \(V(\boldsymbol{r})\) as an independent variable, we define the total energy by the functional
where \(T_\mathrm{s}\) is the kinetic energy, \(E_\mathrm{xc}\) the exchange-correlation energy, and \(\varepsilon(\boldsymbol{r})\) a position-dependent relative permittivity. The variation with respect to \(V\) is
The boundary conditions make the surface integral vanish, so the stationarity condition \(\delta E/\delta V=0\) gives the generalized Poisson equation and the corresponding equation for the Green's function:
The electrostatic potential is then
Substituting the solution of the Poisson equation into Eq. \(\eqref{eq:functional}\) and integrating by parts with the identity \(\nabla\cdot(\varepsilon V\nabla V)=V\nabla\cdot(\varepsilon\nabla V)+\varepsilon|\nabla V|^2\), the electrostatic terms reduce to \(\frac{1}{2}\int d\boldsymbol{r}\,\rho_\mathrm{tot}V\), and the total energy recovers its familiar form:
The electrostatic term splits into electron–electron (Hartree), electron–ion, and ion–ion contributions:
In other words, the ESM method keeps the form of the DFT total energy and replaces the electrostatic interaction \(1/|\boldsymbol{r}-\boldsymbol{r}'|\) by the Green's function \(G(\boldsymbol{r},\boldsymbol{r}')\).
Varying Eq. \(\eqref{eq:functional}\) with respect to the orbitals gives the Kohn-Sham equation:
The two Euler-Lagrange equations are handled separately.
- Poisson equation: solved with a Green's function obtained analytically in advance for the boundary condition of the model. The Green's function is hard-coded in the program, so no differential equation is solved at run time and the extra cost is minimal.
- Kohn-Sham equation: solved with the standard periodic solver. Because the resulting electrostatic potential is not periodic in \(z\), it is connected smoothly near the cell boundary, where the electron density is negligible. This treatment is valid as long as the electron density of the slab has decayed sufficiently before reaching the cell boundary.
Pseudopotentials are generated for isolated atoms in vacuum, so the nuclei are placed in the region where \(\varepsilon=1\).
Laue representation¶
Consider the case where \(\varepsilon\) depends only on \(z\). Because periodic boundary conditions are imposed in the \(xy\) directions, the Poisson equation can be treated in the Laue representation, Fourier-transformed in the in-plane directions:
Transforming the charge density and the Green's function in the same way gives the Poisson equation in the Laue representation (\(g_\parallel=|\boldsymbol{g}_\parallel|\)):
The electrostatic potential is
The problem has been reduced to a one-dimensional second-order ordinary differential equation in \(z\) for each \(g_\parallel\), which can always be solved under appropriate boundary conditions.
In what follows, the electrons of the slab are confined to \(-z_0\le z\le z_0\) (a computational cell of length \(L=2z_0\)), and the boundary of the screening medium is at \(z=\pm z_1\) (\(z_1\ge z_0\)).
Construction of the Green's function¶
For the inhomogeneous Sturm-Liouville equation
the Green's function is fixed by the following conditions.
- For \(z\ne z'\), it satisfies the homogeneous equation \(\hat{L}G=0\).
- At \(z=a\) and \(z=b\), it satisfies the boundary conditions imposed on \(y\).
- It is continuous at \(z=z'\).
- Its derivative jumps at \(z=z'\) by \(\left.\partial_z G\right|_{z'+0}-\left.\partial_z G\right|_{z'-0}=-1/p(z')\).
Let \(\psi_\mathrm{L}\) be the homogeneous solution that satisfies the left boundary condition and \(\psi_\mathrm{R}\) the one that satisfies the right boundary condition. The Wronskian combination \(p(z)\left[\psi_\mathrm{L}\psi_\mathrm{R}'-\psi_\mathrm{L}'\psi_\mathrm{R}\right]=A\) is a constant, and conditions 3 and 4 give
This construction yields the Green's function for any boundary condition.
In the region where \(\varepsilon=1\), Eq. \(\eqref{eq:poisson-laue}\) becomes
which has the form of Eq. \(\eqref{eq:sturm}\) with \(p=1/4\pi\), \(q=-g_\parallel^2/4\pi\), and \(f=\rho_\mathrm{tot}\). The homogeneous solutions are linear combinations of \(e^{\pm g_\parallel z}\), so the boundary conditions differ only in the choice of \(\psi_\mathrm{L}\) and \(\psi_\mathrm{R}\). A region of infinite permittivity represents an ideal metal (electrode) with a constant potential inside, which we take as the energy reference (zero).
Boundary conditions and Green's functions¶
The three representative boundary conditions of the ESM method are:
| Boundary condition | Setup | Typical use |
|---|---|---|
| bc1 | vacuum/slab/vacuum | neutral surfaces, polar surfaces, work functions |
| bc2 | metal/slab/metal | slab in a capacitor, voltage applied between electrodes |
| bc3 | vacuum/slab/metal | STM tip or gate electrode, electrochemical interface with a counter electrode |
bc1: vacuum/slab/vacuum¶
The boundary condition is
The homogeneous solutions \(\psi_\mathrm{L}=a\,e^{g_\parallel z}\) and \(\psi_\mathrm{R}=b\,e^{-g_\parallel z}\) give \(A=-2abg_\parallel/4\pi\), and
This is simply the Laue representation of the bare Coulomb interaction \(1/|\boldsymbol{r}-\boldsymbol{r}'|\). Nothing screens the charge, so the Green's function diverges at \(g_\parallel=0\). Consequently bc1 applies only to neutral systems: for a charged slab, \(\partial_zV(0,z)\) does not vanish for at least one of \(z\to\pm\infty\), the potential grows linearly with \(z\), and the electrostatic energy diverges.
With bc1 the electrostatic potential is always referenced to the vacuum level (at infinity), so the work function follows directly. Even when the slab carries a dipole moment and the potential has a long tail, no dipole correction of the kind used in repeated-slab calculations is required.
bc2: metal/slab/metal¶
Ideal metals occupy \(|z|\ge z_1\) on both sides. Without an applied voltage, the boundary condition is
The homogeneous solutions \(\psi_\mathrm{L}=a\sinh[g_\parallel(z+z_1)]\) and \(\psi_\mathrm{R}=b\sinh[g_\parallel(z-z_1)]\) give \(A=abg_\parallel\sinh(2g_\parallel z_1)/4\pi\), and
for \(|z|,|z'|<z_1\). The first term is the bare Coulomb term; the second is the contribution of the image charges induced on the two electrodes. Expanding the second term gives
which is an infinite series of image charges of the same sign at \(z=z'\pm(4m+4)z_1\) and of opposite sign at \(z=-z'\pm(4m+2)z_1\).
bc3: vacuum/slab/metal¶
An ideal metal occupies only \(z\ge z_1\). The boundary condition is
The homogeneous solutions \(\psi_\mathrm{L}=a\,e^{g_\parallel z}\) and \(\psi_\mathrm{R}=b\sinh[g_\parallel(z-z_1)]\) give \(A=abg_\parallel e^{g_\parallel z_1}/4\pi\), and
for \(z,z'<z_1\). Unlike bc2, there is only one image charge, at \(z=2z_1-z'\). This boundary condition models an electrochemical cell consisting of an electrode (the slab) and a counter electrode (the ESM metal), and it is used together with the Constant-μ method.
The zero in-plane wave-vector component¶
At \(g_\parallel=0\), Eq. \(\eqref{eq:green-laue}\) becomes the one-dimensional Poisson equation \(\partial_z^2G=-4\pi\delta(z-z')\), whose solution is \(-2\pi|z-z'|\) plus a linear function. Imposing each boundary condition gives
The bc2 and bc3 expressions also agree with the \(g_\parallel\to0\) limits of Eqs. \(\eqref{eq:g-bc2}\) and \(\eqref{eq:g-bc3}\). In bc2 and bc3 the metal electrodes carry the counter charge, so the \(g_\parallel=0\) component stays finite even for a charged slab. The constant in bc1 is arbitrary but does not contribute to the total energy of a neutral system.
Real-space form and image charges¶
Writing the in-plane distance as \(R=|\boldsymbol{r}_\parallel-\boldsymbol{r}'_\parallel|\) and using the fact that \(\frac{4\pi}{2g_\parallel}e^{-g_\parallel|\zeta|}\) is the two-dimensional Fourier transform of \(1/\sqrt{R^2+\zeta^2}\), we obtain the Green's functions in real space:
Each has the form "bare Coulomb term plus image-charge terms", which shows that the electrostatics of the ESM method is the classical method of images. Field lines do not enter the ESM metal.
Terms of the total energy¶
In a plane-wave pseudopotential code, the terms that contain the Green's function must be modified: the ion–ion interaction (Ewald sum), the local pseudopotential, and the Hartree potential. The nonlocal pseudopotential, exchange-correlation, and kinetic terms are unchanged. In codes such as OpenMX and SIESTA, which compute the electrostatic potential from the total charge density in one step, only that step needs to be modified.
Hartree potential¶
Expand the electron density in the cell as \(\rho_\mathrm{e}(\boldsymbol{g}_\parallel,z)=\sum_{g_z}\rho_\mathrm{e}(\boldsymbol{g}_\parallel,g_z)e^{ig_zz}\). The Hartree potential is
The Green's function is a sum of exponentials in \(z'\), so the integral can be done analytically. For bc1 (\(g_\parallel\ne0\)),
The first term has the same form as the usual three-dimensional Hartree term \(4\pi\rho_\mathrm{e}(\boldsymbol{G})/G^2\); the second is a boundary term arising from the finite integration range. For bc2 and bc3 an image-charge term is added. For bc3, for example, it reads
The \(g_\parallel=0\) component is obtained in the same way from Eqs. \(\eqref{eq:g0-bc1}\)–\(\eqref{eq:g0-bc3}\).
Local pseudopotential¶
The electron–ion interaction is split into a long-range part generated by Gaussian pseudo-ion charges \(\rho_\mathrm{g}\), a short-range local part, and a nonlocal part:
The short-range and nonlocal parts are nonzero only near the nuclei and do not reach the ESM region, so they induce no image charges and need no modification. Only the long-range part \(V_\mathrm{loc}^\mathrm{long}(\boldsymbol{r})=\int d\boldsymbol{r}'\,G(\boldsymbol{r},\boldsymbol{r}')\rho_\mathrm{g}(\boldsymbol{r}')\) has to be changed.
Let the Gaussian charge of atom \(\nu\) (valence \(Z_\nu\), position \(\boldsymbol{\tau}_\nu\)) be \(\rho_\mathrm{g}^\nu(\boldsymbol{r})=-Z_\nu(\alpha/\pi)^{3/2}e^{-\alpha|\boldsymbol{r}-\boldsymbol{\tau}_\nu|^2}\) (negative because the electron charge is taken as positive). It is a product of an in-plane Gaussian and a Gaussian in \(z\), so the \(z'\) integral can be done analytically. For \(g_\parallel\ne0\),
where \(S_0\) is the area of the in-plane unit cell. \(F_\alpha\) is the bare Coulomb term \(G^\mathrm{bc1}\) averaged over the Gaussian charge, and \(G^\mathrm{i}\) is the image-charge term of the Green's function (zero for bc1, the second term of Eq. \(\eqref{eq:g-bc2}\) for bc2, and the second term of Eq. \(\eqref{eq:g-bc3}\) for bc3). As long as the Gaussian charge does not overlap the ESM region, its image coincides with that of a point charge, so the image-charge term need not be averaged over the Gaussian.
The short-range part \(v_\mathrm{loc}^\mathrm{short}=v_\mathrm{loc}-v_\mathrm{loc}^\mathrm{long}\) is Fourier-transformed numerically as before. With ESM, the analytic Fourier transform of the long-range part, which a standard calculation adds in reciprocal space, is omitted; Eq. \(\eqref{eq:vloc-long}\) is added in the Laue representation instead.
Ion–ion interaction (Ewald sum)¶
Treating the ions as point charges, the ion–ion interaction is
where \(\boldsymbol{R}_{\parallel\mu}=\boldsymbol{R}_\parallel+\boldsymbol{\tau}_\mu\) is the position of atom \(\mu\) in unit cell \(\boldsymbol{R}_\parallel\) and \(E_\mathrm{self}\) is the self-interaction. As in the standard Ewald method, adding and subtracting auxiliary Gaussian charges splits the energy into real-space and reciprocal-space sums, \(E_\mathrm{I\text{-}I}=E_\mathrm{R}+E_\mathrm{G}-E_\mathrm{self}\).
The real-space sum \(E_\mathrm{R}\) involves the difference "point charge minus Gaussian charge". Both have the same image charges, so the image contributions cancel and only the bare Coulomb term remains. The real-space sum is therefore identical to the standard Ewald sum:
The prime excludes the term with \(\boldsymbol{R}_\parallel+\boldsymbol{\tau}_\mu-\boldsymbol{\tau}_\nu=0\). The reciprocal-space sum is expressed with in-plane reciprocal lattice vectors \(\boldsymbol{G}_\parallel\) and real space in \(z\); its \(\boldsymbol{G}_\parallel\ne0\) part is
with the same functions \(F_\alpha\) and \(G^\mathrm{i}\) as for the local pseudopotential. For the \(\boldsymbol{G}_\parallel=0\) part, \(|z-z'|\) in Eqs. \(\eqref{eq:g0-bc1}\)–\(\eqref{eq:g0-bc3}\) is averaged over the Gaussian,
and the result is \(E_\mathrm{G}^{(\boldsymbol{G}_\parallel=0)}=\frac{1}{2S_0}\sum_{\mu\nu}Z_\mu Z_\nu\,G(0,\tau_{\mu z},\tau_{\nu z})\). The self-interaction term is the usual \(E_\mathrm{self}=\sum_\mu Z_\mu^2\sqrt{\alpha/\pi}\).
Forces on atoms¶
The forces on atoms are the analytic derivatives of the total energy, \(\boldsymbol{F}_\mu=-\partial E/\partial\boldsymbol{\tau}_\mu\). ESM changes only the contributions from the Ewald term and the long-range local pseudopotential term, which contain the Green's function. The cost of the force calculation is almost the same as in a repeated-slab calculation, and the forces are suitable for molecular dynamics that conserves the total energy accurately.
Applying excess charge and bias¶
In the ESM method, an excess charge is introduced by changing the number of electrons from its neutral value. With boundary conditions that include metal electrodes (bc2, bc3), the counter charge appears as image charge on the electrode surface, and the whole system including the medium remains neutral automatically. No uniform background charge is needed, and the reference of the electrostatic potential (electrode potential = 0) is preserved.
The bias corresponding to an excess charge is the difference between the Fermi levels of the charged and neutral systems:
where \(\Delta N\) is the number of excess electrons. In the Al/Si(111) calculations of Otani & Sugino (2006), \(\Delta\mu\) was nearly proportional to the induced charge. For bc2 the charge is induced on both faces of the slab, so the slope is about half that for bc3. For a surface about 10 Å from the electrode, a 1 V bias induced roughly one electron per 100 surface sites.
With bc2, a potential difference \(V_0\) can also be applied directly between the two electrodes. With the boundary conditions \(V(\boldsymbol{g}_\parallel,z_1)=0\) and \(V(\boldsymbol{g}_\parallel,-z_1)=V_0\,\delta_{\boldsymbol{g}_\parallel,0}\), the solution is the solution \(V^\mathrm{bc2}\) obtained with Eq. \(\eqref{eq:g-bc2}\) plus a uniform-field term:
Calculations that fix the number of electrons and obtain the bias as a result (constant-charge calculations) converge quickly and stably in the self-consistent cycle. Otani et al. (2008) used this approach, adding excess electrons in a bc3-type setup, to perform first-principles molecular dynamics of the Pt/water interface under bias. Calculations that vary the number of electrons at a target potential are covered in Constant-μ method.
Computational cost and scope¶
The main features of the ESM method are:
- Various boundary conditions can be imposed on surfaces and interfaces.
- Polar surfaces and surfaces with a dipole are handled without corrections.
- It is easy to add to existing plane-wave pseudopotential codes. Only three parts change: the Hartree potential, the local pseudopotential, and the ion–ion interaction.
- The computational time is comparable to that of the conventional method. The usual Hartree term costs \(O(N^3)\) operations over three-dimensional reciprocal lattice vectors; in ESM, each in-plane reciprocal lattice vector (\(N^2\) of them) requires a sum over \(g_z\), a one-dimensional FFT, and a loop over \(z\), giving \(N^2(N+N\log N+N)=(2+\log N)N^3\) operations. The order is unchanged; only the prefactor increases.
The screening medium is not limited to an ideal metal or vacuum. Otani & Sugino (2006) also presented a boundary condition that combines a solvent of finite relative permittivity with an electrolyte ion distribution from modified Poisson-Boltzmann theory. Hamada et al. (2013) generalized the original model, in which the permittivity changes stepwise, to a smooth permittivity profile with Debye-Hückel-type screening. Treating the electrolyte with an integral-equation theory of liquids is covered in ESM-RISM method.
References¶
- M. Otani and O. Sugino, "First-principles calculations of charged surfaces and interfaces: A plane-wave nonrepeated slab approach," Phys. Rev. B 73, 115407 (2006). doi
- M. Otani, I. Hamada, O. Sugino, Y. Morikawa, Y. Okamoto, and T. Ikeshoji, "Structure of the water/platinum interface—a first principles simulation under bias potential," Phys. Chem. Chem. Phys. 10, 3609 (2008). doi
- N. Bonnet, T. Morishita, O. Sugino, and M. Otani, "First-Principles Molecular Dynamics at a Constant Electrode Potential," Phys. Rev. Lett. 109, 266101 (2012). doi
- I. Hamada, O. Sugino, N. Bonnet, and M. Otani, "Improved modeling of electrified interfaces using the effective screening medium method," Phys. Rev. B 88, 155427 (2013). doi