Exact exchange plane-wave-pseudopotential calculations for slabs: Extending the width of the vacuum

Standard plane-wave pseudopotential (PWPP) calculations for slabs such as graphene become extremely demanding, as soon as the exact exchange (EXX) of density functional theory is applied. Even if the Krieger-Li-Iafrate (KLI) approximation for the EXX potential is utilized, such EXX-PWPP calculations suffer from the fact that an accurate representation of the occupied states throughout the complete vacuum between the replicas of the slab is required. In this contribution, a robust and efficient extension scheme for the PWPP states is introduced, which ensures the correct exponential decay of the slab states in the vacuum for standard cutoff energies and therefore facilitates EXX-PWPP calculations for very wide vacua and rather thick slabs. Using this scheme, it is explicitly verified that the Slater component of the EXX/KLI potential decays as $-1/z$ over an extended region sufficiently far from the surface (assumed to be perpendicular to the $z$ direction) and from the middle of the vacuum, thus reproducing the asymptotic behavior of the exact EXX potential of a single slab. The calculations also reveal that the orbital-shift component of the EXX/KLI potential is quite sizable in the asymptotic region. In spite of the long-range exchange potential, the replicas of the slab decouple rather quickly with increasing width of the vacuum. Relying on the identity of the work function with the Fermi energy obtained with a suitably normalized total potential, the present EXX/KLI calculations predict work functions for both graphene and the Si(111) surface which are substantially larger than the corresponding experimental data. Together with the size of the orbital-shift potential in the asymptotic region, the very large EXX/KLI work functions indicate a failure of the KLI approximation for nonmetallic slabs.

Figure 1

Densities $|{\varphi }_{\mathit{k}\alpha }{\left(\mathit{r}\right)|}^2$ of occupied states of graphene as functions of $z$ for $x=\frac{9a}{128}$ and $y=\frac{\sqrt{3}a}{32}$: (a) $\mathit{k}={\mathit{k}}_{\mathrm{a}}$, Eq. (45); (b) $\mathit{k}={\mathit{k}}_{\mathrm{b}}$, Eq. (46). Both the original numerical state (23) and its extension by the SE approach are shown. $z_0$ is defined by (44) with $s={10}^{-8}{\text{Bohr}}^{-3}$ ($d=10a\approx 46.5\mathrm{Bohr}, E_{\mathrm{cut}}=1880\mathrm{Ry}$).

Figure 2

Density $|{\varphi }_{\mathit{k}\alpha }{\left(\mathit{r}\right)|}^2$ of state $\alpha =3$ at $\mathit{k}={\mathit{k}}_{\mathrm{b}}$, Eq. (46), for different transition points from the numerical Fourier expanded state (23) to the analytical asymptotic form (27), using the SE scheme: comparison of transition points 4.02, 5.42, 6.89, and 9.04 Bohr resulting from $s={10}^{-n}{\text{Bohr}}^{-3}$ with $n=6,8,10,\text{and}12$, respectively ($x, y, d$, and $E_{\mathrm{cut}}$ as in Fig. 1).

Figure 3

Logarithmic derivative (47) of state $\alpha =4$ at $\mathit{k}={\mathit{k}}_{\mathrm{a}}$, Eq. (45), for different transition points from the numerical Fourier expansion (23) to the analytical asymptotic form (27), using the SE scheme: comparison of transition points 6.57, 8.72, 10.84, and 12.88 Bohr resulting from $s={10}^{-n}{\text{Bohr}}^{-3}$ with $n=6,8,10,\text{and}12$, respectively ($x, y, d$, and $E_{\mathrm{cut}}$ as in Fig. 1). Also shown is the LO extension for $s={10}^{-8}{\text{Bohr}}^{-3}$.

Figure 4

