Calculations of Al dopant in α-quartz using a variational implementation of the Perdew–Zunger self-interaction correction

While Kohn–Sham density functional theory (KS-DFT) has been extremely successful and is now widely used in studies of molecules and condensed matter [1, 2], there are several shortcomings of this approach of which the description of localized electronic states in semiconductors and insulators is particularly problematic [3]. An important issue is the energy of the defect state with respect to the valence and conduction band edges as well as structural changes due to the introduction of the defect. A widely studied and well documented case is the Al-atom substitutional defect in α-quartz. Since the Al-atom has three valence electrons while the Si-atom it replaces has four, an electron hole is introduced in the electronic structure.

Extensive electron paramagnetic resonance (EPR) experiments have been interpreted in terms of a localized hole on one of the O-atoms neighboring the Al-atom and a 12% lengthening of the corresponding Al–O bond [4, 5]. An experimental estimate of the defect energy level has been obtained from optical absorption measurements where an electron is excited from a lone pair of an O-atom adjacent to the Al-atom into the localized defect state, a light-induced transfer of the hole between oxygen atoms of the AlO₄ tetrahedron. The estimate obtained for the energy difference between the defect level and valence band edge is 1.96 eV [16]. The band gap of α-quartz is estimated experimentally to be 9 eV so the defect level is much closer to the valence band edge than the conduction band edge. Both the energetics as well as atomic structure related to the Al substitutional defect in α-quartz are, therefore, well established from experimental measurements.

Theoretical calculations of clusters using the Hartree–Fock (HF) approximation have supported the interpretation of the EPR data while DFT calculations have predicted a more delocalized hole. Pacchioni et al [5] have presented a review of calculations of the atomic structure using various theoretical methods as well as systematic computational studies of finite clusters of varying size. Comparison of the calculated hyperfine splitting with the EPR data has shown conclusively that the original interpretation of the EPR measurements was correct while the DFT calculations were in error. HF calculations and DFT calculations based on the generalized gradient approximation (GGA) such as the PBE and PW91 functionals [7] tend to have errors of opposite sign. For example, HF calculations predict the band gap of α-quartz to be 17 eV [8], while the energy difference between conduction and valence band edges is 5.9 eV in results obtained using the PBE functional [9]. By mixing different theoretical approaches that have errors of opposite sign to form a hybrid functional, the errors can be made to cancel to some extent and a number of such mixtures have been proposed for various different purposes. Calculations of small clusters using the B3LYP functional which includes 20% exact exchange have given a delocalized hole spread over two O-atoms, while calculations using the BB1K functional which includes 42% exact exchange have given localization on a single O-atom [6]. The construction of hybrid functionals, however, does not address the root of the problem inherent in KS-DFT functionals.

The basic reason for the instability of the localized hole is the self-interaction error in commonly used KS-DFT functionals. Lægsgaard and Stokbro illustrated how the self-interaction destabilizes the localized hole using a model calculation [10]. Even though Perdew and Zunger proposed a method for self-interaction correction (SIC) more than 30 years ago [11], there have been few attempts to make use of this approach in calculations of defect states in solids. Several perfect crystals have been studied with a linear-muffin-tin-orbital implementation of SIC and the local density approximation (LDA) functional, were the correction was applied selectively to d- or f-electrons, and good results in terms of band gaps, f- and d-state occupation and magnetic properties obtained [12–14]. The full implementation of SIC, including structural optimization, is challenging in many respects, as described in the following section. A limited SIC procedure has been applied to Al-doped α-quartz where the correction was only applied to the spin density, i.e. a single unpaired electron, while the rest of the electrons were described with the PBE functional in a spin restricted formulation [15]. Such a limited version of the Perdew–Zunger correction was shown to help stabilize the localized hole to some extent, but made it only metastable. The valence band energy is still too high in such calculations because of the remaining self-interaction error in the rest of the occupied orbitals. The method cannot give an accurate estimate of the energy of the defect state with respect to the valence band edge since the two are treated at a different level of theory.

