Asymptotic Behavior of Exact Exchange for Slabs: Beyond the Leading Order

Asymptotic Behavior of Exact Exchange for Slabs: Beyond the Leading Order

Computation 2018, 6(2), 35; doi:10.3390/computation6020035

Article

Asymptotic Behavior of Exact Exchange for Slabs: Beyond the Leading Order

Eberhard EngelOrcid

Center for Scientific Computing, J. W. Goethe-Universität Frankfurt, Max-von-Laue-Str. 1, D-60438 Frankfurt/Main, Germany

Received: 26 February 2018 / Accepted: 21 April 2018 / Published: 29 April 2018

Abstract:

Far outside the surface of slabs, the exact exchange (EXX) potential falls off as , if z denotes the direction perpendicular to the surface and the slab is localized around . Similarly, the EXX energy density behaves as , where n is the electron density. Here, an alternative proof of these relations is given, in which the Coulomb singularity in the EXX energy is treated in a particularly careful fashion. This new approach allows the derivation of the next-to-leading order contributions to the asymptotic and . It turns out that in both cases, the corrections are proportional to in general.

Keywords:

density functional theory; exact exchange; slabs; asymptotic behavior

PACS:

71.15.Mb; 71.15.-m

1. Introduction and Summary of Results

The application of Kohn–Sham density functional theory (DFT) to surfaces was initiated by the seminal work of Lang and Kohn [1,2,3] (see also [4]). While the initial focus was on metals and the jellium model, soon, also semi-conducting materials were investigated [5]. In practice, DFT calculations for surfaces often rely on slabs, in which a limited number of atomic layers simulates the bulk material [5,6]. The number of layers required for an accurate description of the surface primarily depends on the electronic structure, in particular on the localization of the states. Obviously, a decoupling of the surface states on the two sides of the slab is desirable. If relaxation of the atomic positions at the surface or surface reconstruction is of interest, one has to make sure that some bulk-like layers remain in the middle of the slab. Nevertheless, often a rather small number of layers is sufficient in the case of non-metallic solids. Particularly accurate results for quantities such as the work function and the surface energy can be obtained by extrapolating data from slabs with different thicknesses to the limit of infinite thickness.

However, slabs are also of interest in themselves. The most prominent example is single- and bi-layer graphene [7,8,9,10], but also other 2D materials have attracted considerable attention recently, such as silicene [11,12,13], hexagonal boron nitride nanosheets [14,15] and mono- and multi-layer transition metal dichalcogenides [16,17,18] (as well as structures obtained by a combination of these materials with graphene).

The most important quantity for the density functional description of surfaces and slabs is the exchange-correlation (xc) potential (for an overview of DFT, see [19]). Starting with Lang and Kohn [1,2,3], the xc-potential at surfaces has been studied extensively [20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43]. While the discussion of the xc-potential turned out to be difficult for semi-infinite matter (see [31,33]), definitive information is available for the exact exchange (EXX) potential of slabs. Both and the EXX energy density of metallic slabs have been investigated by Horowitz and collaborators [28,29,32,36] and by Qian [33] (based on extensive work by Sahni and collaborators [23,24,25,26] on semi-infinite matter), relying on the jellium model (i.e., on full translational symmetry in the directions parallel to the surface). If the slab is parallel to the -plane and has its center at , the exact behaves as far outside the surface (atomic units are used throughout this work). Similarly, the EXX energy density of jellium slabs falls off as , where n is the electron density [28]. This latter result is of particular interest, since defines the Slater potential, , which is the core contribution to the often used Krieger–Li–Iafrate approximation [44] for the exact and its refinement, the localized Hartree–Fock/common energy denominator approximation [45,46]. Both results have subsequently been generalized [34,35] to non-jellium slabs, for which one has a Bravais lattice in the -directions,

(1)

(2)

Here, L characterizes the position of the surface, i.e., the slab extends from to L.