Ratio of extension (27) to numerical wave function (23) in the vicinity of the transition point for Bloch state $\alpha =4$ at $\mathit{k}={\mathit{k}}_{\mathrm{a}}$, Eq. (45): (a) LO, SE, and UI extensions including only the star of the lowest nonvanishing $\mathit{G}$-vectors vs SE including the stars of the four lowest $\mathit{G}\ne \mathit{0}$; (b) UI versus CI and CR extensions for differently many $\mathit{G}\ne \mathit{0}$ ($z_0=8.72\mathrm{Bohr}, x, y, d$, and $E_{\mathrm{cut}}$ as in Fig. 1).

Figure 5

Logarithmic derivative (47) of state $\alpha =3$ at $\mathit{k}={\mathit{k}}_{\mathrm{b}}$, Eq. (46) for different transition points from the numerical Fourier expansion (23) to the analytical asymptotic wave function (27) obtained with the SE scheme: comparison of transition points 4.02, 5.42, 6.89, and 9.04 Bohr resulting from $s={10}^{-n}{\text{Bohr}}^{-3}$ with $n=6,8,10,\text{and}12$, respectively ($x, y, d$, and $E_{\mathrm{cut}}$ as in Fig. 1). Also shown is the UI scheme for $s={10}^{-8}{\text{Bohr}}^{-3}$.

Figure 6

Derivative of $x-y$-averaged $v_{\mathrm{shift}}$, Eq. (34), for graphene: comparison of different transition points, resulting from UI scheme and density thresholds $s={10}^{-n}{\text{Bohr}}^{-3}$. (a) vacuum width $d=7a, n=6-13$; (b) vacuum width $d=10a, n=6-12$. The position of the number indicates the associated transition point ($E_{\mathrm{cut}}=1880\mathrm{Ry}, 8\times 8\mathit{k}$-point grid).

Figure 7

Derivative of $x-y$-averaged $v_{\mathrm{shift}}$, Eq. (34), for graphene: comparison of different extension schemes for (a) the transition region ($d=10a, E_{\mathrm{cut}}=1880\mathrm{Ry}$) and (b) very large $z$ ($d=40a, E_{\mathrm{cut}}=520\mathrm{Ry}$) (in all cases $s={10}^{-8}{\text{Bohr}}^{-3}, 8\times 8\mathit{k}$-point grid).

Figure 8

$x-y$-averaged $v_{\mathrm{shift}}$, Eq. (34), for graphene: comparison of different extension schemes ($s={10}^{-8}{\text{Bohr}}^{-3}, E_{\mathrm{cut}}=1880\mathrm{Ry}, 8\times 8\mathit{k}$-point grid).

Figure 9

First derivative of $x-y$-averaged $v_{\mathrm{Slater}}$, Eq. (29), for a 6-layer Si(111) slab with $d=82a/\sqrt{3}$ ($32\times 32\mathit{k}$-point grid, LO extension).

Figure 10

Convergence of (a) $x-y$-averaged $v_{\mathrm{Slater}}$, Eq. (29), and (b) its first derivative with respect to $z$ with $\mathit{k}$-point sampling. The $\mathit{k}$-point grids used are characterized by the Monkhorst-Pack parameters [57] ($d=10a=46.5\mathrm{Bohr}, E_{\mathrm{cut}}=520\mathrm{Ry}$, SE scheme).

Figure 11

As Fig. 10 for $d=40a=186\mathrm{Bohr}$. For $z<10\mathrm{Bohr}$, the potentials of panel (a) can not be distinguished from the corresponding results in Fig. 10.

Figure 12

Convergence of (a) $x-y$-averaged $v_{\mathrm{shift}}$, Eq. (34), and (b) its first derivative with respect to $z$ with $\mathit{k}$-point sampling. The $\mathit{k}$-point grids used are characterized by the Monkhorst-Pack parameters [57] ($d=10a=46.5\mathrm{Bohr}, E_{\mathrm{cut}}=520\mathrm{Ry}$, SE scheme).

Figure 13

As Fig. 12 for $d=40a=186\mathrm{Bohr}$. For $z<20\mathrm{Bohr}$, both the potentials and their derivatives can not be distinguished from the corresponding results in Fig. 12.

Figure 14