A semi-empirical DFT + U approach has been used to estimate the energy of the defect state. There, a penalty is added to inhibit fractional occupation of atomic orbitals and, thereby, stabilizing localized states on atoms [17]. By tuning the U parameter to 7.0 eV, a localized hole was produced 1.1 eV above the valence band. A larger value of the U parameter would presumably shift the defect state even higher in energy. But, a significant correction is not made to the energy of the valence band and the band gap obtained from DFT + U calculations using a U value of 7.0 eV is only 6.5 eV [9].

We present here results of variational, self-consistent calculation using the Perdew–Zunger SIC applied equally to all the valence electrons in the system. The electronic structure as well as the atomic structure is analyzed and found to be in close agreement with the experimental estimates, while calculations with commonly used GGA (PBE, PW91) and hybrid (HSE) functionals are found to give results that are not consistent with the experimental results. In the following section, section 2, the SIC and its implementation are discussed. Section 3 describes the simulated system and in section 4 the results are presented. The article concludes with a discussion and summary in section 5.

In KS-DFT [1], the energy of the ground state of an electronic system is estimated as a functional of the total electron density for each spin, ${\rho }_{\uparrow }({\bf{r}})$ and ${\rho }_{\downarrow }({\bf{r}})$, as

$$\begin{eqnarray}&&{E}^{\mathrm{KS}}[{\rho }_{\uparrow },{\rho }_{\downarrow }]={T}_{s}[{\rho }_{\uparrow },{\rho }_{\downarrow }]+{V}_{\mathrm{ext}}[\rho ]+{E}_{{\rm{C}}}^{\mathrm{KS}}[\rho ]+{E}_{\mathrm{xc}}[{\rho }_{\uparrow },{\rho }_{\downarrow }],\end{eqnarray} \tag{ 1 }$$

where $\rho ={\rho }_{\uparrow }+{\rho }_{\downarrow }$. Here, T_s is the kinetic energy of an independent electron reference system having the same electron density as the real system, V_ext is the external potential representing the interaction of the electrons with the nuclei, ${E}_{{\rm{C}}}^{\mathrm{KS}}$ is an estimate of the Coulomb repulsion between the electrons, and E_xc is the exchange and correlation energy. While a set of orthonormal orbitals is introduced in order to obtain a better estimate of the kinetic energy, the rest of the terms in the functional involve just the total electron density of each spin component. The total energy only depends on the space spanned by the orbitals i.e. the energy is invariant to unitary transformations of the orbitals. A mathematically convenient choice of orbitals can be made so as to diagonalize the orthonormalization constraint matrix and this set of orbitals is referred to as the Kohn–Sham orbitals.

The electron–electron Coulomb repulsion is estimated in KS-DFT from the total electron density as

$$\begin{eqnarray}&&{E}_{{\rm{C}}}^{\mathrm{KS}}[\rho ]=\displaystyle \frac{1}{2}\iint {{\rm{d}}}^{3}{\bf{r}}{{\rm{d}}}^{3}{{\bf{r}}}^{\prime }\;\displaystyle \frac{\rho ({\bf{r}})\rho ({{\bf{r}}}^{\prime })}{| {\bf{r}}-{{\bf{r}}}^{\prime }| }.\end{eqnarray} \tag{ 2 }$$

Written in terms of orbital densities, the total electron density is $\rho ({\bf{r}})={\displaystyle \sum }_{{}_{m\sigma }}{\rho }_{m\sigma }({\bf{r}})$, where ${\rho }_{m\sigma }({\bf{r}})=| {\varphi }_{m\sigma }({\bf{r}}){| }^{2}$ is the density associated with orbital m with spin σ, and ${E}_{{\rm{C}}}^{\mathrm{KS}}$ can be expressed as

$$\begin{eqnarray}&&{E}_{{\rm{C}}}^{\mathrm{KS}}=\displaystyle \frac{1}{2}\iint {{\rm{d}}}^{3}{\bf{r}}{{\rm{d}}}^{3}{{\bf{r}}}^{\prime }\displaystyle \sum _{m\sigma }\displaystyle \sum _{{m}^{\prime }{\sigma }^{\prime }}\displaystyle \frac{{\rho }_{m\sigma }({\bf{r}}){\rho }_{{m}^{\prime }{\sigma }^{\prime }}({{\bf{r}}}^{\prime })}{| {\bf{r}}-{{\bf{r}}}^{\prime }| },\end{eqnarray} \tag{ 3 }$$

