How is the derivative discontinuity related to steps in the exact Kohn-Sham potential?

TL;DR

This paper analytically links the derivative discontinuity and steps in the exact Kohn-Sham potential, demonstrating their common origin and proposing a method to extract the fundamental gap from these features.

Contribution

It establishes the fundamental relationship between derivative discontinuity and steps in the exact xc potential through analytical derivation and numerical validation.

Findings

Derived the common origin of potential steps and derivative discontinuity.

Validated the relationship with exact numerical solutions for simple systems.

Proposed a method to extract the fundamental gap from step structures.

Abstract

The reliability of density-functional calculations hinges on accurately approximating the unknown exchange-correlation (xc) potential. Common (semi-)local xc approximations lack the jump experienced by the exact xc potential as the number of electrons infinitesimally surpasses an integer, and the spatial steps that form in the potential as a result of the change in the decay rate of the density. These features are important for an accurate prediction of the fundamental gap and the distribution of charge in complex systems. Although well-known concepts, the exact relationship between them remained unclear. In this Letter, we establish the common fundamental origin of these two features of the exact xc potential via an analytical derivation. We support our result with an exact numerical solution of the many-electron Schroedinger equation for a single atom and a diatomic molecule in one dimension. Furthermore, we propose a way to extract the fundamental gap from the step structures in the potential.

Equations

$$\Delta =I-A-({\epsilon }^{\mathrm{lu}}-{\epsilon }^{\mathrm{ho}}),$$

$$S=I_R-I_L+{\epsilon }_R^{\mathrm{ho}}-{\epsilon }_L^{\mathrm{ho}}$$

$$n(x;N)=(1-\delta )\cdot n(x;N_0)+\delta \cdot n(x;N_0+1),$$

Full text

Density functional theory (DFT) Hohenberg and Kohn (1964), in the Kohn-Sham (KS) approach Kohn and Sham (1965), is widely used for simulating many-electron systems Parr and Yang (1989); R.M. Dreizler and E.K.U. Gross (1990); Fiolhais et al. (2003); Engel and Dreizler (2011); Burke (2012); Martin (2004); Kaxiras (2003); Cramer (2004); Sholl and Steckel (2011). The accuracy of density-functional calculations hinges on approximating the unknown exchange-correlation (xc) energy term, $E_{\mathrm{xc}}[n]$. While numerous successful approximations exist Perdew and Wang (1992); Becke (1988); Lee et al. (1988); Perdew et al. (1996a, 2008); Sun et al. (2015); Schimka et al. (2011); Becke (1993); Stephens et al. (1994); Perdew et al. (1996b); Ernzerhof and Scuseria (1999); Adamo and Barone (1999), they often lack the discontinuous nature of the derivative of $E_{\mathrm{xc}}[n]$ with respect to electron number, $N$, at integer $N$ (derivative discontinuity (DD) Perdew et al. (1982); Cohen et al. (2012); Baerends et al. (2013); Mori-Sánchez and Cohen (2014); Mosquera and Wasserman (2014a, b)).

The DD is essential for exactly describing the fundamental gap – a feature of central importance for any material. When relying on the KS eigenvalues, approximate functionals (e.g., Perdew and Wang (1992); Becke (1988); Lee et al. (1988); Perdew et al. (1996a)) lacking the DD underestimate this quantity by $\sim 50\%$ Tran et al. (2007); Tran and Blaha (2009); Eisenberg and Baer (2009); Schimka et al. (2011); Chan and Ceder (2010); Lucero et al. (2012). Furthermore, these functionals may qualitatively fail in the dissociation limit, by predicting spurious fractional charges on the atoms of a stretched diatomic molecule  Perdew et al. (1982); Ossowski et al. (2003); Dutoi and Head-Gordon (2006); Gritsenko and Baerends (2006); Mori-Sánchez et al. (2006); Ruzsinszky et al. (2006); Vydrov and Scuseria (2006); Perdew et al. (2007); Vydrov et al. (2007); Kraisler and Kronik (2015), violating the principle of integer preference Perdew (1990). This indicates that common approximations may also fail to describe charge transfer in molecules and materials Perdew and Smith (1984); Tozer (2003); Maitra (2005); Toher et al. (2005); Koentopp et al. (2006); Ke et al. (2007); Hofmann and Kümmel (2012); Nossa et al. (2013); Fuks (2016).