Convergence of (a) $x-y$-averaged $v_{\mathrm{Slater}}$, Eq. (29), and (b) its first derivative with respect to $z$ with the width of vacuum $d$. All results have been obtained with the SE scheme. For $z<5\mathrm{Bohr}$, the potentials of panel (a) can not be distinguished from the result in Fig. 10.

Figure 15

Convergence of (a) $x-y$-averaged $v_{\mathrm{shift}}$, Eq. (34), and (b) its first derivative with respect to $z$ with the width of vacuum $d$. All results have been obtained with the $24\times 24\mathit{k}$-point grid and the SE scheme. For $z<10\mathrm{Bohr}$, both the potentials and their derivatives can not be distinguished from the corresponding results in Fig. 12.

Figure 16

$x-y$-averaged $v_{\mathrm{x}}^{\mathrm{KLI}}$, Eq. (28), of graphene in comparison to $-1/z$ ($32\times 32\mathit{k}$-point grid, CR scheme).

Figure 17

$x-y$-averaged $v_{\mathrm{x}}^{\mathrm{KLI}}$, Eq. (28), of a 6-layer Si(111) slab with $d=82a/\sqrt{3}$ in comparison to $-1/z$ ($32\times 32\mathit{k}$-point grid, CI scheme).

Figure 18

Fermi energy ${\epsilon }_{\mathrm{F}}$ of graphene: convergence with vacuum width $d$ for different $\mathit{k}$-point grids. Lines are drawn to guide the eye.

Figure 19

EXX/KLI band energies ${\epsilon }_{\mathit{k}\alpha }$ of graphene at $\mathrm{\Gamma }$-point: convergence with vacuum width $d$ for different $\mathit{k}$-point grids. The eigenvalues are normalized to the highest valence-band energy. Lines are drawn to guide the eye.

Figure 20

EXX/KLI band structure of graphene: convergence with vacuum width $d$.

Figure 21

$xy$-averaged exchange potential: EXX/KLI result versus LDA (x-only) for (a) graphene and (b) a 6-layer Si(111) slab.

Figure 22

Densities $|{\varphi }_{\mathit{k}\alpha }{\left(\mathit{r}\right)|}^2$ of valence states of a 6-layer Si(111) slab as functions of $z$ for $x=\frac{9a}{40\sqrt{2}}$ and $y=\frac{\sqrt{6}a}{20}$; $\mathit{k}=\frac{2\pi }{a}\left(1/\sqrt{32},0\right)$. Both the original numerical state (23) and its extension by the UI approach are shown ($s={10}^{-8}{\text{Bohr}}^{-3}, d=12a/\sqrt{3}\approx 70.9\mathrm{Bohr}$, and $E_{\mathrm{cut}}=320\mathrm{Ry}$).

Figure 23

Derivative of $x-y$-averaged $v_{\mathrm{shift}}$, Eq. (34), for a 6-layer Si(111) slab with $d=12a/\sqrt{3}$: comparison of LO, SE, UI, and CR results for $s={10}^{-8}{\text{Bohr}}^{-3}$ ($E_{\mathrm{cut}}=80\mathrm{Ry}, 8\times 8\mathit{k}$-point grid) with LO potential for $s={10}^{-12}{\text{Bohr}}^{-3}$ ($E_{\mathrm{cut}}=640\mathrm{Ry}, 8\times 8\mathit{k}$-point grid). The result of the CI scheme is indistinguishable from the UI data.

Figure 24

As Fig. 23 for $v_{\mathrm{shift}}$ itself.

Figure 25

(a) $x-y$-averaged $v_{\mathrm{shift}}$, Eq. (34), for a 6-layer Si(111) slab with $d=82a/\sqrt{3}$ and (b) its first derivative: comparison of LO, SE, UI, and CR extensions ($s={10}^{-8}{\text{Bohr}}^{-3}, E_{\mathrm{cut}}=80\mathrm{Ry}, 8\times 8\mathit{k}$-point grid). The result of the CI scheme is indistinguishable from the UI data.