In this contribution, a particularly careful derivation of the relations (1) and (2) is given, in which the crucial Coulomb singularity in the EXX energy is handled in a mathematically unambiguous fashion. In contrast to the approaches used in [34,35], this derivation allows one to extract the next-to-leading order contributions to both and . Moreover, unlike the argument in [35], the present proof of (2) does not employ the ultimate asymptotic form of the density from the very outset, but rather is based on a completely general variation of the asymptotic n as a function of .

The derivation is prepared by a brief review of the asymptotic form of the Kohn–Sham (KS) states (in Section 2), relying on the analysis of [47]. In particular, the coupling of the Fourier amplitudes for the states is discussed, in order to show which amplitudes are asymptotically dominant. On this basis, the asymptotic behavior of (Section 3), (Section 4) and (Section 5) is analyzed. In this analysis, the Coulomb singularity is kept under control by addition and subtraction of a suitably-screened exchange kernel. One finds that the next-to-leading order term in has the form where p is a generalized density with weights determined by the first moments of the transition matrix. p can be explicitly expressed in terms of n, if the asymptotic density is dominated by the states in the vicinity of a single -point in the interior of the first Brillouin zone (which is the standard situation for non-metallic slabs [43]). In this case, one finds , where denotes the first moment of the most weakly-decaying state at . In the same fashion, the next-to-leading order contribution to is obtained as , and this correction also applies to . It seems worthwhile to remark that all results can be generalized to the situation that there are several degenerate most weakly-decaying states at . Moreover, the derivation of higher order corrections can basically follow the same route.

2. Asymptotic Behavior of States

In the case of slabs, the total KS potential is invariant under translation by an arbitrary two-dimensional (2D) Bravais vector in the -directions, but confines the electrons in the z-direction,

(3)

Here, , and is a vector of the 2D reciprocal lattice in the -directions. The corresponding KS states are given by:

(4)

where is the 2D crystal momentum. The normalization is chosen so that integrates up to one for a single 2D unit cell with area A,

(5)

The coefficients satisfy the KS equations on the reciprocal lattice,

(6)

In [35], it has been shown that the exact behaves as for . In the following, is therefore assumed to have the asymptotic form:

(7)

with u being allowed to vanish. For the Fourier component with , one then has:

(8)

while all other decay at least as fast as (typically even faster). It seems worthwhile to emphasize that the next-to-leading order contributions to not included in Equations (7) and (8) are irrelevant for the derivation of the asymptotic forms of , and in Section 3, Section 4 and Section 5.

The asymptotic behavior of the solutions of the differential equation (6) with the potential (8) has been analyzed in detail in [47]. For large z, the solutions can be expressed as:

(9)

(the quantum numbers are omitted for brevity in the rest of this section). Here, denotes the general solution of the homogeneous differential equation:

(10)

For the asymptotic potential (8), this solution is to leading order given by:

(11)

with:

(12)

On the other hand, accounts for the coupling of all amplitudes in (6),

(13)

( must be chosen sufficiently large, so that the -integration only extends over the asymptotic region).

The asymptotic behavior of the solutions is primarily controlled by the exponent , Equation (12). As long as for all , i.e., for states inside the first Brillouin zone (BZ), there is a unique most slowly-decaying amplitude, . For states on the boundary of the first BZ, on the other hand, one has a second amplitude with , which shows the same asymptotic decay as . In fact, for large z, these two amplitudes only differ by the constant in Equation (11), since .

For all states and all , the function falls off faster than . The contribution of is largest for states on the boundary of the first BZ, for which a second amplitude decays as slowly as . For these states, one finds that is suppressed by:

compared to (and an analogous statement holds for the relation between and ). However, vanishes at least as fast as , so that the second term on the right-hand side of (9) is irrelevant for the asymptotically-leading amplitude(s) in the case of all states. One thus has:

(14)