One manifestation of the DD is the emergence of a spatially uniform ‘jump’, $\Delta$, in the level of the KS potential, $v_{\mathrm{s}}(\mathbf{r})$, as $N$ infinitesimally surpasses an integer value of $N_{0}$ by $\delta$: $\Delta=\lim_{\delta\rightarrow 0^{+}}v_{\mathrm{s}}(\mathbf{r};N_{0}+\delta)-v_{\mathrm{s}}(\mathbf{r};N_{0}-\delta)$ Perdew and Levy (1983); Sham and Schlüter (1983); Kohn (1986); Godby et al. (1987, 1988); Chan (1999); Allen and Tozer (2002); Mundt and Kümmel (2005); Teale et al. (2008); Sagvolden and Perdew (2008); Gori-Giorgi and Savin (2008); Yang et al. (2012); Gould and Toulouse (2014); Görling (2015). $\Delta$ originates from the piecewise-linearity of the total energy, $E(N)$ Perdew et al. (1982); Cohen et al. (2008); Mori-Sánchez et al. (2008, 2009); Stein et al. (2012); Kronik et al. (2012); Atalla et al. (2016), which implies a stair-step structure of the highest occupied (ho) KS eigenvalue, $\varepsilon^{\mathrm{ho}}(N)$, with discontinuities at integer $N$. In particular, for a neutral system of $N_{0}$ electrons, infinitesimally below $N_{0}$, $\varepsilon^{\mathrm{ho}}(N_{0}^{-})$ equals $-I$, the negative of the ionization potential (IP), and infinitesimally above $N_{0}$, $\varepsilon^{\mathrm{ho}}(N_{0}^{+})$ equals $-A$, the negative of the electron affinity (EA) Perdew et al. (1982); Levy et al. (1984); Perdew and Levy (1997); Harbola (1998, 1999); Yang et al. (2012). To enforce this behavior, the exact $v_{\mathrm{s}}(\mathbf{r})$ has to experience a discontinuous jump Perdew and Levy (1983); Sham and Schlüter (1983):

[TABLE]

where $\varepsilon^{\mathrm{ho}}$ and $\varepsilon^{\mathrm{lu}}$ are the ho and the lowest unoccupied (lu) KS eigenvalues at and infinitesimally below $N_{0}$ (from here on, the argument $N_{0}^{-}$ is suppressed for brevity). While this constant shift in $v_{\mathrm{s}}(\mathbf{r})$ does not affect the electron density, $n(\mathbf{r})$, it is vital to accurately predict the fundamental gap Perdew and Levy (1983); Perdew (1985a); Bylander and Kleinman (1996); Städele et al. (1997, 1999); Magyar et al. (2004); Rinke et al. (2005); Becke and Johnson (2006); Grüning et al. (2006a, b); Tran et al. (2007); Rinke et al. (2008); Tran and Blaha (2009); Kuisma et al. (2010); Chan and Ceder (2010); Zheng et al. (2011); Kronik et al. (2012); Chai and Chen (2013); Armiento and Kümmel (2013); Laref et al. (2013); Kraisler and Kronik (2013, 2014); Kraisler et al. (2015); Atalla et al. (2016).

In addition, the existence of the DD implies that the KS potential may form a ‘plateau’ – a constant increase in the level of $v_{\mathrm{s}}(\mathbf{r})$ in a given region – to correctly distribute charge throughout the system Perdew (1985b, 1990); Krieger et al. (1992); Sagvolden and Perdew (2008); Karolewski et al. (2009); Tempel et al. (2009); Makmal et al. (2011); Fuks et al. (2011); Hofmann and Kümmel (2012); Nafziger and Wasserman (2015); Gould and Hellgren ; Li et al. (2015); Komsa and Staroverov (2016); Kohut et al. (2016). At the edge of a plateau, $v_{\mathrm{s}}(\mathbf{r})$ forms a spatial step.

In a stretched diatomic molecule $L\cdots R$, the height of the step, $S$, at the interface between the two atoms can be deduced Hodgson et al. (2016) from the change in the exponential decay of the density. When the density decay from the left atom, which has the analytic form $n(\mathbf{r})\propto e^{-2\sqrt{2I_{L}}|\mathbf{r}|}$ Perdew et al. (1982); Perdew and Levy (1997), meets the density decay from the right atom ($e^{-2\sqrt{2I_{R}}|\mathbf{r}|}$), a step of height