where the double sum includes the diagonal, self-interaction terms

$$\begin{eqnarray}&&{E}_{{\rm{C}}}^{\mathrm{SI}}=\displaystyle \frac{1}{2}\iint {{\rm{d}}}^{3}{\bf{r}}{{\rm{d}}}^{3}{{\bf{r}}}^{\prime }\displaystyle \sum _{m\sigma }\displaystyle \frac{{\rho }_{m\sigma }({\bf{r}}){\rho }_{m\sigma }({{\bf{r}}}^{\prime })}{| {\bf{r}}-{{\bf{r}}}^{\prime }| }.\end{eqnarray} \tag{ 4 }$$

These diagonal terms represent Coulomb interaction of the orbital densities with themselves. This is the Coulomb part of the so-called self-interaction error in KS-DFT. The term ${E}_{{\rm{C}}}^{\mathrm{KS}}$ is often referred to as the 'Hartree energy', even though Hartree in his early calculations of electronic structure of atoms did not include the diagonal, self-interaction terms [18].

One of the tasks of the ${E}_{\mathrm{xc}}[{\rho }_{\uparrow },{\rho }_{\downarrow }]$ term in the KS-DFT functional is to cancel out these unphysical self-interaction terms in the estimate of the electron–electron Coulomb energy. Such a cancellation does occur in HF theory where exchange is calculated exactly. However, functionals developed for KS-DFT typically only give partial cancellation and, therefore, include a self-interaction error which is the source of several inaccuracies in results obtained from such calculations, in particular the destabilization of localized electronic states. The more localized an orbital is, the larger the repulsive self-interaction Coulomb energy is and, thereby, an artificial tendency towards delocalization.

Perdew and Zunger proposed a procedure where a KS-DFT functional (originally the LDA) is corrected by explicitly subtracting an orbital based estimate of the self-interaction [11]. The corrected energy functional, E^SIC, is

$$\begin{eqnarray}&&{E}^{\mathrm{SIC}}[\{{\rho }_{m\sigma }\}]={E}^{\mathrm{KS}}[{\rho }_{\uparrow },{\rho }_{\downarrow }]-\displaystyle \sum _{m\sigma }\left({E}_{{\rm{C}}}^{\mathrm{KS}}[{\rho }_{m\sigma }]+{E}_{\mathrm{xc}}[{\rho }_{m\sigma },0]\right),\end{eqnarray} \tag{ 5 }$$

where the spin density for spin component σ is ${\rho }_{\sigma }={\sum }_{m}{\rho }_{m\sigma }$. This correction procedure cancels out the diagonal terms in ${E}_{{\rm{C}}}^{\mathrm{KS}}$ and gives a self-interaction free functional for one-electron systems where E_xc vanishes, but for many-electron systems this orbital based estimate of the self-interaction in ${E}_{\mathrm{xc}}[{\rho }_{\uparrow },{\rho }_{\downarrow }]$ is only approximate. Unlike a KS-DFT functional, the total energy given by E^SIC depends on the orbital densities, ${\rho }_{m\sigma }$. Because of this, the total energy is not unitary invariant and the variational minimization of the energy is mathematically more challenging than for KS-DFT. No longer can a convenient linear combination be chosen for the orbitals and a new inner loop needs to be added to the self-consistency procedure to find the set of optimal orbitals. The variational minimization is, therefore, more challenging. As a result, most calculations involving some form of SIC have not been variational and in many cases not even self-consistent. Rather, the PZ-SIC has mostly been applied only to selected electrons (d- and f-electrons, or only unpaired electron), and in many cases assuming a priori some form for the optimal orbitals (Foster-Boys or Wannier localized orbitals), applied perturbatively without correcting the electron density, etc For a recent review of PZ-SIC calculations, see [19].

The present calculations make use of a variational, self-consistent implementation of PZ-SIC applied equally to all occupied orbitals. The method has been described elsewhere [21, 22] and will only be reviewed here briefly. Two orthonormal basis sets are used in the variational minimization of the energy. The first set consists of the optimal orbitals $\{{\varphi }_{m\sigma }\}$ that minimize the self-interaction corrected energy functional, equation (5), and define the SIC potential

