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

Abstract

Far outside the surface of slabs, the exact exchange (EXX) potentialvxfalls off as−1/z, if z denotes the direction perpendicular to the surface and the slab is localized aroundz=0. Similarly, the EXX energy densityexbehaves as−n/(2z), 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 asymptoticvxandex. It turns out that in both cases, the corrections are proportional to1/z2in general.

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 _{vxc vxc} (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 _{vx vx} of slabs. Both _{vx vx} and the EXX energy density _{ex ex} 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 xy xy -plane and has its center at z=0 z=0 , the exact _{vx vx} behaves as −1/z −1/z far outside the surface (atomic units are used throughout this work). Similarly, the EXX energy density of jellium slabs falls off as −n/(2z) −n/(2z) , where n is the electron density [28]. This latter result is of particular interest, since _{ex ex} defines the Slater potential, _vSlater=2_ex/n _vSlater=2_ex/n , which is the core contribution to the often used Krieger–Li–Iafrate approximation [44] for the exact _{vx vx} 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 xy xy-directions,

_ex(r)→z≫L−n(r)2z

_vx(r)→z≫L−1z.

Here, L characterizes the position of the surface, i.e., the slab extends from −L −Lto 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 _{ex ex} and _{vx vx} . 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 r r.

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 _{ex ex} (Section 3), _{vSlater vSlater} (Section 4) and _{vx vx} (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 _{ex ex} has the form −p(r)/(2^z2) −p(r)/(2^z2) 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 k k -point q q in the interior of the first Brillouin zone (which is the standard situation for non-metallic slabs [43]). In this case, one finds p(r)→z→∞_mqβn(r) p(r)→z→∞_mqβn(r) , where _{mqβ mqβ} denotes the first moment of the most weakly-decaying state β β at q q . In the same fashion, the next-to-leading order contribution to _vx(r) _vx(r) is obtained as −_mqβ/^z2 −_mqβ/^z2 , and this correction also applies to _{vSlater vSlater} . It seems worthwhile to remark that all results can be generalized to the situation that there are several degenerate most weakly-decaying states at q q. 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 _{vs vs} is invariant under translation by an arbitrary two-dimensional (2D) Bravais vector in the xy xy-directions, but confines the electrons in the z-direction,

_vs(r)=∑G^eiG·_r∥v(G,z).

Here, _r∥=(x,y) _r∥=(x,y) , and G G is a vector of the 2D reciprocal lattice in the xy xy-directions. The corresponding KS states are given by:

ϕkα(r)=eik·r∥A∑GeiG·r∥ckα(G,z),

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

_δα,^α′=_∫−∞∞dz∑G_ckα*(G,z)_ck^α′(G,z).

The coefficients _{ckα ckα}satisfy the KS equations on the reciprocal lattice,

^(G+k)2−_∂z22_ckα(G,z)+∑^G′v(G−^G′,z)_ckα(^G′,z)=_ϵkαckα(G,z).

In [35], it has been shown that the exact _{vx vx} behaves as −1/z −1/z for z→∞ z→∞ . In the following, _vs(r) _vs(r)is therefore assumed to have the asymptotic form:

_vs(r)→z≫L−uz+…,

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

v(0,z)→z≫L−uz+…,

while all other v(G,z) v(G,z) decay at least as fast as 1/^z2 1/^z2 (typically even faster). It seems worthwhile to emphasize that the next-to-leading order contributions to v(0,z) v(0,z) not included in Equations (7) and (8) are irrelevant for the derivation of the asymptotic forms of _{ex ex} , _{vSlater vSlater} and _{vx vx} 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:

c(G,z)=_ch(G,z)+∑^G′≠G_ci(G,^G′,z)

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

_{∂z2 ch}(G,z)=^(G+k)2−2ϵ+2v(0,z)_ch(G,z).

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

_ch(G,z)→z≫L_f0(G)^{zu/γ(G)e−γ(G)z},

with:

γ(G)=^{(G+k)2−2ϵ1/2}.shi

On the other hand, _{ci ci} accounts for the coupling of all amplitudes in (6),

_ci(G,^G′,z)=−1γ(G)_∫z1zd^{z′zz′u/γ(G) eγ(G)(z′−z)}v(G−^G′,^z′)c(^G′,^z′)−1γ(G)_∫z∞d^{z′z′zu/γ(G) eγ(G)(z−z′)}v(G−^G′,^z′)c(^G′,^z′)

( _z1<z _z1<z must be chosen sufficiently large, so that the ^{z′ z′}-integration only extends over the asymptotic region).

The asymptotic behavior of the solutions is primarily controlled by the exponent γ(G) γ(G) , Equation (12). As long as ^(G+k)2>^{k2 (G+k)2}>^k2 for all G≠0 G≠0 , i.e., for states inside the first Brillouin zone (BZ), there is a unique most slowly-decaying amplitude, _ch(0,z) _ch(0,z) . For states on the boundary of the first BZ, on the other hand, one has a second amplitude _ch(G¯,z) _ch(G¯,z) with ^(G¯+k)2=^{k2 (G¯+k)2}=^k2 , which shows the same asymptotic decay as _ch(0,z) _ch(0,z) . In fact, for large z, these two amplitudes only differ by the constant _{f0 f0} in Equation (11), since γ(G¯)=γ(0) γ(G¯)=γ(0).

For all states and all ^G′≠0 ^G′≠0 , the function _ci(0,^G′,z) _ci(0,^G′,z) falls off faster than _ch(0,z) _ch(0,z) . The contribution of _ci(0,^G′,z) _ci(0,^G′,z) is largest for states on the boundary of the first BZ, for which a second amplitude _ch(G¯,z) _ch(G¯,z) decays as slowly as _ch(0,z) _ch(0,z) . For these states, one finds that _ci(0,G¯,z) _ci(0,G¯,z)is suppressed by:

_∫z∞d^z′v(−G¯,^z′)

compared to _ch(0,z) _ch(0,z) (and an analogous statement holds for the relation between _ci(G¯,0,z) _ci(G¯,0,z) and _ch(G¯,z) _ch(G¯,z) ). However, _∫z∞d^z′v(−G¯,^z′) _∫z∞d^z′v(−G¯,^z′) vanishes at least as fast as 1/z 1/z , 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:

c(0,z)→z≫L_ch(0,z)

(and an analogous relation for c(G¯,z) c(G¯,z)for states on the boundary of the first BZ).

In the case of all other c(G,z) c(G,z) , the contribution of _ci(G,^G′,z) _ci(G,^G′,z) cannot be neglected, since, at least in general, _ch(G,z) _ch(G,z) decays faster than the product v(G,z)_ch(0,z) v(G,z)_ch(0,z) . However, all _ci(G,^G′,z) _ci(G,^G′,z) decay faster than the most weakly-decaying amplitude _ch(0,z) _ch(0,z) , so that Equation (9) can be solved iteratively for large z. A single iteration of this equation, i.e., replacement of the full solution c(^G′,z) c(^G′,z) on the right-hand side of (13) by _ch(^G′,z) _ch(^G′,z) , 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 v(G,z) v(G,z)involved.

In summary, the asymptotic behavior the states well inside the first BZ is determined by the amplitude with G=0 G=0 , which, in turn, approaches the homogeneous solution (11) for large z. All amplitudes with G≠0 G≠0 are suppressed relative to _ch(0,z) _ch(0,z) , with the suppression factor c(G,z)/_ch(0,z) c(G,z)/_ch(0,z) 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 v(_G1,z) v(_G1,z) with non-zero _{G1 G1} or by _ch(G,z)/_ch(0,z) _ch(G,z)/_ch(0,z) , if this ratio decays more slowly than v(_G1,z) v(_G1,z).

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 γ(0) γ(0) , is found for a single k k -point well inside the first BZ; (ii) γ(0) γ(0) is Taylor expandable at this k k -point; and (iii) all states in a complete neighborhood of this k k -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, _ex(r) _ex(r) and n(r) n(r)are given by:

_ex(r)=−^A2 _∫1BZ^d2k^{(2π)2d2 k′(2π)2}∑α,^α′_{ΘkαΘ^{k′ α′}}∫^{d3 r′}_ϕkα†(r)_{ϕ^{k′ α′}}(r)_{ϕ^{k′ α′}†}(^r′)_ϕkα(^r′)|r−^r′|

n(r)=2A_∫1BZ^d2k^(2π)2∑α_Θkα ^|_ϕkα(r)|2,shift

where _{Θkα Θkα} denotes the occupation of the state kα kα and the k k -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:

ex(r)=∑GeiG·r∥ex(G,z)

_ex(G,z)=−_∫1BZ^d2k^{(2π)2d2 k′(2π)2}∑α,^α′_{Θkα Θ^{k′ α′}}∑^G′_∫−∞∞d^z′2π^{e−|k−k′+G′||z−z′|}|k−^k′+^G′|×∑^G′′_{c^{k′ α′}*}(^G′′−^G′,z)_ckα(^G′′+G,z)∑^G′′′_ckα*(^G′′′,^z′)_{c^{k′ α′}}(^G′′′−^G′,^z′).shif

In order to extend the ^{k′ k′} -integration to the complete 2D k-space, one uses the fact that the coefficients _ckα(G) _ckα(G) only depend on k+G k+G,

_ckα(G,z)=_ck+G,α(0,z).

If one extends the occupation function accordingly,

_Θk+G,α:=_Θkα,

the exchange energy density can be expressed as:

_ex(G,z)=−_∫1BZ^d2k^(2π)2∫^{d2 k′(2π)2}∑α,^α′_{Θkα Θ^{k′ α′}∫−∞∞}d^z′2π^{e−|k′−k||z−z′|}|^k′−k|×∑^G′_{c^{k′ α′}*}(^G′,z)_ckα(^G′+G,z)∑^G′′_ckα*(^G′′,^z′)_{c^{k′ α′}}(^G′′,^z′).

At this point, one can choose the origin of the ^{k′ k′} -integration to be at k k,

_ex(G,z)=−∫^{d2 k′(2π)2}_∫−∞∞d^z′2π^{e−|k′||z−z′|}|^k′|f(^k′,G,z,^z′),shift

with:

f(^k′,G,z,^z′):=_∫1BZ^d2k^(2π)2∑α,^α′_{Θkα Θk+^{k′ α′}}∑^G′_{ck+^{k′ α′}*}(^G′,z)_ckα(^G′+G,z)×∑^G′′_ckα*(^G′′,^z′)_{ck+^{k′ α′}}(^G′′,^z′).shif

In the following, the behavior of _{ex ex} in the regime z≫L z≫Lis 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,

_ex(G,z)=_exL(G,z)+_exS(G,z)

_exL(G,z)=−∫^{d2 k′}2π_∫−∞∞d^{z′e−|k′|(|z−z′|+μ)}|^k′|f(0,G,z,^z′)

_exS(G,z)=−∫^{d2 k′}2π_∫−∞∞d^{z′e−|k′||z−z′|}|^k′|f(^k′,G,z,^z′)−f(0,G,z,^z′)^e−μ|k′|.

The ^{k′ k′}-integral in (25) can be performed directly,

_exL(G,z)=−_∫−∞zd^z′f(0,G,z,^z′)z−^z′+μ−_∫z∞d^z′f(0,G,z,^z′)^z′−z+μ.shi

For large z, the overlap function f(^k′,G,z,^z′) f(^k′,G,z,^z′) vanishes exponentially with z. This is immediately clear for the integrand in (23), since, in view of Equation (14), _ckα(0,z) _ckα(0,z) has the form (11) for large z, and all _ckα(G,z) _ckα(G,z) decay even faster than _ckα(0,z) _ckα(0,z) . Moreover, this exponential decay is preserved by the k k -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,

_exL(G,z)→z≫L−_∫−∞zd^z′f(0,G,z,^z′)z−^z′+μ.

The function f(^k′,G,z,^z′) f(^k′,G,z,^z′) also vanishes exponentially for ^z′→∞ ^z′→∞ , 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 ^{z′ z′} are large. The exponential decay of f(^k′,G,z,^z′) f(^k′,G,z,^z′) with large ^{z′ z′} allows one to expand the denominator in powers of ^z′/(z+μ) ^z′/(z+μ) . Subsequently, one can extend the ^{z′ z′}-integration to infinity, without introducing additional power-law corrections,

_exL(G,z)→z≫L−_∫−∞∞d^z′f(0,G,z,^z′)z+μ1+^z′z+μ+….shif

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

_∫−∞∞d^z′f(0,G,z,^z′)=_∫1BZ^d2k^(2π)2∑α_Θkα∑^G′_ckα*(^G′,z)_ckα(^G′+G,z).

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

n(r)=2∑GeiG·r∥ ∫1BZd2k(2π)2∑αΘkα∑G′ckα*(G′,z)ckα(G′+G,z).shift

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

_exL(G,z)→z≫L−n(G,z)2(z+μ)−p(G,z)2^(z+μ)2.shif

with:

p(G,z)=2_∫−∞∞d^z′z′f(0,G,z,^z′).shif

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

_exL(G,z)→z≫L−n(G,z)2z1−μz−p(G,z)2^z2+….shif

Now, consider the short-wavelength component of _{ex ex} , Equation (26). In order to analyze its behavior for large z, the ^{k′ k′} -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 Q Q , 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 Q Qare the reciprocal lattice vectors of the associated 2D simple cubic lattice,

_exS(G,z)=−∑Q_∫V^{d2 k′}2π_∫−∞∞d^{z′e−|k′−Q||z−z′|}|^k′−Q|f(^k′−Q,G,z,^z′)−f(0,G,z,^z′)^{e−μ|k′−Q|}.

As soon as Q≠0 Q≠0 , there is no Coulomb singularity in the ^{k′ k′} -integral: in this case, the ^{k′ k′} -integral is well defined for arbitrary z,^z′ z,^z′ . For z→∞ z→∞ , these contributions decay exponentially faster than the term with Q=0 Q=0, so that only the latter term remains to be analyzed,

_exS(G,z)→z≫L−_∫V^{d2 k′}2π_∫−∞∞d^{z′e−|k′||z−z′|}|^k′|f(^k′,G,z,^z′)−f(0,G,z,^z′)^e−μ|k′|.shif

For large z, an expansion of the expression in square brackets about ^k′=0 ^k′=0 is legitimate, since: (i) the ^{k′ k′}-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,

_exS(G,z)→z≫L−_∫V^{d2 k′}2π_∫−∞∞d^{z′e−|k′||z−z′|}|^k′|×_{∂f(^k′,G,z,^z′)∂^k′k′=0}·^k′+f(0,G,z,^z′)μ|^k′|+….shif

Due to the inversion symmetry of V, one has:

_∫V^d2k2πk|k|^{e−|k||z−z′|}=0.

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

_∫V^d2k2π^{e−|k||z−z′|}=_∫02πdφ2π_∫0S(φ)kdk^{e−|z−z′|k}=_∫02πdφ2π1|z−^{z′ |2}−^{e−|z−z′|S}|z−^{z′ |2}1+|z−^z′|S,shift

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

S(φ)=dcos(φ)0≤φ≤π4dsin(φ)π4≤φ≤3π4−dcos(φ)3π4≤φ≤5π4−dsin(φ)5π4≤φ≤7π4dcos(φ)7π4≤φ≤2π.

The leading contribution to _{exS exS}is thus obtained as

_exS(G,z)→z≫L−μ_∫−∞∞d^z′f(0,G,z,^z′)|z−^{z′ |2}_∫02πdφ2π1−^{e−|z−z′|S}1+|z−^z′|S.shift

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

1−^{e−|z−z′|S}1+|z−^z′|S=1−1−|z−^z′|S+(|z−^{z′ |S)2}2+…1+|z−^z′|S=(|z−^{z′ |S)2}2+….

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

_exS(G,z)→z≫L−μn(G,z)2^z2,

since S>0 S>0 for all φ φ and f(0,G,z,^z′) f(0,G,z,^z′) decays exponentially with ^{z′ z′} . Upon addition of (33) and (39), one finds for the leading two orders of the asymptotic energy density,

_ex(G,z)→z≫L−n(G,z)2z−p(G,z)2^z2.

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

_ex(r)→z≫L−n(r)2z−p(r)2^z2.

4. Asymptotic Behavior of Slater Potential

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

_vSlater(r)=2_ex(r)n(r)→z≫L−1z−p(r)n(r)^z2.

In order to evaluate the next-to-leading order term further, the k k -integrations involved in p(r) p(r) and n(r) n(r) need to be analyzed. Asymptotically, the density (30) is dominated by the coefficients _{ckα ckα} with the lowest exponents _γkα(G) _γkα(G) for which _Θkα>0 _Θkα>0 . If the lowest value of _γkα(G) _γkα(G) is found well inside the first BZ (which is usually the case [43]), the coefficients _ckα(G=0,z) _ckα(G=0,z) decay most slowly, as discussed in Section 2,

n(r)→z≫L2_∫1BZ^d2k^(2π)2∑α_Θkα^{|_ckα(0,z)|2},

with _ckα(0,z) _ckα(0,z) given by (11), due to (14). The same argument can be applied to the asymptotic p(r) p(r),

p(r)=∑G^eiG·_r∥p(G,z)→z≫L2_∫−∞∞d^z′z′_∫1BZ^d2k^(2π)2∑α,^α′_{Θkα Θk^α′ck^α′*}(0,z)_ckα(0,z)∑^G′′_ckα*(^G′′,^z′)_ck^α′(^G′′,^z′).shi

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

p(r)→z≫L2_∫1BZ^d2k^(2π)2∑α_Θkαmkα^{|_ckα(0,z)|2}.

Here, _{mkα mkα} denotes the first moment of |_ckα ^(G,z)|2 |_ckα ^(G,z)|2,

_mkα=_∫−∞∞dzz∑G^{|_ckα(G,z)|2},

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 _{ckα ckα} 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 _{Akα Akα} in Equation (18) of [43] corresponds to 2|_f0,kα ^|2 2|_f0,kα ^|2 , for the first moment _{Akα Akα} corresponds to 2_mkα ^|_f0,kα|2 2_mkα ^|_f0,kα|2 . Both (43) and (45) are asymptotically dominated by the state(s) with the minimum exponent _γkα(G=0) _γkα(G=0) , as discussed in [43]. If there is only a single k k -point q q and band β β for which _γkα(G=0) _γkα(G=0)reaches its minimum value inside the integration region, one finds:

p(r)→z≫L_mqβ2π^{|_f0,qβ(G=0)|2z2u/_γqβ−1e−2_γqβz detΓ−1/2}shift

n(r)→z≫L12π^{|_f0,qβ(G=0)|2z2u/_γqβ−1e−2_γqβz detΓ−1/2},shift

where Γ Γis defined by:

Γ=^∂2 _γkβ(G=0)∂k∂k(q).

p(r) p(r) therefore has the same asymptotic behavior as n(r) n(r) , so that p/n→z→∞_mqβ p/n→z→∞_mqβ . The ratio p/n p/n also approaches a constant for large z, if _f0,qβ(G=0)=0 _f0,qβ(G=0)=0 , 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 k k -point and/or more than one band for which _γkα(G=0) _γkα(G=0) assumes its minimum value (see [43]). One thus finally arrives at:

_vSlater(r)=2_ex(r)n(r)→z≫L−1z−_mqβ ^z2+…,

with _{mqβ mqβ} 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 k=q k=q under the transformation z→−z z→−z , _{mqβ mqβ} can vanish. This is the case in particular, if one has reflection symmetry with respect to the xy xy -plane, so that |_cqβ(G,−z)|=|_cqβ(G,z)| |_cqβ(G,−z)|=|_cqβ(G,z)|.

5. Asymptotic Behavior of Exact Exchange Potential

For a discussion of the exact _{vx vx} , one starts with the total exchange energy per unit cell, expressed in terms of the amplitudes _ckα(G,z) _ckα(G,z),

_Ex=−A_∫1BZ^d2k^{(2π)2d2 k′(2π)2}∑α,^α′_{ΘkαΘ^{k′ α′}}∑G_∫−∞∞dz_∫−∞∞d^z′2π^{e−|k−k′+G||z−z′|}|k−^k′+G|×∑^G′_{c^{k′ α′}*}(^G′,z)_ckα(^G′+G,z)∑^G′′_ckα*(^G′′,^z′)_{c^{k′ α′}}(^G′′−G,^z′)

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

δn(r)=2∑GeiG·r∥ ∫1BZd2k(2π)2∑αΘkα∑G′δckα*(G′,z)ckα(G′+G,z)+c.c..shi

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

Now, consider the variation of _{Ex Ex} under the variation of _ckα(G,z) _ckα(G,z),

δ_Ex=−2A_∫1BZ^d2k^{(2π)2d2 k′(2π)2}∑α,^α′_{ΘkαΘ^{k′ α′}}∑G_∫−∞∞dz_∫−∞∞d^z′2π^{e−|k−k′+G||z−z′|}|k−^k′+G|×∑^G′′_{c^{k′ α′}*}(^G′′,^z′)_ckα(^G′′+G,^z′)∑^G′δ_ckα*(^G′,z)_{c^{k′ α′}}(^G′−G,z)+c.c..shif

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

δ_Ex=−2A_∫1BZ^d2k^(2π)2∫^{d2 k′(2π)2}∑α,^α′_{ΘkαΘk+^k′,^α′∫−∞∞}dz_∫−∞∞d^z′2π^{e−|k′||z−z′|}|^k′|×∑^G′′_{ck+^k′,^α′*}(^G′′,^z′)_ckα(^G′′,^z′)∑^G′δ_ckα*(^G′,z)_{ck+^k′,^α′}(^G′,z)+c.c..shif

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

δ_ExL=−2A_∫1BZ^d2k^(2π)2∫^{d2 k′(2π)2}∑α,^α′_{ΘkαΘk^α′∫−∞∞}dz_∫−∞∞d^z′2π^{e−|k′|(|z−z′|+μ)}|^k′|×∑^G′′_ck^α′*(^G′′,^z′)_ckα(^G′′,^z′)∑^G′δ_ckα*(^G′,z)_ck^α′(^G′,z)+c.c.

δ_ExS=δ_Ex−δ_ExL.shif

The evaluation of δ_ExL δ_ExL proceeds in straightforward fashion. One first performs the ^{k′ k′}-integral,

δ_ExL=−2A_∫1BZ^d2k^(2π)2∑α,^α′_{ΘkαΘk^α′∫−∞∞}dz_∫−∞∞d^z′1|z−^z′|+μ×∑^G′′_ck^α′*(^G′′,^z′)_ckα(^G′′,^z′)∑^G′δ_ckα*(^G′,z)_ck^α′(^G′,z)+c.c..shif

One can then follow the argument leading from (27) to (29), since _ck^α′(^G′,z) _ck^α′(^G′,z) falls off exponentially with z, irrespective of δ_ckα*(^G′,z) δ_ckα*(^G′,z),

δ_ExL→z≫L−2A_∫1BZ^d2k^(2π)2∑α,^α′_{ΘkαΘk^α′∫−∞∞}dz_∫−∞∞d^z′1z+μ1+^z′z+μ+…×∑^G′′_ck^α′*(^G′′,^z′)_ckα(^G′′,^z′)∑^G′δ_ckα*(^G′,z)_ck^α′(^G′,z)+c.c..shif

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

δ_ExL0→z≫L−2A_∫1BZ^d2k^(2π)2∑α_{Θkα ∫−∞∞}dzz+μ∑^G′δ_ckα*(^G′,z)_ckα(^G′,z)+c.c..shif

One can now re-introduce the x-y-integrations and identify δn(r) δn(r),

δExL0→z≫L−∫A d2 r∥ ∫−∞∞dzz+μ×2∑GeiG·r∥ ∫1BZd2k(2π)2∑αΘkα∑G′δckα*(G′,z)ckα(G′+G,z)+c.c.=−∫A d2 r∥ ∫−∞∞dzδn(r)z+μ.shift

Finally, an expansion in powers of μ/z μ/z is legitimate, since μ μ can be chosen to be much smaller than the z-values for which δ_ckα(G,z) δ_ckα(G,z)is non-zero,

δ_ExL0→z≫L−_∫A ^d2 _{r∥ ∫−∞∞}dzδn(r)z1−μz+….

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

δ_ExS=−2A_∫1BZ^d2k^(2π)2∫^{d2 k′(2π)2}∑α,^α′_{Θkα∫−∞∞}dz_∫−∞∞d^z′2π^{e−|k′||z−z′|}|^k′|×[_{Θk+^k′,^α′}∑^G′′_{ck+^k′,^α′*}(^G′′,^z′)_ckα(^G′′,^z′)∑^G′δ_ckα*(^G′,z)_{ck+^k′,^α′}(^G′,z)×−_Θk^α′∑^G′′_ck^α′*(^G′′,^z′)_ckα(^G′′,^z′)∑^G′δ_ckα*(^G′,z)_ck^α′(^G′,z)^e−|k′|μ]+c.c..shif

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

δ_ExS→z≫L−2Aμ_∫1BZ^d2k^(2π)2_∫V^{d2 k′(2π)2}∑α,^α′_{ΘkαΘk^α′∫−∞∞}dz_∫−∞∞d^z′2π^{e−|k′||z−z′|}×∑^G′′_ck^α′*(^G′′,^z′)_ckα(^G′′,^z′)∑^G′δ_ckα*(^G′,z)_ck^α′(^G′,z)+c.c..shif

One can then employ (37) for the ^{k′ k′} -integration and expand in powers of 1/z 1/z , since the remaining ^{z′ z′}-integral is free of any singularities,

δ_ExS→z≫L−2Aμ_∫1BZ^d2k^(2π)2∑α,^α′_{ΘkαΘk^α′∫−∞∞}d^z′∑^G′′_ck^α′*(^G′′,^z′)_ckα(^G′′,^z′)×_∫−∞∞dz^z2∑^G′δ_ckα*(^G′,z)_ck^α′(^G′,z)+c.c..shif

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),

δ_ExL1→z≫L−2A_∫1BZ^d2k^(2π)2∑α,^α′_{ΘkαΘk^α′∫−∞∞}dz^z2∑^G′δ_ckα*(^G′,z)_ck^α′(^G′,z)×_∫−∞∞d^z′z′∑^G′′_ck^α′*(^G′′,^z′)_ckα(^G′′,^z′)+c.c..shif

A completely general representation of δ_ExL1 δ_ExL1 in terms of δn δn , as achieved for δ_ExL0 δ_ExL0 , 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 k k -point(s) for which the exponent _γkα(G) _γkα(G) , 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 _γkα(G) _γkα(G) . Assume, for simplicity, that the lowest value of the exponent is obtained only for a single state, denoted by qβ qβ . If only the amplitudes _{ckβ ckβ} in the vicinity of q q are varied, while all other _{ckα ckα} are kept unchanged, δn δnis given by:

δn(r)=2∑GeiG·r∥ ∫1BZd2k(2π)2∑G′δckβ*(G′,z)ckβ(G′+G,z)+c.c..shif

Now, consider δ_ExL1 δ_ExL1 for this variation. For any given δ_ckα*(^G′,z) δ_ckα*(^G′,z) , i.e., any given state kα kα , the sum over ^{α′ α′} in (65) is dominated by the most weakly-decaying occupied amplitude among the _ck^α′(^G′,z) _ck^α′(^G′,z) . In the case of the variation (66), δ_ExL1 δ_ExL1thus reduces to:

δ_ExL1→z≫L−2A_∫1BZ^d2k^(2π)2_{mkβ ∫−∞∞}dz^z2∑^G′δ_ckβ*(^G′,z)_ckβ(^G′,z)+c.c.,shif

with _{mkβ mkβ} defined by Equation (46). Since δ_ckβ(^G′,z) δ_ckβ(^G′,z) is non-zero only for large z and _ckβ(^G′,z) _ckβ(^G′,z) decays exponentially for these z-values, the k k -integration in (67) can be simplified by a Taylor expansion of the smooth parts of the integrand about k=q k=q,

δ_ExL1→z≫L−2A_{mqβ ∫1BZ}^d2k^(2π)2_∫−∞∞dz^z2∑^G′δ_ckβ*(^G′,z)_ckβ(^G′,z)+c.c..shif

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

δ_ExL1→z≫L−_{mqβ ∫A} ^d2 _{r∥ ∫−∞∞}dzδn(r)^z2.

Thus, for the relevant class of variations δ_ExL1 δ_ExL1 can be explicitly expressed in terms of δn δn. 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 1/^z2 1/^z2,

_vx(r)=δ_Exδn(r)→z≫L−1z−_mqβ ^z2+…,

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 k k -point q qin the interior of the first BZ, one finds:

_ex(r)→z≫Ln(r)2−1z−_mqβ ^z2+O^z−3

_vSlater(r)→z≫L−1z−_mqβ ^z2+O^z−3

_vx(r)→z≫L−1z−_mqβ ^z2+O^z−3,

where _{mqβ mqβ} denotes the first moment (46) of the most weakly-decaying occupied state. If there are several degenerate most weakly-decaying states at the k k -point q q or states at several k k -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, ^z′=z+Δ ^z′=z+Δ . At the same physical point far from the surface, Equation (73) then yields in the two coordinate systems,

−1^z′−_mqβ′ ^(z′)2+O^(z′)−3=!−1z−_mqβ ^z2+O^z−3,

where _{mqβ′ mqβ′} denotes the moment of the state qβ qβevaluated with respect to the primed coordinate system,

_mqβ′=_∫−∞∞d^{z′ z′}∑G^{|_cqβ′(G,z′)|2}=_∫−∞∞d^{z′ z′}∑G^{|_cqβ(G,z′−Δ)|2}=_∫−∞∞dz(z+Δ)∑G^{|_cqβ(G,z)|2}=_mqβ+Δ.

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

−1^z′−_mqβ′ ^(z′)2=−1z1−Δz−_mqβ+Δ^z2+O^z−3=−1z−_mqβ ^z2+O^z−3.

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 _{mqβ mqβ} ; a highly localized most weakly-decaying state right at one of the two surfaces of the slab would yield |_mqβ|≈L |_mqβ|≈L . However, as Equation (75) shows, _{mqβ mqβ} 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.

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

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

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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

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.

16. Mak, K.F.; Lee, C.; Hone, J.; Shan, J.; Heinz, T.F. Atomically thin MoS2: A new direct-gap semiconductor. Phys. Rev. Lett. 2010, 105, 136805.

17. Radisavljevic, B.; Radenovic, A.; Brivio, J.; Giacometti, V.; Kis, A. Single-layer MoS2 transistors. Nat. Nanotechnol. 2011, 6, 147.

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.

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

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

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.

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.

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.

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

25. Qian, Z.; Sahni, V. Quantum mechanical image potential theory. Phys. Rev. B 2002, 66, 205103.

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.

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

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.

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.

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.

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.

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.

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

34. Engel, E. Exact exchange plane-wave-pseudopotential calculations for slabs. J. Chem. Phys. 2014, 140, 18A505.

35. Engel, E. Asymptotic behavior of exact exchange potential of slabs. Phys. Rev. B 2014, 89, 245105.

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.

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

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.

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

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.

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.

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.

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

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.

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.

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.

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