Finite-size corrections for defect-involving vertical transitions in supercell calculations

TL;DR

This paper introduces a correction method for defect-related vertical transition levels in supercell calculations, significantly reducing errors from about 1 eV to below 0.12 eV across multiple materials.

Contribution

A general correction scheme for defect-involving vertical transition levels in supercell calculations is formulated and validated on various materials.

Findings

Errors reduced from ~1 eV to <0.12 eV with correction

Method verified on ten different defects in multiple materials

Potential for wide application in defect level calculations

Abstract

A correction method for vertical transition levels (VTLs) involving defect states calculated with a supercell technique is formulated and its effectiveness is systematically verified with ten defects in prototypical materials: cubic-BN, GaN, MgO, and 3C-SiC. Without any corrections, the absolute errors are around 1 eV with moderate size supercells in most cases. In contrast, when our correction method is adopted, the absolute errors are reduced and become less than 0.12 eV in all the cases. Our correction scheme is general and will have the potential for wide application as it is adaptive for evaluating various quantities at fixed geometry, as represented by those relevant to the generalized Koopmans' theorem.

Equations

$$E_{\mathrm{corr}}=E_{\mathrm{PC}}-Q\Delta V_{\mathrm{PC},Q/b}{|}_{\mathrm{far}}.$$

$$\begin{array}{rl}\mu [Q/Q+\Delta Q;S(Q)]=& \left\{E\left[Q+\Delta Q;S(Q)\right]-E\left[Q;S(Q)\right]\right\}\\ & +\Delta Q\cdot {\epsilon }_{\mathrm{V}(\mathrm{C})\mathrm{BM}},\end{array}$$

$$\Delta E=\int d^3\mathbf{r}\Delta q\left(\mathbf{r}\right)V^Q\left(\mathbf{r}\right)+\frac{1}{2}\int d^3\mathbf{r}\Delta q\left(\mathbf{r}\right)V^{\Delta Q}\left(\mathbf{r}\right),$$

$$V_{\mathrm{sr}}^Q\left(\mathbf{r}\right)$$

$$V_{\mathrm{sr}}^{\Delta Q}\left(\mathbf{r}\right)$$

$$V_{\mathrm{corr}}^Q\left(\mathbf{r}\right)$$

$$V_{\mathrm{corr}}^{\Delta Q}\left(\mathbf{r}\right)$$

$$\Delta E_{\mathrm{corr}}=\int d^3\mathbf{r}\Delta q\left(\mathbf{r}\right)V_{\mathrm{corr}}^Q\left(\mathbf{r}\right)+\frac{1}{2}\int d^3\mathbf{r}\Delta q\left(\mathbf{r}\right)V_{\mathrm{corr}}^{\Delta Q}\left(\mathbf{r}\right).$$

$$V_{\mathrm{corr}}^Q\left(\mathbf{r}\right)=\frac{\alpha Q}{{\epsilon }_{0}L}-\frac{2\pi \int d^3\mathbf{r}q\left(\mathbf{r}\right)r^2}{3{\epsilon }_0{L}^3}-\frac{2\pi Q{r}^2}{3{\epsilon }_0{L}^3}.$$

$$C^{\Delta Q}=\frac{2\pi \int d^3\mathbf{r}{\Delta \rho }^{\Delta Q}\left(\mathbf{r}\right)r^2}{3{{\epsilon }_{\infty }L}^3}$$

$$\int d^3\mathbf{r}{\Delta \rho }^{\Delta Q}{V}_{\mathrm{corr}}^Q\left(\mathbf{r}\right)={\frac{{\epsilon }_{\infty }}{{\epsilon }_0}QC}^{\Delta Q}.$$

$$\begin{array}{rl}\Delta E_{\mathrm{corr}}=& \frac{\alpha Q\Delta Q}{{\epsilon }_{0}L}+\frac{\alpha {\left(\Delta Q\right)}^2}{2{\epsilon }_{\infty }L}\\ & -\left(\Delta Q{C}^{\Delta Q}+\Delta Q{C}^Q+\frac{{\epsilon }_{\infty }}{{\epsilon }_0}Q{C}^{\Delta Q}\right).\end{array}$$

$$\begin{array}{rl}\Delta E_{\mathrm{corr}}=& \frac{2\Delta Q}{Q}{E}_{\mathrm{PC}}^Q\left({\bar{\epsilon }}_0\right)+E_{\mathrm{PC}}^{\Delta Q}\left({\bar{\epsilon }}_{\infty }\right)\\ & -\left(\Delta Q{C}^{\Delta Q}+\Delta Q{C}^Q+\frac{⟨{\epsilon }_{\infty }⟩}{⟨{\epsilon }_0⟩}Q{C}^{\Delta Q}\right),\end{array}$$