[TABLE]

forms in $v_{\mathrm{s}}(\mathbf{r})$ at this point Almbladh and von Barth (1985); van Leeuwen et al. (1995); Gritsenko and Baerends (1996); Helbig et al. (2009); Hellgren and Gross (2012); Yang et al. (2014); Hodgson et al. (2014, 2016); Benítez and Proetto (2016).

The exact relationship between the two aforementioned manifestations of the DD, $\Delta$ and $S$, remains unclear. These two quantities are usually treated as unrelated, because the jump in $v_{\mathrm{s}}(\mathbf{r})$ is a function of $N$, whereas a plateau is a function of space. Furthermore, the EA and the lu energy contribute to $\Delta$, while they are absent from $S$.

A complete understanding of how the presence of steps in the KS potential account for the DD is essential for our ability to predict the fundamental gap of many-electron systems within DFT. Hence, in this Letter we directly address this problem by identifying the general mechanism that gives rise to $\Delta$ and $S$. We explore the properties of the plateau in the exact KS potential as a function of $N$ and show how $\Delta$ can be deduced from the step structure of the KS potential as $N$ decreases to an integer. This is done first via an analytical derivation, relying on fundamental properties of many-electron systems, and supported with an exact numerical solution of the many-electron Schrödinger equation in one dimension, for a single atom and a diatomic molecule. Modelling in 1D is necessary to solve the Schrödinger equation exactly; yet, the principles demonstrated may be generalized to 3D systems  111See the Supplemental Material at [added by journal] for technical details.

.

To establish the relationship between $\Delta$ and $S$, we propose initially to study a single atom with a fractional $N=N_{0}+\delta$. Then, the density is piecewise-linear Perdew et al. (1982):

[TABLE]

being a combination of the density of a neutral atom and an anion ($N_{0}$ and $N_{0}+1$ electrons, respectively). When $\delta$ is small and positive, the density has two regions of exponential decay: as $|x|\rightarrow\infty$, $n(x;N)\propto n(x;N_{0}+1)\propto e^{-2\sqrt{2A}|x|}$ (the decay rate is governed by the IP of the anion, which equals the EA of the neutral). We term this the region of ‘$A$-decay’. However, for small $\delta$, $n(x;N_{0})\propto e^{-2\sqrt{2I}|x|}$ starts to dominate the density at some point approaching the nucleus (‘$I$-decay’). Based on Ref. Hodgson et al. (2016), we expect a step in the KS potential to form at the crossover between these two decays.

From the KS perspective, in the $A$-decay region $v_{\mathrm{s}}(x)$ reaches the asymptotic value $v^{\prime}$, and in the $I$-decay region it has the value $v$. In general, $v$ and $v^{\prime}$ may differ, forming the step $S:=v-v^{\prime}$.

In the limit $\delta\rightarrow 0^{+}$, in the $A$-decay region, $n(x;N)\propto|\varphi^{\mathrm{lu}}(x)|^{2}\propto e^{-2\sqrt{2(v^{\prime}-\varepsilon^{\mathrm{lu}})}|x|}$, where $\{\varphi^{i}(x)\}$ are the KS orbitals. In the $I$-decay region, $n(x;N)\propto|\varphi^{\mathrm{ho}}(x)|^{2}\propto e^{-2\sqrt{2(v-\varepsilon^{\mathrm{ho}})}|x|}$. Therefore, $A=v^{\prime}-\varepsilon^{\mathrm{lu}}$ and $I=v-\varepsilon^{\mathrm{ho}}$. Thus, $S=v-v^{\prime}=I+\varepsilon^{\mathrm{ho}}-(A+\varepsilon^{\mathrm{lu}})=\Delta$, exactly as in Eq. (1).

To summarize, for an atom, when $N$ infinitesimally surpasses an integer, the plateau that forms in the region of the atom elevates the level of the KS potential by that required to obtain the exact fundamental gap.