(and an analogous relation for for states on the boundary of the first BZ).

In the case of all other , the contribution of cannot be neglected, since, at least in general, decays faster than the product . However, all decay faster than the most weakly-decaying amplitude , so that Equation (9) can be solved iteratively for large z. A single iteration of this equation, i.e., replacement of the full solution on the right-hand side of (13) by , is sufficient to obtain the correct asymptotic behavior of the next-to-most weakly-decaying amplitudes. For all yet faster decaying amplitudes, further iteration can (but need not) be required, depending on the asymptotic behavior of the involved.

In summary, the asymptotic behavior the states well inside the first BZ is determined by the amplitude with , which, in turn, approaches the homogeneous solution (11) for large z. All amplitudes with are suppressed relative to , with the suppression factor depending on the structure of the potential. In the case of the next-to-most weakly-decaying amplitudes, the suppression factor is either given by the most weakly-decaying Fourier component of the potential with non-zero or by , if this ratio decays more slowly than .

In the following, it is assumed for simplicity that: (i) the most weakly-decaying occupied state in the first BZ, i.e., the minimum exponent , is found for a single -point well inside the first BZ; (ii) is Taylor expandable at this -point; and (iii) all states in a complete neighborhood of this -point are also occupied. In this situation, the statements of the previous paragraph apply to all asymptotically-relevant states, and the simplest version of the BZ integration scheme of [43] can be used.

3. Asymptotic Behavior of Exchange Energy Density

In the case of spin-saturated slabs, and are given by:

(15)

(16)

where denotes the occupation of the state and the -integrations extend over the first BZ. Insertion of (4) into the exchange energy density (15) and use of the 2D Fourier transform of the Coulomb interaction yields:

(17)

(18)

In order to extend the -integration to the complete 2D k-space, one uses the fact that the coefficients only depend on ,

(19)

If one extends the occupation function accordingly,

(20)

the exchange energy density can be expressed as:

(21)

At this point, one can choose the origin of the -integration to be at ,

(22)

with:

(23)

In the following, the behavior of in the regime is analyzed (with L characterizing the extension of the slab whose center is at the origin). In order to isolate the asymptotically-dominating term without any mathematical ambiguity, the integrand is first split into a short- and a long-wavelength component,

(24)

(25)

(26)

The -integral in (25) can be performed directly,

(27)

For large z, the overlap function vanishes exponentially with z. This is immediately clear for the integrand in (23), since, in view of Equation (14), has the form (11) for large z, and all decay even faster than . Moreover, this exponential decay is preserved by the -integration, which may be performed using the approach of [43]. Consequently, the second term in (27) is exponentially smaller than the first term in the asymptotic regime,

(28)

The function also vanishes exponentially for , the arguments being the same as in the case of z. In fact, the BZ integration scheme of [43] is particularly appropriate if both z and are large. The exponential decay of with large allows one to expand the denominator in powers of . Subsequently, one can extend the -integration to infinity, without introducing additional power-law corrections,

(29)

The leading term of this expansion can be evaluated by use of (5),

This expression is easily identified as the Fourier component of the density (16),

(30)

The long-range component of the exchange energy density therefore asymptotically behaves as:

(31)

with:

(32)

Since can be chosen reasonably small, one obtains for the leading two orders:

(33)

Now, consider the short-wavelength component of , Equation (26). In order to analyze its behavior for large z, the -integration is decomposed into integrals over an elementary (simply connected) cell V, which is periodically repeated to fill the complete 2D momentum space without leaving any void. V could, for instance, have the shape of the first BZ, but a reduced spatial extension. In this case, the vectors , which connect the centers of the cells, are correspondingly scaled reciprocal lattice vectors. However, for the present purpose, V could also be a small square, so that the vectors are the reciprocal lattice vectors of the associated 2D simple cubic lattice,