$${\Delta }_{\mathrm{XC}}=\{E\left[Q+\Delta Q;S(Q)\right]-E\left[Q;S(Q)\right]\}+\Delta Q\cdot \epsilon \left[Q;S(Q)\right],$$

Full text

Point defects determine various physical properties in solids. One of the most important properties is the optical property. For instance, ZnO is known to show green luminescence, which has been attributed to point defects Rodnyi (2011). Another example is SrTiO3, known to show blue-light emission at room temperature after Ar*+* irradiation, whose origin is presumably the emergence of oxygen vacancies Kan et al. (2005). In addition, we can investigate defects with deep states, which often degrade the device performance, from the photoabsorption and photoemission spectra. It is generally considered that absorption (emission) of a photon by a defect promotes (demotes) an electron to the excited (ground) state, most probably without altering the atomic configuration, based on the Franck-Condon principle Franck and Dymond (1926). This is a consequence of the fact that electrons are much lighter than nuclei. Thus, we can represent the optical transition by a vertical arrow in the configuration coordination diagram, and its transition energy is given as an optical transition level or vertical transition level (VTL) with respect to the valence band maximum (VBM) or the conduction band minimum (CBM) Van de Walle and Neugebauer (2004); Freysoldt et al. (2014).

First-principles calculations have become a powerful tool to understand and predict the defect properties. In the calculations of extended systems with defects, they are almost always evaluated by a supercell approach nowadays, where a charged defect interacts with its periodic images and background charge, which erroneously modifies the total energies of charged defect supercells Lany and Zunger (2008); Freysoldt et al. (2009); Komsa et al. (2012); Kumagai and Oba (2014). Methods to correct the energies to the dilute limit are well established, as represented by the scheme proposed by Freysoldt, Neugebauer, and Van de Walle (FNV) Freysoldt et al. (2009), and its extension to anisotropic systems and/or relaxed geometries (eFNV) Kumagai and Oba (2014). The correction energy of the (e)FNV scheme is written as

[TABLE]

The first term is the point charge (PC) correction energy and the second is an alignmentlike term, where $\Delta V_{\mathrm{PC,}Q\mathrm{/b}}|_{\mathrm{far}}$ is the potential difference between defect-induced and PC potentials at a region outside of the defect in the supercell, and Q is the defect charge state. The PC correction energy and potential are evaluated using the static dielectric tensor because a charged defect in its ground state is screened by both electrons and ions. The charged defect formation energy in the dilute limit is accurately calculated using Eq. (1) as long as the defect charge is enclosed in a supercell Freysoldt et al. (2009); Komsa et al. (2012); Kumagai and Oba (2014).

However, conventional corrections are not applicable to the calculation of a VTL as this involves a non-trivial combination of ionic and electronic screening, as we will explain in a moment. In general, a vertical transition starts from some charge state Q in the equilibrium atomic configuration S of that charge state, S(Q), which includes the ionic screening response to $Q$. The VTL to final charge state $Q+\Delta Q$, involving the VBM (CBM), is calculated as

[TABLE]

where $\Delta Q$ is the additional charge, i.e., $+1$ or $-1$, and ${\varepsilon}_{\mathrm{V(C)BM}}$ denotes the energy level of the VBM (CBM). Let us illustrate the issues of the VTL calculations from supercells through the example of the $+1$ to neutral charge transition of the oxygen vacancy (VO) in MgO. One would naively expect the VTL to be well estimated by applying the conventional FNV scheme to the initial state using the static dielectric tensor since the final state is charge neutral. Unfortunately, this is incorrect. As shown in Fig. 1(a), the eFNV-corrected VTLs vary even more with cell size than the uncorrected ones. As the total energy of the initial state should be well corrected by the FNV method, the cell size dependence is attributed to the energy of the final state although it is neutral. Indeed, the formation energies of V${}_{\mathrm{O}}^{0}$ at the V${}_{\mathrm{O}}^{+1}$ atomic configurations increase with increasing supercell size as shown in Fig. 1(b). This is an artifact purely caused by a spurious electrostatic potential. The reason is that the interaction of the added charge ($\Delta Q=-1$) with the periodic images of the initial state ($Q=+1$) are statically screened (electrons + ions), but the interactions with its own periodic images ($\Delta Q=-1$) involve electronic screening only, and the two contributions do not cancel.