$$\begin{eqnarray}&&{v}_{m\sigma }^{\mathrm{SIC}}({\bf{r}})=-{v}_{{\rm{C}}}[{\rho }_{m\sigma }]-{v}_{\mathrm{xc}}[{\rho }_{m\sigma }]\end{eqnarray} \tag{ 6 }$$

which enters the Hamiltonian as shown below. The second set consists of canonical orbitals ${\psi }_{n}({\bf{r}})$ that are eigenfunctions of the Hamiltonian operator

$$\begin{eqnarray}&&{\hat{H}}_{n\sigma }={\hat{H}}^{0}+{\hat{V}}_{n\sigma }\end{eqnarray} \tag{ 7 }$$

which is now comprised of a part ${\hat{{\text{}}H}}^{0}$ that is unitary invariant and an orbital density dependent term ${\hat{V}}_{n\sigma }$. The ${\hat{{\text{}}H}}^{0}$ Hamiltonian is of KS-DFT form

$$\begin{eqnarray}&&{\hat{H}}^{0}=-\displaystyle \frac{1}{2}{\nabla }^{2}+{v}_{\mathrm{ext}}+{v}_{{\rm{C}}}[\rho ]+{v}_{\mathrm{xc}}[{\rho }_{\uparrow },{\rho }_{\downarrow }],\end{eqnarray} \tag{ 8 }$$

and depends only on the total electron density for each spin component $[{\rho }_{\uparrow },{\rho }_{\downarrow }]$, while the orbital density dependent part ${\hat{V}}_{n\sigma }$ is defined by the SIC potential ${v}_{m\sigma }^{\mathrm{SIC}}({\bf{r}})$ of the optimal orbitals as

$$\begin{eqnarray}&&{\hat{V}}_{n\sigma }{\psi }_{n\sigma }({\bf{r}})=\displaystyle \sum _{m}{W}_{{nm}\sigma }{v}_{m\sigma }^{\mathrm{SIC}}({\bf{r}}){\varphi }_{m\sigma }({\bf{r}}).\end{eqnarray} \tag{ 9 }$$

The matrix, ${\bf{W}}$, gives the unitary transformation relating the two basis sets

$$\begin{eqnarray}&&{\varphi }_{m}=\displaystyle \sum _{n}{W}_{{nm}}{\psi }_{n}.\end{eqnarray} \tag{ 10 }$$

In order to minimize E^SIC, an iterative procedure is followed where the eigenvalue problem

$$\begin{eqnarray}&&{\hat{H}}_{n\sigma }{\psi }_{n\sigma }={\epsilon }_{n\sigma }{\psi }_{n\sigma }\end{eqnarray} \tag{ 11 }$$

is solved and a unitary optimization [21] carried out to find the ${\bf{W}}$ that minimizes the energy. It has recently been realized that it is important to allow the orbitals be complex valued functions [23, 24].

Previous calculations of molecules and solids have shown that the accuracy of E^SIC depends strongly on the KS-DFT functional it is applied to. While the correction is too small when applied to LDA, it is too large when applied to the PBE [7] and PW91 [25] functionals, as judged from atomization energy of molecules and band gaps [24, 26]. Ideally, a new exchange enhancement factor in the generalized gradient functional form should be developed for use with PS-SIC, that has a weaker enhancement factor than PBE or PW91 [27]. Here, a scaling with a factor of 1/2 will be used as has been justified based on the adiabatic connection formula and has been shown to work well for a wide range of properties and various different types of systems [24, 26]

$$\begin{eqnarray}&&{E}_{s}^{\mathrm{SIC}}[\{{\rho }_{m\sigma }\}]={E}^{\mathrm{PW}91}[{\rho }_{\uparrow },{\rho }_{\downarrow }]-\displaystyle \frac{1}{2}\displaystyle \sum _{m\sigma }\left({E}_{{\rm{C}}}^{\mathrm{KS}}[{\rho }_{m\sigma }]+{E}_{\mathrm{xc}}[{\rho }_{m\sigma },0]\right).\end{eqnarray} \tag{ 12 }$$

We choose to use the PW91 functional rather than the PBE functional because the exchange enhancement factor goes to zero for large values of the reduced gradient. Since the corrected energy functional depends on orbital densities which have nodal planes, the reduced gradient of the orbital densities can become large even if the total electron density only has a small reduced gradient.