As soon as , there is no Coulomb singularity in the -integral: in this case, the -integral is well defined for arbitrary . For , these contributions decay exponentially faster than the term with , so that only the latter term remains to be analyzed,

(34)

For large z, an expansion of the expression in square brackets about is legitimate, since: (i) the -integral is dominated by the vicinity of the Coulomb singularity; (ii) the integrand falls off rapidly for large z; and (iii) V can be chosen much smaller than the first BZ,

(35)

Due to the inversion symmetry of V, one has:

(36)

The -integral in the second term on the right-hand side of (35) can be evaluated further in polar coordinates,

(37)

where defines the boundary of the integration area V. If, for instance, V is chosen to be a square with its corners at , is given by:

The leading contribution to is thus obtained as

(38)

It seems worthwhile to emphasize that the -integral is well defined, since for close to z, one has:

Taking this into account, the leading term for large z is given by:

(39)

since for all and decays exponentially with . Upon addition of (33) and (39), one finds for the leading two orders of the asymptotic energy density,

(40)

The result is independent of , as it should be. In real space, one thus obtains:

(41)

4. Asymptotic Behavior of Slater Potential

Equation (41) shows that the Slater potential asymptotically behaves as:

(42)

In order to evaluate the next-to-leading order term further, the -integrations involved in and need to be analyzed. Asymptotically, the density (30) is dominated by the coefficients with the lowest exponents for which . If the lowest value of is found well inside the first BZ (which is usually the case [43]), the coefficients decay most slowly, as discussed in Section 2,

(43)

with given by (11), due to (14). The same argument can be applied to the asymptotic ,

(44)

In general, one particular band falls off most slowly for each . On the other hand, if there are two degenerate highest occupied states and for some , the asymptotically-leading terms in the corresponding coefficients and only differ by some constant, since . It is thus sufficient to restrict the expression (44) to the terms with ,

(45)

Here, denotes the first moment of ,

(46)

in the non-degenerate case and a combination of the moments of the degenerate states otherwise.

The combination of Equations (43) and (45) with the asymptotic form (11) of shows that the asymptotic behavior of both the density (30) and the first moment (32) is controlled by a BZ integral of the form discussed in [43], i.e., Equation (18) of this reference. In the case of the density, the generic integrand in Equation (18) of [43] corresponds to , for the first moment corresponds to . Both (43) and (45) are asymptotically dominated by the state(s) with the minimum exponent , as discussed in [43]. If there is only a single -point and band for which reaches its minimum value inside the integration region, one finds:

(47)

(48)

where is defined by:

(49)

therefore has the same asymptotic behavior as , so that . The ratio also approaches a constant for large z, if , since this affects the z-dependence of p and n in the same way. Moreover, the same statement applies, if there is more than one -point and/or more than one band for which assumes its minimum value (see [43]). One thus finally arrives at:

(50)

with being replaced by a superposition of the moments of all asymptotically-relevant states in the more general situation. Depending on the symmetry of the relevant states at under the transformation , can vanish. This is the case in particular, if one has reflection symmetry with respect to the -plane, so that .

5. Asymptotic Behavior of Exact Exchange Potential

For a discussion of the exact , one starts with the total exchange energy per unit cell, expressed in terms of the amplitudes ,

(51)

(which is obtained by integration of (17) over the complete unit cell). In the following, the variation of this expression under a variation of the asymptotic density is studied. However, a variation of the asymptotic density implies a variation of , since there exists a one-to-one correspondence between the KS states and the density (30),

(52)

In addition, since only the asymptotic density is to be varied (but not the complete ), is non-zero only for large z [ can be assumed to have the standard shape of a test function].

Now, consider the variation of under the variation of ,

(53)

Using (19) and (20), the -integration can be combined with the sum over . If the origin of the resulting -integration over the complete k-space is chosen to be at , one obtains:

(54)

Next, is split into a long- and a short-wavelength component, in analogy to (24),

(55)

(56)