We now model atoms in real space comprised of $N$ same-spin electrons in 1D, where $1\leqslant N\leqslant 2$, using the iDEA code Hodgson et al. (2013)  ††footnotemark: . Solving the Schrödinger equation for $N\geqslant 2$ poses an immense computational challenge, which scales exponentially with $N$. Therefore, by enforcing the same spin for all electrons, we ensure that for a two-electron system two KS orbitals are occupied (in contrast to a spin singlet – one orbital occupied by two electrons with opposite spins), which is necessary for the concepts demonstrated in this Letter to be general 222These concepts do not rely on the electrons forming a singlet, but still apply in such a scenario.. Our electrons interact via the appropriately softened, 1D Coulomb interaction $(\left|x-x^{\prime}\right|+1)^{-1}$ (atomic units) Gordon et al. (2005) ††footnotemark: . The external potential is $v_{\mathrm{ext}}(x)=-2/(0.4\cdot\left|x\right|+1)$; it tends to zero as $|x|\rightarrow\infty$, as appropriate. We initially calculate $n(x)$ for $N=1$ and separately for $N=2$. Then, $n(x)$ for the $(1+\delta)$-electron system is calculated via Eq. (3). Finally, we reverse engineer the exact KS potential from the exact ($1+\delta$)-electron density for varying $\delta$.

Figure 1(a) shows the natural log of the density, $\ln{[n(x)]}$, for $\delta=0,10^{-8},10^{-6}$, and $10^{-4}$, as a function of $x$. Plotting the log helps to recognize the regions of exponential decay. For $\delta=0$, there is only one such region ($I$-decay). For $\delta>0$, we clearly recognize the $I$- and the $A$-decay regions. The above analysis indicates that the points where the decay rate changes are also the points where the steps in $v_{\mathrm{s}}(x)$ occur. This is demonstrated numerically in Fig. 1(b) and (c): for all finite values of $\delta$, the exact KS potential has two spatial steps. The steps act to elevate the level of the KS potential for the central region, where most of the electron density resides [cf. $v_{\mathbf{s}}(x=0)$ for $\delta=0$ and $\delta>0$]. For all small values of $\delta$ the step height is the same and equals $\Delta$, which we calculate from total energy differences and the exact KS eigenvalues; see the Hartree-exchange-correlation (Hxc) potential [$v_{\mathrm{Hxc}}(x)=v_{\mathrm{s}}(x)-v_{\mathrm{ext}}(x)$] in Fig. 1(c). The positions of the steps vary with $\delta$ – the smaller $\delta$ is, the further from the atom the steps form. Therefore, as $\delta\rightarrow 0^{+}$, the plateau becomes a spatially uniform shift in the potential, as required for the exact fundamental gap.

The presented example clearly establishes the relationship between the derivative discontinuity and the steps’ height in the atomic case, namely, that $\Delta=S$. It suggests that approximate xc functionals that are sensitive to the change in the decay rate of the density, and respond by forming steps in the xc potential (see suggestions in Refs. Armiento and Kümmel (2013); Hodgson et al. (2014); Mori-Sánchez and Cohen (2014)), are theoretically capable of producing the expected jump in the KS potential as $N$ surpasses an integer, and therefore yielding an accurate fundamental gap.

We now consider a stretched diatomic molecule, which is one system consisting of two atoms, $L$ and $R$, separated by a large distance, $d$. As $d\rightarrow\infty$, the number of electrons on each atom, $N_{L}$ and $N_{R}$, can be defined. We now imagine transferring an infinitesimal electronic charge, $\delta$, from $L$ to $R$. Therefore, $N_{L}=N^{0}_{L}-\delta$ and $N_{R}=N^{0}_{R}+\delta$, so the total number of electrons, $N_{L\cdots R}=N_{L}+N_{R}$ is constant. Relying on Eq. (1), one may expect a plateau to form in the region of atom $R$, whose height is $\Delta_{R}=I_{R}-A_{R}-(\varepsilon^{\mathrm{lu}}_{R}-\varepsilon^{\mathrm{ho}}_{R})$. Similarly, for a charge transfer from $R$ to $L$, one may expect a plateau $\Delta_{L}$ around $L$. However, Eq. (2) demonstrates that the height of the step in $v_{\mathrm{s}}(\mathbf{r})$ that forms between $L$ and $R$ is independent of the EA of either atom.