While variational, self-consistent calculations of ${E}_{s}^{\mathrm{SIC}}$ using complex optimal orbitals have been used to study molecules and even clusters of molecules [24, 28, 29], the present work represents the first application of this implementation to a defect in an extended solid.

The calculations were performed using a 72 atom cell subject to periodic boundary conditions. Only the gamma point for this extended cell was included in the k-point sampling. For the perfect α-SiO₂ crystal, optimization of the lattice constants using the SIC-PW91 functional gave a = b = 4.972 Å and c = 5.473 Å which is in good agreement with measurements, the former being 1.19% and the latter 1.08% larger than experimental values [30]. Without the SIC, the PW91 functional gives lattice constants that are 1.43% and 2.36% larger than the experiment. The calculations were performed with a real space grid with a spacing of 0.2 Å, and an orthogonal cell of size 9.94 Å × 8.61 Å × 10.93 Å using the GPAW software [31]. The cell size was not changed after substituting one of the Si-atoms with an Al-atom.

The PAW method was used to represent the ion core [32]. The optimal orbitals are then written as

$$\begin{eqnarray}&&{\varphi }_{m}({\bf{r}})={\tilde{\varphi }}_{m}({\bf{r}})+\displaystyle \sum _{a}\displaystyle \sum _{i}\left({\phi }_{i}^{a}({\bf{r}})-{\tilde{\phi }}_{i}^{a}({\bf{r}})\right)\langle {\tilde{p}}_{i}^{a}| {\tilde{\varphi }}_{m}\rangle ,\end{eqnarray} \tag{ 13 }$$

where the smooth projector function $| {\tilde{p}}_{i}^{a}\rangle$ determines the expansion coefficients for the orbitals inside an augmentation sphere around each atom a. The smooth optimal orbital ${\tilde{\varphi }}_{m}({\bf{r}})$ is given by a rotation of the smooth canonical orbitals through the transformation ${\bf{W}}$, and the same goes for the projections $\langle {\tilde{p}}_{i}^{a}| {\tilde{\varphi }}_{m}\rangle$. A more detailed description of the PAW implementation for SIC calculations is given elsewhere [33].

The PZ-SIC calculations were started by constructing first highly localized orbitals. This is accomplished by including initially only the Coulomb term in the SIC. Then, the contribution of E_xc is gradually increased. The reason is that highly delocalized orbitals can also represent a local minimum for the E^SIC functional when applied to periodic systems. In order to converge on the global minimum corresponding to localized orbitals, the iterative self-consistency calculation needs to be started with sufficiently localized orbitals [33].

For comparison with the PW91 and SIC-PW91 results, calculations using the HSE hybrid functional [34] were carried out using the VASP software [35].

The band gap estimated from the calculated eigenvalues of the perfect α-SiO₂ crystal is 8.87 eV, in good agreement with the experimental measurements [36]. The gap increases when the SIC is applied mainly because of a lowering of the valence band edge. The energy of valence band orbitals is too high in KS-DFT because of the repulsive self-interaction error. Similar accuracy of band gap estimates has been obtained in SIC calculations for several semiconductors and insulators [26].

The properties of an Al-doped crystal were calculated by replacing one of the Si-atoms with an Al-atom and minimizing the energy given by the ${E}_{s}^{\mathrm{SIC}}$ functional with respect to both orbitals and atomic structure. The resulting atomic structure and spin density are shown in figure 1. The spin density, ${\rho }_{\uparrow }-{\rho }_{\downarrow }$, is highly localized on a single O-atom and the distance between that atom and the Al-atom becomes 11% longer than the short Al–O distances, as listed in table 1. This is in good agreement with the EPR experiments from which the lengthening has been estimated to be 12% [4, 5]. Embedded HF calculations using large clusters have given slightly larger value, 14% [5]. Bader analysis of the spin density [37, 38] gives integrated spin of 0.77 unpaired electrons on the O-atom, showing that the hole is mainly localized on this single atom. The shape of the spin density shows that the hole largely corresponds to a p-orbital. As has been discussed in detail previously, GGA functionals give a delocalized hole spread over all four O-atoms neighboring the Al-atom, the PW91 functional being no exception, as shown in figure 1 and listed in table 1.