The evaluation of proceeds in straightforward fashion. One first performs the -integral,

(57)

One can then follow the argument leading from (27) to (29), since falls off exponentially with z, irrespective of ,

(58)

First, focus on the leading order term of the expansion (L0). This term can be evaluated further by use of the orthonormality relation (5),

(59)

One can now re-introduce the x-y-integrations and identify ,

(60)

Finally, an expansion in powers of is legitimate, since can be chosen to be much smaller than the z-values for which is non-zero,

(61)

Next, consider the short-wavelength component (56),

(62)

The arguments that lead from (26) to (35) also apply to (62). Using (36), one arrives at:

(63)

One can then employ (37) for the -integration and expand in powers of , since the remaining -integral is free of any singularities,

(64)

Use of the orthonormality relation (5) and the density variation (52) then shows that this expression cancels exactly with the first order contribution to (61).

Finally, consider the first order contribution (L1) to the expansion in (58),

(65)

A completely general representation of in terms of , as achieved for , is not obvious. In the following, therefore, only the most important class of variations will be considered. For very large z, the BZ integral in the density (30) is dominated by the vicinity of the -point(s) for which the exponent , Equation (12), assumes its lowest value in the integration region [43]. If one is only interested in a variation of the asymptotic density, it is therefore sufficient to vary only the amplitudes of the states in the vicinity of the states with minimal . Assume, for simplicity, that the lowest value of the exponent is obtained only for a single state, denoted by . If only the amplitudes in the vicinity of are varied, while all other are kept unchanged, is given by:

(66)

Now, consider for this variation. For any given , i.e., any given state , the sum over in (65) is dominated by the most weakly-decaying occupied amplitude among the . In the case of the variation (66), thus reduces to:

(67)

with defined by Equation (46). Since is non-zero only for large z and decays exponentially for these z-values, the -integration in (67) can be simplified by a Taylor expansion of the smooth parts of the integrand about ,

(68)

At this point, one can finally re-introduce the -integrations and utilize (66),

(69)

Thus, for the relevant class of variations can be explicitly expressed in terms of . A completely general variation of n includes this particular class.

One thus concludes that, in general, the next-to-leading order contribution to the asymptotic exact exchange potential is proportional to ,

(70)

where (61), (64) and (69) have finally been combined.

6. Concluding Remarks

In this work, the next-to-leading order contributions to the EXX energy density, the associated Slater potential and the exact EXX potential have been calculated for arbitrary slabs, including semi-conducting and insulating materials. If the asymptotic density is dominated by the states of the band in the vicinity of a single -point in the interior of the first BZ, one finds:

(71)

(72)

(73)

where denotes the first moment (46) of the most weakly-decaying occupied state. If there are several degenerate most weakly-decaying states at the -point or states at several -points show the same asymptotic behavior (due to symmetry), this coefficient has to be suitably generalized (compare [43]).

The results (71)–(73) show exactly the required behavior under a translation of the coordinate system. Consider two coordinate systems whose origin differs by a shift in the z-direction, . At the same physical point far from the surface, Equation (73) then yields in the two coordinate systems,

(74)

where denotes the moment of the state evaluated with respect to the primed coordinate system,

(75)

Insertion into (74) demonstrates that this relation is in fact satisfied in the leading two orders,

Equation (75) reflects the fact that the value of a moment always depends on the coordinate system in which it is calculated. As long as the origin of the coordinate system is chosen to be in the middle of the slab, its extension L essentially provides an upper limit for the size of ; a highly localized most weakly-decaying state right at one of the two surfaces of the slab would yield . However, as Equation (75) shows, can also be made to vanish in a suitably chosen coordinate system. This approach could be particularly attractive for the prime application of Equations (72) and (73), the normalization of numerically-generated EXX potentials.

Funding

This research received no external funding.

Conflicts of Interest

The author declares no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