Instead, one might think the (e)FNV scheme using the high-frequency (electronic contribution only) dielectric tensor works better as the energy caused by the ionic displacement is canceled out when taking the energy difference between the final and initial states. Again, this is not correct as shown in Fig. 1(a), where much worse cell size dependence remains.

Similar issues are of course present for any VTL. Although there have been numerous theoretical reports on the calculations of VTLs Lyons et al. (2012); Lyons and Van de Walle (2017); Frodason et al. (2017); Ho et al. (2018), surprisingly, none has referred to this issue. Therefore, we aim to formulate a correction method for the VTLs within the total energy framework.

To simplify the derivation, we will assume the cubic system with isotropic screening. However, the final correction energy expression is common to the anisotropic systems. A VTL or an additional electrostatic energy, arising from adiabatically introducing $\Delta$Q with a spatial distribution $\Delta$q($\bm{r}$) that excludes the screening charge, becomes for the isolated defect

[TABLE]

where $V^{Q}\left(\bm{r}\right)$ is the initial defect-induced electrostatic potential of a single defect, and $V^{\Delta Q}\left(\bm{r}\right)$ an additional contribution to the electrostatic potential, which is caused by the addition of ${\Delta Q}$. The factor 1/2 arises only for the second term as $V^{\Delta Q}\left(\bm{r}\right)$ builds up while $\Delta$Q is introduced, whereas $V^{Q}\left(\bm{r}\right)$ is unchanged throughout the process. Since Q and $\Delta$Q are screened with static (${\varepsilon}_{0}$) and electronic (${\varepsilon}_{\infty}$) dielectric constants, respectively, the asymptotic potential behaves like $\frac{Q/{\varepsilon}_{0}+\Delta Q/{\varepsilon}_{\infty}}{r}$ when the centers of Q and $\Delta$Q are assumed to be the same. Therefore, standard finite-size supercell corrections with a single charge and a single dielectric constant, cannot capture this effect. Also, it is not a difference of simple electrostatic corrections of initial and final state as demonstrated above.

In the dilute limit, $V^{Q}\left(\bm{r}\right)$ is written as $V^{Q}_{\mathrm{lr}}\left(\bm{r}\right)+V^{Q}_{\mathrm{sr}}\left(\bm{r}\right)$, where $V^{Q}_{\mathrm{lr}}\left(\bm{r}\right)$ is the long-range potential caused by point charge Q and $V^{Q}_{\mathrm{sr}}\left(\bm{r}\right)$ the remaining part representing the short-range potential. Similarly, $V^{\Delta Q}\left(\bm{r}\right)$ is written as $V^{\Delta Q}_{\mathrm{lr}}\left(\bm{r}\right)+V^{\Delta Q}_{\mathrm{sr}}\left(\bm{r}\right)$. Following the discussion in Ref. Freysoldt et al. (2009), the short-range potential can be obtained directly from the density functional theory (DFT) electrostatic potentials $\tilde{V}^{\mathrm{DFT}}$ as

[TABLE]

The symbol $\tilde{V}$ is used for periodic potentials defined within the DFT supercell. The initial state is the defect-containing supercell in charge state Q, while the final state is the one in charge state $Q+\Delta Q$. $C^{Q}$ and $C^{\Delta Q}$ are alignment constants chosen such that the short-range potentials decay to zero far from the defect Freysoldt et al. (2009, 2011). The correction potentials for charge Q and $\Delta$Q are therefore written as

[TABLE]

where $V^{Q,\mathrm{corr}}_{\mathrm{lr}}\left(\bm{r}\right)=V^{Q}_{\mathrm{lr}}\left(\bm{r}\right)-{\tilde{V}}^{Q}_{\mathrm{lr}}\left(\bm{r}\right)$ and $V^{\Delta Q,\mathrm{corr}}_{\mathrm{lr}}\left(\bm{r}\right)=V^{\Delta Q}_{\mathrm{lr}}\left(\bm{r}\right)-{\tilde{V}}^{\Delta Q}_{\mathrm{lr}}\left(\bm{r}\right)$ mean the sign-reversed Madelung potentials. The correction energy is then obtained by replacing the electrostatic potentials to the correction potentials

[TABLE]