Zoom In Zoom Out Reset image size

The calculated DOS is shown in figure 2. After applying some broadening, the results of the PW91 calculation do not show a distinct peak corresponding to the defect state in the band gap. Furthermore, the energy difference between the valence and conduction band edges is only 5.7 eV, much smaller than the experimentally estimated band gap. When PZ-SIC is applied, however, a band gap of 8.87 eV is obtained, in good agreement with measured value of 9 eV [36]. The DOS is quite similar in shape to that obtained with GW calculations [39], where the band gap was estimated to be 9.4 eV. In the SIC calculation, a peak corresponding to the defect state appears 2.2 eV above the valence band edge. This is in good correspondence with the optical absorption experiments which gave an estimate of 1.96 eV [16]. Without the SIC, the valence band is too high in energy and the defect state calculated with the PW91 functional is only 0.17 eV above the valence band edge, as shown in figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows the calculated electronic energy levels using the HSE hybrid functional, which contains 20% exact exchange, as well as the PW91 and PZ-SIC results. The inclusion of exact exchange both lowers the valence band edge and raises the conduction band edge as compared with the PW91 results, but the band gap obtained, 7.61 eV, is still significantly smaller than the experimental value. Furthermore, the defect state is only 0.93 eV over the valence band edge, significantly lower than the experimental estimate. The spin density obtained from the HSE functional is delocalized over two of the O-atoms neighboring the Al-atom. The corresponding Al–O bonds are 1.72 and 1.80 Å while the remaining two are 1.69 Å (see table 1). This is quite consistent with results of previously reported calculations of small clusters using the B3LYP hybrid functional [10]. In order to increase the band gap sufficiently and move the defect energy level high enough above the valence band edge, a larger fraction of exact exchange would need to be included in order to get agreement with the experimental results on this system. For other systems, 20% exact exchange can work well. This mixing factor is increasingly being treated as a free parameter in hybrid calculations, making the approach semi-empirical.

The calculations presented here, as well as other calculations using the same approach [24, 26], illustrate that an orbital density dependent form of the energy functional, of which the PZ-SIC is an example, can significantly improve the accuracy of density functional calculations. The electron hole introduced by the Al-dopant atom is properly localized, the band gap is estimated accurately from the orbital energies and the energy of the defect state above the valence band edge is consistent with optical absorption experiments. Elimination of the self-interaction terms in the electron–electron Coulomb interaction is clearly an important improvement over the KS-DFT functional form. Interestingly, the early calculations of Hartree [18] also excluded these self-interaction terms, but the interaction was calculated using canonical orbitals. In the variational, self-consistent implementation of PZ-SIC used here, the Coulomb interaction between the electrons is estimated using optimal orbitals, which in the present case resemble sp³ hybrid orbitals rather than the more delocalized canonical orbitals. The optimization of the orbital densities introduces additional computational effort, but has the advantage that meaningful, energy optimal orbitals are obtained from the calculation. The scaling of the computational effort of PZ-SIC calculations is the same as for GGA functional calculations, i.e. it scales as the system size cubed, but there is a large prefactor due to the fact that an effective potential needs to be evaluated for each orbital. These calculations can, however, readily be carried out in parallel and the wall clock time can, thereby, be reduced quite easily by using a larger number of processors.

The exchange and correlation term of the PW91 functional is, however, not optimal for PZ-SIC. A smaller exchange enhancement factor can be used [27], bringing the exchange estimate closer to LDA, as illustrated by the fact that PZ-SIC gives an overcorrection when applied to PBE while it is an undercorrection when applied to LDA [24]. By developing an optimal exchange enhancement factor tailored to PZ-SIC, the scaling factor of 1/2 used here as well as in previous calculations based on the PBE functional [24, 26] could be dropped. The promising results that have been obtained here with the variational, self-consistent implementation of PZ-SIC show that this orbital density dependent functional form can give highly accurate results for both energetics as well as structure, the first time this is achieved for the Al–SiO₂ system.

This work was supported by the Icelandic Research Fund, Nordic Energy Research and the Academy of Finland (grants 263294 and 278260). The calculations were carried on the Nordic High Performance Computer (NHPC) in Iceland.