EXX | exact exchange
xc | exchange-correlation
BZ | Brillouin zone
DFT | density functional theory
KS | Kohn–Sham
2D | two-dimensional

References

1. Lang, N.D.; Kohn, W. Theory of metal surfaces: Charge density and surface energy. Phys. Rev. B 1970, 1, 4555. [Google Scholar] [CrossRef]

2. Lang, N.D.; Kohn, W. Theory of metal surfaces: Work function. Phys. Rev. B 1971, 3, 1215. [Google Scholar] [CrossRef]

3. Lang, N.D.; Kohn, W. Theory of metal surfaces: Induced surface charge and image potential. Phys. Rev. B 1973, 7, 3541. [Google Scholar] [CrossRef]

4. Smith, J.R. Self-consistent many-electron theory of electron work functions and surface potential characteristics for selected metals. Phys. Rev. 1969, 181, 522. [Google Scholar] [CrossRef]

5. Schlüter, M.; Chelikowsky, J.R.; Louie, S.G.; Cohen, M.L. Self-Consistent Pseudopotential Calculations on Si(111) Unreconstructed and (2 × 1) Reconstructed Surfaces. Phys. Rev. Lett. 1975, 34, 1385. [Google Scholar] [CrossRef]

6. Chelikowsky, J.R.; Schlüter, M.; Louie, S.G.; Cohen, M.L. Self-Consistent Pseudopotential Calculation for the (111) Surface of Aluminum. Solid State Commun. 1975, 17, 1103. [Google Scholar] [CrossRef]

7. Castro Neto, A.H.; Guinea, F.; Peres, N.M.R.; Novoselov, K.S.; Geim, A.K. The electronic properties of graphene. Reviews of modern physics. Rev. Mod. Phys. 2009, 81, 109. [Google Scholar] [CrossRef]

8. Morozov, S.V.; Novoselov, K.; Katsnelson, M.; Schedin, F.; Elias, D.; Jaszczak, J.; Geim, A. Giant intrinsic carrier mobilities in graphene and its bilayer. Phys. Rev. Lett. 2008, 100, 016602. [Google Scholar] [CrossRef] [PubMed]

9. Zhang, Y.; Tang, T.; Girit, C.; Hao, Z.; Martin, M.C.; Zettl, A.; Crommie, M.F.; Shen, Y.R.; Wang, F. Direct observation of a widely tunable bandgap in bilayer graphene. Nature 2009, 459, 820. [Google Scholar] [CrossRef] [PubMed]

10. Min, L.; Hovden, R.; Huang, P.; Wojcik, M.; Muller, D.A.; Park, J. Twinning and Twisting of Tri- and Bilayer Graphene. Nano Lett. 2012, 12, 1609–1615. [Google Scholar]

11. Cahangirov, S.; Topsakal, M.; Aktürk, E.; Sahin, H.; Ciraci, S. Two- and one-dimensional honeycomb structures of silicon and germanium. Phys. Rev. Lett. 2009, 102, 236804. [Google Scholar] [CrossRef] [PubMed]

12. Aufray, B.; Kara, A.; Vizzini, S.B.; Oughaddou, H.; LéAndri, C.; Ealet, B.; Le Lay, G. Graphene-like silicon nanoribbons on Ag (110): A possible formation of silicene. Appl. Phys. Lett. 2010, 96, 183102. [Google Scholar] [CrossRef]

13. Lalmi, B.; Oughaddou, H.; Enriquez, H.; Kara, A.; Vizzini, S.B.; Ealet, B.N.; Aufray, B. Epitaxial growth of a silicene sheet. Appl. Phys. Lett. 2010, 97, 223109. [Google Scholar] [CrossRef]

14. Li, L.H.; Chen, Y. Atomically thin boron nitride: Unique properties and applications. Adv. Funct. Mater. 2016, 26, 2594–2608. [Google Scholar] [CrossRef]