Here, we split the additional charge $\Delta q\left(\bm{r}\right)$ to ${\rho}^{\Delta Q}_{\mathrm{PC}}$ and ${\Delta\rho}^{\Delta Q}$, where ${\rho}^{\Delta Q}_{\mathrm{PC}}$ is the point charge $\Delta Q$ at the center of the defect and ${\Delta\rho}^{\Delta Q}$ is the finite-size distribution. Based on the conventional FNV scheme, the second term in Eq. (8) is written as $\frac{\alpha{\left(\Delta Q\right)}^{2}}{2{\varepsilon}_{\infty}L}-{\Delta QC}^{\Delta Q}$ Freysoldt et al. (2009); Komsa et al. (2012), where L and $\alpha$ are the size of the supercell and the Madelung constant, respectively. In the same manner, the first term is rewritten as $\frac{\alpha Q\Delta Q}{{\varepsilon}_{0}L}-{\Delta QC}^{Q}-\int{d^{3}\bm{r}}{\Delta\rho}^{\Delta Q}V^{Q}_{\mathrm{corr}}\left(\bm{r}\right)$. Here, the third term indicates the interaction between the finite size of the additional charge and the correction potential for the initial charge. When we assume Coulomb interactions are screened by a dielectric constant independently of distance, $V^{Q}_{\mathrm{corr}}\left(\bm{r}\right)$ can be expanded as,

[TABLE]

Using

[TABLE]

as derived by Komsa et al. Komsa et al. (2012),

[TABLE]

By summing all the contributions, the correction energy can be written as

[TABLE]

For anisotropic systems, one needs to use the averaged dielectric constant, i.e., $\langle\varepsilon\rangle=\frac{1}{3}\left(\varepsilon_{\mathrm{11}}+\varepsilon_{\mathrm{22}}+\varepsilon_{\mathrm{33}}\right)$, where $\varepsilon_{ij}$ is a component of dielectric tensor $\bar{\varepsilon}$. Then, the correction energy can be rewritten as

[TABLE]

where $E^{Q}_{\mathrm{PC}}\left(\bar{\varepsilon}\right)$ is the point-charge corrections for charge Q screened by $\bar{\varepsilon}$. The first two contributions mean first-order point-charge correction, while the third term (in parentheses) corresponds to the alignment term.

To test our correction method, we investigate the supercell size dependence of the VTLs involving defects in prototypical nonmetallic materials, namely cubic-BN (c-BN), GaN, MgO, and 3C-SiC, with respect to the CBM and VBM [Eq. (2)] for donor- and acceptor-type defects, respectively. All the VTLs in this study were calculated with the Heyd-Scuseria-Ernzerhof hybrid functional Heyd et al. (2003), and the Fock-exchange parameters were tuned to reproduce each experimental band gap Varley et al. (2012); Gordon et al. (2015); Lyons and Van de Walle (2017); Weston et al. (2017). We also confirmed that the formation energies of all the defects considered in this study were accurately corrected with eFNV (see Fig. S1 in the Supplemental Material Not ), which indicated these defects were mostly enclosed by relatively smaller, e.g., 64-atom supercells.

Figures 2(a) and 2(b) show both the corrected and uncorrected VTLs of VB in c-BN as a function of the supercell size. Without any corrections, one can see the VTLs strongly depend on the supercell size, and the absolute error is larger than 0.7 eV even with the largest 1000-atom supercell in our test set. Conversely, the PC-level correction based on Eq. (13) drastically reduces the cell size dependences, and the absolute errors are only 0.1 eV at the smallest 64-atom supercell. Such a small cell size dependence indicates the V${}_{\mathrm{B}}^{-2}$ and V${}_{\mathrm{B}}^{-3}$ defect charges are very much localized. Indeed, their defect formation energies are satisfactorily corrected with the PC corrections as well.

Figure 2 also shows the rest of the VTLs involving vacancies in the four materials as a function of the supercell size. Again, our PC corrections work well with reasonable accuracy on average. Regarding cation vacancies [Figs. 2(e) and 2(g)], the corrected values are almost constant within 0.1 eV even at the smallest supercells, as in the case of VB in c-BN [Figs. 2(a) and 2(b)]. In contrast, the VTLs of anion vacancies are not so well corrected especially at the smallest supercells [Figs. 2(c), 2(d), 2(f), and 2(h)–2(j)]: the worst case is VN($+2$/$+3$) in c-BN with a deviation of 0.8 eV at the 64-atom supercell compared with the corrected VTL at the 1000-atom supercell [Fig. 2(c)]. This means the PC corrections are not appropriate for the anion vacancies, which is consistent with the fact that their total energies are not corrected by the PC corrections (see Fig. S1 in the Supplemental Material Not ).