To consolidate these seemingly opposing viewpoints, we consider our atomic example above. We realize that for $\delta>0$, we expect the density of atom $R$ to have two regions of exponential decay ($I_{R}$- and $A_{R}$-decay), while the density of atom $L$ will have one such region ($I_{L}$-decay only). For $\delta<0$, the reverse picture applies.

Figure 2(a) shows a diagram of $\ln{[n(x)]}$ very far from, and between, the atoms. In this region the decaying density from atom $R$ changes its decay rate (because of the transferred charge) at point (2), then meets the decaying density from atom $L$ at point (1). Therefore, between the atoms there are two changes in the decay, and hence two steps that manifest in the KS potential (Fig. 2(b)). We emphasize that Step (1) arises because $L\cdots R$ is one system, despite $d$ being large. Notably, the height of the step does not depend on the magnitude of $n(x)$ at the point where it forms, and therefore the step is expected to appear at any $d$, as long as the atoms may be considered one system 333This situation is qualitatively different from the case of a dissociated molecule, originally introduced in Ref. Perdew et al. (1982) and recently discussed in Ref. Kraisler and Kronik (2015), where one considers two completely separated atoms with varying $N_{L}$ and $N_{R}$, at constant $N_{L\cdots R}$. Although $N_{L}$ and $N_{R}$ are varied in a concerted manner, strictly speaking this is no longer one molecule.. The heights of Steps (1) and (2) are derived as before: $S^{(1)}=A_{R}+\varepsilon^{\mathrm{lu}}_{R}-(I_{L}+\varepsilon^{\mathrm{ho}}_{L})$ and $S^{(2)}=I_{R}+\varepsilon^{\mathrm{ho}}_{R}-(A_{R}+\varepsilon^{\mathrm{lu}}_{R})$. They depend both on the IP, the EA, and the ho and lu KS eigenvalues. $S^{(2)}=\Delta_{R}$, whereas $S^{(1)}$ is not so recognizable a quantity. We term this the negative of the ‘charge-transfer derivative discontinuity’, $\Delta_{L\rightarrow R}^{\mathrm{CT}}$, being the difference between the energy it costs to move an electron from $L$ to $R$ and the corresponding difference in the KS energy levels Tozer (2003); Dreuw and Head-Gordon (2004); Hellgren and Gross (2012); Fuks et al. (2013); Fuks and Maitra (2014); Fuks (2016).

The overall difference in the level of $v_{\mathrm{s}}(x)$ in the region of $L$ and the region of $R$, which is the determining characteristic of $v_{\mathrm{s}}(x)$ regarding the distribution of charge between the atoms, is $S=S^{(1)}+S^{(2)}=I_{R}-I_{L}+\varepsilon^{\mathrm{ho}}_{R}-\varepsilon^{\mathrm{ho}}_{L}$, exactly as given by Eq. (2), and independent of $A_{R}$ and $\varepsilon^{\mathrm{lu}}_{R}$, in contrast to the individual gaps $S^{(1)}$ and $S^{(2)}$ 444For $\delta<0$, one obtains $S^{(2)}=-\Delta_{L}$, $S^{(1)}=\Delta_{R\rightarrow L}^{\mathrm{CT}}$, with the overall gap, $S$, remaining the same..

We now model a one-dimensional molecule $L\cdots R$ consisting of two atoms and two same-spin electrons ($N_{L\cdots R}=2$) separated by a large distance $d=40$ a.u. The external potential, $v_{\mathrm{ext}}(x)=-4/(0.8\cdot|x-\tfrac{1}{2}d|+1)-2/(0.4\cdot|x+\tfrac{1}{2}d|+1)$, consists of two wells that represent the left and right nuclei. The potential is chosen so that two interacting electrons occupying this potential localize such that there is one electron’s worth of charge in each well; see Fig. 3(a). To reproduce such a density in the KS system, $v_{\mathrm{s}}(x)$ must form a spatial step between the atoms Hodgson et al. (2016). In the absence of the step, $\varepsilon^{\mathrm{lu}}_{R}$ would be lower than $\varepsilon^{\mathrm{ho}}_{L}$, which would cause both electrons to artificially localize on atom $R$. Figure 3(b) shows a clear step in $v_{\mathrm{s}}(x)$, whose height is given by Eq. (2) and its position is at the point where the decay rate of $n(x)$ changes [cf. Fig. 3(a) and (b)]. To the far right of the atom $R$ there is another change in the decay rate of $n(x)$, which causes a step down (not shown).