15. Bhimanapati, G.R.; Glavin, N.R.; Robinson, J.A. 2D Materials. In Semiconductors and Semimetals; Iacopi, F., Boeckl, J.J., Jagadish, C., Eds.; Elsevier: Amsterdam, The Netherlands, 2016; Volume 95, p. 101. [Google Scholar]

16. Mak, K.F.; Lee, C.; Hone, J.; Shan, J.; Heinz, T.F. Atomically thin MoS₂: A new direct-gap semiconductor. Phys. Rev. Lett. 2010, 105, 136805. [Google Scholar] [CrossRef] [PubMed]

17. Radisavljevic, B.; Radenovic, A.; Brivio, J.; Giacometti, V.; Kis, A. Single-layer MoS₂ transistors. Nat. Nanotechnol. 2011, 6, 147. [Google Scholar] [CrossRef] [PubMed]

18. Yu, W.J.; Liu, Y.; Zhou, H.; Yin, A.; Li, Z.; Huang, Y.; Duan, X. Highly efficient gate-tunable photocurrent generation in vertical heterostructures of layered materials. Nat. Nanotechnol. 2013, 8, 952. [Google Scholar] [CrossRef] [PubMed]

19. Engel, E.; Dreizler, R.M. Density Functional Theory: An Advanced Course; Springer: Berlin, Germany, 2011. [Google Scholar]

20. Sham, L.J. Density-functional theory of the band gap. Phys. Rev. B 1985, 32, 3876. [Google Scholar] [CrossRef]

21. Eguiluz, A.G.; Hanke, W. Evaluation of the exchange-correlation potential at a metal surface from many-body perturbation theory. Phys. Rev. B 1989, 39, 10433. [Google Scholar] [CrossRef]

22. Eguiluz, A.G.; Heinrichsmeier, M.; Fleszar, A.; Hanke, W. First-principles evaluation of the surface barrier for a Kohn–Sham electron at a metal surface. Phys. Rev. Lett. 1992, 68, 1359. [Google Scholar] [CrossRef] [PubMed]

23. Solomatin, A.; Sahni, V. Analytical asymptotic structure of the exchange and correlation potentials at a metal surface. Phys. Lett. A 1996, 212, 263–269. [Google Scholar] [CrossRef]

24. Solomatin, A.; Sahni, V. Analytical asymptotic structure of the Kohn–Sham exchange potential at a metal surface. Phys. Rev. B 1997, 56, 3655. [Google Scholar] [CrossRef]

25. Qian, Z.; Sahni, V. Quantum mechanical image potential theory. Phys. Rev. B 2002, 66, 205103. [Google Scholar] [CrossRef]

26. Qian, Z.; Sahni, V. Exact electronic properties in the classically forbidden region of a metal surface. Int. J. Quantum Chem. 2005, 104, 929–945. [Google Scholar] [CrossRef]

27. Horowitz, C.M.; Proetto, C.R.; Rigamonti, S. Kohn–Sham exchange potential for a metallic surface. Phys. Rev. Lett. 2006, 97, 026802. [Google Scholar] [CrossRef] [PubMed]

28. Horowitz, C.M.; Proetto, C.R.; Pitarke, J.M. Exact-exchange Kohn–Sham potential, surface energy, and work function of jellium slabs. Phys. Rev. B 2008, 78, 085126. [Google Scholar] [CrossRef]

29. Horowitz, C.M.; Constantin, L.A.; Proetto, C.R.; Pitarke, J.M. Position-dependent exact-exchange energy for slabs and semi-infinite jellium. Phys. Rev. B 2009, 80, 235101. [Google Scholar] [CrossRef]

30. Horowitz, C.M.; Proetto, C.R.; Pitarke, J.M. Localized versus extended systems in density functional theory: Some lessons from the Kohn–Sham exact exchange potential. Phys. Rev. B 2010, 81, 121106. [Google Scholar] [CrossRef]