As shown in Fig. 2, corrections with the alignment term successfully estimate the VTLs in the dilute limit even with relatively smaller supercells and work better than the PC corrections in most cases. Notable examples are VN($+2$/$+3$) and VN($+3$/$+2$) in c-BN [Figs. 2(c) and 2(d)]; the PC corrections lead to errors of 0.8 and 0.7 eV at the smallest supercells, respectively, while the errors are reduced to 0.0 eV when the alignment terms are incorporated. Exceptions are VB($-2$/$-3$) in c-BN [Fig. 2(a)] and VGa($-2$/$-3$) in GaN [Fig. 2(e)], where the alignment terms slightly increase the errors.

To see the overall performance of our correction methods, we align the VTLs calculated with the second-smallest supercells (128 atoms for GaN and 216 atoms for the others) to the same scale in Fig. 3; such moderate size supercells are regularly adopted for practical defect calculations nowadays. Although the VTLs are reasonably well corrected even at the PC correction level, the addition of the alignment term reduces the errors to be within 0.12 eV in all the cases.

Another important merit to consider for the alignmentlike term is that it removes the projector augmented wave (PAW) or pseudopotential core radius dependences from the VTLs. Bruneval and co-workers have reported that the charged defect formation energies depend on the PAW or pseudopotential radii, because the standard of the electrostatic potential depends on the radii, and this modifies the total energy of a charged system Bruneval et al. (2014). The VTLs depend on these radii as well because the charge state is different by $+1$ or $-1$ between the final and initial states, but such an energy dependence is removed when an alignmentlike term, more specifically the second term in the parentheses in Eq. (13), is considered.

Finally, we emphasize our correction method is adaptive for evaluating any quantities involving different charge states at fixed geometry. One of the most important examples is the quantity relevant to the generalized Koopmans’ theorem (gKT) Perdew et al. (1982). On the basis of the gKT, the deviation from linearity is evaluated at a fixed geometry as $\Delta_{\mathrm{XC}}$:

[TABLE]

where $\varepsilon[Q;S(Q)]$ is a single-particle level of a localized orbital such as a defect-induced state. $\Delta_{\mathrm{XC}}$ has recently been actively used to evaluate the accuracy of the exchange-correlation functional for studying point defects and small polarons Miceli et al. (2018); Elmaslmane et al. (2018); Gake et al. (2019). The finite size correction to the eigenvalue for the second term in Eq. (14) has been well established by many authors Jain et al. (2011); Komsa et al. (2012); Chen and Pasquarello (2013). Conversely, the correction to the first bracketed term has not been discussed so far. This is indeed the same as Eq. (2), and therefore our newly developed correction is applicable for it. We note here that the eigenvalues are not corrected as accurately as the total energies or the VTLs as shown in Fig. S2 in the Supplemental Material Not . Therefore, we recommend the use of larger supercells when one needs to accurately determine $\Delta_{\mathrm{XC}}$.

In summary, we derived correction schemes for calculating VTLs using defect supercells, and tested the effectiveness with ten vacancies in prototypical materials: cubic-BN, GaN, MgO, and 3C-SiC. At the smallest supercells, the PC corrections satisfactorily evaluate the VTLs in the dilute limit involving cation vacancies but do not work well for those involving anion vacancies. We then verified that potential alignment terms reduce the absolute errors to less than 0.12 eV at the second-smallest supercells in our test set. Since our method is routinely and automatically applied, systematic calculations of VTLs are feasible. Furthermore, it is adaptive for evaluating any quantities involving different charge states at fixed geometry, as represented by the gKT-relevant quantity.

Computational methods. First-principles calculations were performed using the PAW method Blöchl (1994), as implemented in VASP Kresse and Furthmüller (1996); Kresse and Joubert (1999). We used PAW data sets with radial cutoffs of 0.91, 1.40, 1.52, 1.31, 0.74, 0.90, and 0.86 Å for B, Ga, Mg, Si, N, O, and C, respectively. In the supercell calculations, a 2 $\times$ 2 $\times$ 2 k-point mesh was applied to the 216-atom or smaller supercells, while $\Gamma$-only sampling was applied to the others, and spin polarization was considered in all the defect calculations. Note that all the VTLs considered did not involve partial occupation at the VBM, CBM or defect localized states.

Acknowledgments. Fruitful discussions with Naoki Tsunoda are deeply appreciated. This work was supported by the Grants-in-Aid for Scientific Research B (Grant No. 19H02416), PRESTO (Grant No. JPMJPR16N4), CREST (Grant No. JPMJCR17J2), and the Support Program for Starting Up Innovation Hub MI2I from JST, Japan. The computing resources of ACCMS at Kyoto University and Research Institute for Information Technology at Kyushu University were used for part of this work.

T.G. and Y.K. contributed equally to this work.