For the above molecule we correctly observe no regions of $A$-decay, as none of the excited states of $L$ or $R$ are occupied. We now consider transferring a small amount of charge from $L$ to $R$ by making the right well deeper. This shifts the density minimum to the right, hence shifting the step in $v_{\mathrm{s}}(x)$ correspondingly, but no $A_{R}$-decay region forms. The height of the step is unaffected.

To further investigate the structure of $v_{\mathrm{s}}(x)$ when the density has three regions of decay in between the nuclei, as shown in Fig. 2, we increase $N_{L\cdots R}$ to be $2+\delta$, with $\delta\rightarrow 0^{+}$. We additionally calculate a three-electron system with the same $v_{\mathrm{ext}}(x)$ and use Eq. (3). In the three-electron case, one electron is localized on the left and two on the right. For $\delta=10^{-8}$, the $A_{R}$-decay region can be observed in Fig. 3(c) and two points of decay change can be recognized in $\ln{[n(x)]}$ between the atoms, similar to Fig. 2(a). Correspondingly, two steps in $v_{\mathrm{s}}(x)$ can now be observed between the atoms [Fig. 3(d)]. The sum $S^{(1)}+S^{(2)}$ is described by Eq. (2) and is independent of $A_{R}$.

As $\delta$ decreases, points (1) and (2) travel towards each other until they meet, forming a density minimum. Likewise, steps (1) and (2) coincide and form the overall step $S$. In parallel, far from both atoms the $A_{R}$-decay prevails (not shown). Hence, far to the left and to the right of the molecule we find a step $S_{L\cdots R}=I_{L}-A_{R}-(\varepsilon_{R}^{\mathrm{lu}}-\varepsilon_{L}^{\mathrm{ho}})$ that increases the level of the KS potential everywhere for the molecule; we note $S_{L\cdots R}=-S^{(1)}$. This is exactly the uniform jump predicted by Eq. (1), noting that the global EA and lu level are those of atom $R$, whereas the global IP and ho level are those of atom $L$. Moreover, the DD of atom $R$, $\Delta_{R}$, can be deduced directly from the step structure of the KS potential at $\delta\rightarrow 0^{+}$, simply by adding $S$ and $S_{L\cdots R}$.

This example demonstrates that the sensitivity of approximate xc functionals to the change in the density decay rate and their ability to form steps in the xc potential is crucial to correctly distribute electronic charge in a system with appreciable spatial separation. Suggestions made, e.g., in Refs. Armiento and Kümmel (2013); Hodgson et al. (2014); Mori-Sánchez and Cohen (2014) are relevant in this context.

To conclude, in this Letter we clarified the relationship between the jump experienced by the Kohn-Sham (KS) potential as the number of electrons in the system infinitesimally surpasses an integer, and the spatial steps that form in the KS potential as a result of a change in the density decay rate.

For a single atom, as the electron number passes an integer, two regions of exponential decay in the density manifest far from the nucleus. As a result, two steps in the exact KS potential form a plateau. In the limit of $\delta\rightarrow 0^{+}$, this plateau elevates the level of the potential everywhere in real space by $\Delta$, which is the amount required to obtain the exact fundamental gap from KS eigenvalues.

For a stretched diatomic molecule with an integer $N$, the KS potential forms a step between the atoms to correctly distribute the electronic charge in the molecule. The height of the step is independent of the atomic EAs. The step forms at the interface between the atoms, where the decay rate of the density changes. Slightly increasing the number of electrons above an integer causes the molecule’s KS potential to develop two steps between the atoms. Each step individually depends both on the atomic IP and EA; however, the overall difference in the level of the KS potential between the atoms remains insensitive to the atomic EAs.

The properties of the exact KS potential outlined in this work depend on the fine details of the asymptotic decay of the electron density. Accounting for these properties with standard approximations to the exchange-correlation (xc) functional poses a great challenge. However, since an accurate step structure in an approximate xc potential is crucial to predict the fundamental gap and provide a correct distribution of the electronic charge in complex systems, it should be taken into account in the development of future approximations.