31. Constantin, L.A.; Pitarke, J.M. Adiabatic-connection-fluctuation-dissipation approach to long-range behavior of exchange-correlation energy at metal surfaces: A numerical study for jellium slabs. Phys. Rev. B 2011, 83, 075116. [Google Scholar] [CrossRef]

32. Luo, H.; Horowitz, C.M.; Flad, H.J.; Proetto, C.R.; Hackbusch, W. Direct comparison of optimized effective potential and Hartree-Fock self-consistent calculations for jellium slabs. Phys. Rev. B 2012, 85, 165133. [Google Scholar] [CrossRef]

33. Qian, Z. Asymptotic behavior of the Kohn–Sham exchange potential at a metal surface. Phys. Rev. B 2012, 85, 115124. [Google Scholar] [CrossRef]

34. Engel, E. Exact exchange plane-wave-pseudopotential calculations for slabs. J. Chem. Phys. 2014, 140, 18A505. [Google Scholar] [CrossRef] [PubMed]

35. Engel, E. Asymptotic behavior of exact exchange potential of slabs. Phys. Rev. B 2014, 89, 245105. [Google Scholar] [CrossRef]

36. Rigamonti, S.; Horowitz, C.M.; Proetto, C.R. Spin-dependent optimized effective potential formalism for open and closed systems. Phys. Rev. B 2015, 92, 235145. [Google Scholar] [CrossRef]

37. Lazar, P.; Otyepka, M. Accurate surface energies from first principles. Phys. Rev. B 2015, 91, 115402. [Google Scholar] [CrossRef]

38. Ye, L.H. Surface calculations with asymptotically long-ranged potentials in the full-potential linearized augmented plane-wave method. Phys. Rev. B 2015, 92, 115132. [Google Scholar] [CrossRef]

39. Ye, L.H. Accurate ionization potential of semiconductors from efficient density functional calculations. Phys. Rev. B 2016, 94, 035113. [Google Scholar] [CrossRef]

40. De Waele, S.; Lejaeghere, K.; Sluydts, M.; Cottenier, S. Error estimates for density-functional theory predictions of surface energy and work function. Phys. Rev. B 2016, 94, 235418. [Google Scholar] [CrossRef]

41. Constantin, L.A.; Fabiano, E.; Pitarke, J.M.; Della Sala, F. Semilocal density functional theory with correct surface asymptotics. Phys. Rev. B 2016, 93, 115127. [Google Scholar] [CrossRef]

42. Ruzsinszky, A.; Constantin, L.A.; Pitarke, J.M. Kernel-corrected random-phase approximation for the uniform electron gas and jellium surface energy. Phys. Rev. B 2016, 94, 165155. [Google Scholar] [CrossRef]

43. Engel, E. Exact exchange potential for slabs: Asymptotic behavior of the Krieger-Li-Iafrate approximation. Phys. Rev. B 2018, 97, 075102. [Google Scholar] [CrossRef]

44. Krieger, J.B.; Li, Y.; Iafrate, G.J. Derivation and application of an accurate Kohn–Sham potential with integer discontinuity. Phys. Lett. A 1990, 146, 256–260. [Google Scholar] [CrossRef]

45. Gritsenko, O.V.; Baerends, E.J. Orbital structure of the Kohn–Sham exchange potential and exchange kernel and the field-counteracting potential for molecules in an electric field. Phys. Rev. A 2001, 64, 042506. [Google Scholar] [CrossRef]

46. Della Sala, F.; Görling, A. Efficient localized Hartree–Fock methods as effective exact-exchange Kohn–Sham methods for molecules. J. Chem. Phys. 2001, 115, 5718–5732. [Google Scholar] [CrossRef]

47. Engel, E. Exact exchange plane-wave-pseudopotential calculations for slabs: Extending the width of the vacuum. Phys. Rev. B 2018, 97, 155112. [Google Scholar] [CrossRef]

© 2018 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).