Can we rely on hybrid-DFT energies of solid-state problems with local-DFT geometries?

TL;DR

This paper assesses the reliability of hybrid-DFT energies in solid-state calculations when using geometries pre-relaxed with local-DFT, highlighting key sources of error and their impact on energy accuracy.

Contribution

It identifies and analyzes the main error sources in single-point hybrid-DFT calculations based on local-DFT geometries for solids, providing insights into their effects.

Findings

Strain-related errors are minimal (~10 meV) due to error cancelation.

State over-mixing causes significant errors without cancelation.

Pre-relaxed geometries are generally reliable for energy difference calculations.

Abstract

Hybrid functionals often improve considerably the accuracy of density-functional calculations, in particular of quantities resulting from the band structure. In plane-wave (PW) calculations this benefit comes at the cost of an increase in computation time by several orders of magnitude. For this reason, large-scale problems addressed within the PW formalism have to rely on pre-relaxed atomistic geometries, obtained with cheaper local or semi-local exchange-correlation functionals. We investigate how suitable these geometries are when plugged into single-point hybrid-DFT calculations. Based on several case studies, we find two important sources of error originating from (i) bond strain and (ii) over-mixing between defect and crystalline states. The first arises from the mismatch between the pre-relaxed geometry and that obtained after a subsequent hybrid-DFT-level relaxation. The second occurs when defect states edging an underestimated band gap artificially mix with crystalline states, affecting the local bonding character of the defect, and therefore leading the spurious hybrid-DFT energies. Due to cancelation effects, the lingering strain contributes little ($\lesssim10$ meV) to the error bar of quantities based on energy differences of pre-relaxed structures. The error from state over-mixing does not benefit from cancelation effects and has to be monitored with caution.

Tables

Table 1. Table 1: Four data sets related to (1) V O in MgO in different charge states, (2) four configurations of O i in Si, (3) two reconstructions for the Si(001) surface, and (4) two configurations of C i in 3 ​ C 3 𝐶 3C -SiC. First and second data rows of each set show displacement magnitudes ( d 𝑑 d ) of selected atoms after full HSE relaxation, relatively to pre-relaxed positions. The third row of each data set reports the relaxation energy Δ ​ E Δ 𝐸 \Delta E or surface relaxation energy Δ ​ σ = Δ ​ E / 2 ​ A Δ 𝜎 Δ 𝐸 2 𝐴 \Delta\sigma=\Delta E/2A of pre-relaxed structures, where A 𝐴 A is the surface unit cell area. The fourth and fifth rows of each data set show HSE results using pre-relaxed and fully-relaxed geometries. These include transition levels E ​ ( q / q + 1 ) − E c 𝐸 𝑞 𝑞 1 subscript 𝐸 c E(q/q\!+\!1)-E_{\mathrm{c}} with respect to the conduction band bottom, relative energies E − E GS 𝐸 subscript 𝐸 GS E-E_{\mathrm{GS}} with respect to the ground state, and surface formation energies, σ 𝜎 \sigma (see text).

MgO:VO	Units	$q=0$	$q=+1$	$q=+2$
$d$(Mg)	Å	0.007	0.010	0.015
$d$(O)	Å	0.001	0.001	0.001
$\Delta E$	eV	$-0.003$	$-0.012$	$-0.038$
$E_{\text{pre}}\left(q/q+1\right)-E_c$	eV	2.872	4.190	—
$E_{\mathrm{full}}\left(q/q+1\right)-E_c$	eV	2.863	4.163	—
Si:Oi	Units	$C_{1h}$	$C_2$	$D_{3d}$	$C_{2v}$
$d$(O)	Å	0.129	0.127	0.000	0.045
$d$(Si)	Å	0.012	0.012	0.014	0.027
$\Delta E$	eV	$-0.019$	$-0.018$	$-0.009$	$-0.013$
$E_{\mathrm{pre}}-E_{\mathrm{pre},\mathrm{GS}}$	eV	0.000	0.002	0.001	2.758
$E_{\mathrm{full}}-E_{\mathrm{full},\mathrm{GS}}$	eV	0.000	0.002	0.010	2.764
Si(001)	Units	$\left(2\times 1\right)$	b$\left(2\times 1\right)$
$d$(Si1)	Å	0.147	0.115
$d$(Si2)	Å	0.147	0.164
$\Delta \sigma$	meV/Å2	$-1.048$	$-0.856$
${\sigma }_{\mathrm{pre}}$	meV/Å2	97.23	95.39
${\sigma }_{\mathrm{full}}$	meV/Å2	96.18	94.53
$3C$-SiC:Ci	Units	$q=0$	$q=+1$
$d$(C1)	Å	0.013	0.000
$d$(C2)	Å	0.539	0.000
$\Delta E$	eV	$-0.101$	$-0.005$
$E_{\text{pre}}\left(q/q+1\right)-E_c$	eV	0.367	—
$E_{\mathrm{full}}\left(q/q+1\right)-E_c$	eV	0.476	—

Equations

$$E_{\mathrm{f}}(q)=E(q)-\sum _i{n}_i{\mu }_i+q(E_{\mathrm{v}}+E_{\mathrm{F}}).$$

Full text

††: Electronic Structure

Keywords: Density functional theory, Hybrid functionals, Local functionals, Energy, Geometry

1 Introduction

The use of hybrid density functionals [1], which include a portion of Hartree-Fock exchange (HF-exchange), has proven to be an invaluable approach to density functional theory (DFT). Within the solid-state community, the one proposed by Heyd, Scuseria and Ernzerhof [2, 3] (HSE) has become a rather popular choice, particularly for diamagnetic and low-spin systems. HSE-like functionals screen long-ranged and computationally demanding HF-exchange interactions, essentially by introducing a fraction of HF-exchange solely within the short-range part of the functional. Long-range contributions are simply accounted for with a local or semi-local functional, often the one prescribed by Perdew, Burke and Ernzerhof (PBE) to the generalised gradient approximation [4].

Besides being able to account for HF-exchange interactions in metals [5], screened hybrid density functionals excel at describing accurately the band structure of finite-gap crystals [6, 7, 8, 9], a virtue with huge importance for the study of defects (including impurities, dislocations and surfaces) with states in the forbidden gap.

While the calculation of the non-local Fock integral is rather efficient when using a Gaussian basis (see for instance Ref. [10]), this is not the case for plane waves (PW), which are a popular and natural choice, particularly among the solid-state community. Plane wave hybrid-DFT calculations can take thousands of times longer than using (semi-)local functionals. For instance, we found that, using the same number of CPUs, one PBE single-point calculation of the oxygen vacancy in MgO took 28 seconds, against over 9 hours after changing the exchange-correlation treatment to HSE (around 1200 times as long). For this reason, the computation of large-scale problems within PW/hybrid-DFT has to rely on pre-relaxed atomic positions obtained using a cheaper local or semi-local approach.

Many examples of this approach have been reported in the literature, including in the study of magnetic materials, two-dimensional materials and surfaces [11, 12, 13], as well as dopants, impurities and radiation-defects in several materials [14, 15, 16, 17, 18]. Significant differences between fully-relaxed PBE- and pre-relaxed HSE-level results were obtained. For instance, quantities like the location of defect transition levels or the height of migration barriers show considerable discrepancies. This raises the following question: if the fully-relaxed HSE ground state structures were employed instead, would the result be the same? Or at least within an acceptable error bar?

In the 1980s, many self-consistent local density calculations were carried out assuming fixed structures obtained on the basis of known bond lengths and bond strengths (see for instance Ref. [19]). Despite being less refined, this approximation is similar to using DFT/pseudopotential pre-relaxed geometries for the calculation of Mössbauer parameters using all-electron full-potential methods [20]. This approach has also been used to obtain the formation energy and electronic structure of defects using Hedin’s $GW$ method [21, 22] or diffusion quantum Monte Carlo [23]. Again, the above practice raises the question: can we actually rely on energies of state-of-the-art calculations that employ local-DFT geometries?

Let us look at this problem with the help of Figure 1, where potential energy curves corresponding to sequential PBE- and HSE-level relaxations are illustrated. A starting (guessed) structure is subjected to a first relaxation step using the PBE functional. This provides the pre-relaxed minimum-energy geometry, $Q_{\mathrm{pre}}$, which is plugged into a single-point HSE calculation to provide the pre-relaxed state $(Q_{\mathrm{pre}},E_{\mathrm{pre}})$. A subsequent relaxation step also using HSE drives the system to the final ground state configuration $(Q_{\mathrm{full}},E_{\mathrm{full}})$. The use of the pre-relaxed $(Q_{\mathrm{pre}},E_{\mathrm{pre}})$ instead of the fully-relaxed $(Q_{\mathrm{full}},E_{\mathrm{full}})$ state is rather tempting. Unfortunately the underlying effects governing the error $\Delta E=E_{\mathrm{full}}-E_{\mathrm{pre}}$ are not obvious.

Below we find some answers based on the magnitude of spurious strain fields and resonances between defect- and crystalline-related states, due to the change of functional.

We address this issue by assessing the accuracy of HSE single-point calculations performed on pre-relaxed structures. These are compared to corresponding HSE fully-relaxed calculations. Four case studies are considered, namely the oxygen vacancy in magnesium oxide, the oxygen interstitial in silicon, the Si(001) surface, and the carbon interstitial in silicon carbide. In each case, we calculate $\Delta E$ values, the displacement of the most relevant atoms (closest to the defect), as well as selected observables. In Section 2 we disclose the technical details of the calculations. In Section 3 we present and discuss our results. Finally, we draw the conclusions in Section 4.

2 Method details

First-principles calculations were performed within hybrid and semi-local density functional theory using the Vienna Ab-initio Simulation Package (VASP) [24, 25, 26, 27]. The HF-exchange mixing fraction and screening parameter for the HSE functional were $a=1/4$ and $\omega=0.2$ Å*-1*, respectively (resulting in the commonly named HSE06 functional) [3]. The projector-augmented wave method was employed to avoid explicit treatment of core electrons in the Kohn-Sham equations [28]. These were solved self-consistently using the PW formalism, until the total energy between two consecutive steps differed by less than $10^{-7}$ eV. The choice of energy cutoff values for PW, as well as the density of Brillouin-zone (BZ) sampling meshes, essentially depend on the problem at hand (chemical species, material of interest and supercell sizes). They were chosen after convergence tests, ensuring that absolute (relative) energies were converged within less than 10 meV/atom (less than 1 meV/atom).

HSE lattice parameters were adopted in all four case studies referred above. Ionic relaxations stopped when the maximum force acting on every atom became smaller than 0.01 eV/Å. For each case, $(Q_{\mathrm{pre}},E_{\mathrm{pre}})$ and $(Q_{\mathrm{full}},E_{\mathrm{full}})$ pairs were obtained by means of the above two-step recipe (see Figure 1), where $\Delta E=E_{\mathrm{full}}-E_{\mathrm{pre}}$.

For the oxygen vacancy in magnesium oxide (MgO:VO), we used a 64-atom cubic supercell (minus an oxygen atom). The PW energy cutoff was converged at $E_{\mathrm{cut}}=500$ eV and a Monkhorst-Pack (MP) $2\times 2\times 2$ grid of $\mathbf{k}$-points proved adequate [29]. The calculated HSE-level (PBE-level) lattice parameter was $a_{0}^{\mathrm{HSE}}=4.200$ Å ($a_{0}^{\mathrm{PBE}}=4.255$ Å), comparing well with the experimental value of $a_{0}^{\mathrm{exp}}=4.216$ Å [30]. The HSE-level (PBE-level) direct band gap was $E_{\mathrm{g}}^{\mathrm{HSE}}=6.64$ eV ($E_{\mathrm{g}}^{\mathrm{PBE}}=4.73$ eV), underestimating the observed gap of $E_{\mathrm{g}}^{\mathrm{exp}}=7.8$ eV by 15% (40%) [31]. The oxygen vacancy (VO) in MgO is a double donor. We investigated the deviation of the transition levels obtained from single-point energies (of pre-relaxed structures) with respect to those from fully-relaxed HSE-level structures. Unwanted Coulomb interactions between periodic replicas of charged defects were removed from the total energy according to the method proposed by Freysoldt, Neugebauer and Van de Walle [32].

For the oxygen interstitial (Oi) defect in Si, several structures close to the minimum ground state and to the saddle point for migration were explored. Supercells with 64 Si atoms (plus one O atom) were used, and the BZ was sampled using an MP $2\times 2\times 2$ grid of $\mathbf{k}$-points. The calculated lattice parameter was $a_{0}^{\mathrm{HSE}}=5.432$ Å ($a_{0}^{\mathrm{PBE}}=5.469$ Å), also in good agreement with the measured value $a_{0}^{\mathrm{exp}}=5.431$ Å [33]. The Kohn-Sham band gap was estimated as $E_{\mathrm{g}}^{\mathrm{HSE}}=1.15$ eV ($E_{\mathrm{g}}^{\mathrm{PBE}}=0.57$ eV), which is to be compared with $E_{\mathrm{g}}^{\mathrm{exp}}=1.17$ eV from optical experiments [34]. The PW energy cutoff was set at $E_{\mathrm{cut}}=400$ eV.

The third case considered was the Si(001) surface, where the energy difference between symmetric $(2\times 1)$ and asymmetric b$(2\times 1)$ reconstructions was investigated [35]. We employed 19-monolayer-thick symmetric slabs (composed of 38 Si atoms), separated by 11 Å of vacuum space. An MP $6\times 6\times 1$ grid of $\mathbf{k}$-points was used to sample the BZ of both reconstructed surfaces and, like in the previous case, $E_{\mathrm{cut}}=400$ eV.

Finally, carbon interstitial (Ci) defects in 3$C$-SiC were investigated on cubic supercells with 64 atoms (plus the interstitial C atom). The energy cutoff was $E_{\mathrm{cut}}=400$ eV and the BZ sampling was carried out using an MP $4\times 4\times 4$ grid of $\mathbf{k}$-points. The calculated lattice parameter and Kohn-Sham gap were $a_{0}^{\mathrm{HSE}}=4.347$ Å ($a_{0}^{\mathrm{PBE}}=4.380$ Å) and $E_{\mathrm{g}}^{\mathrm{HSE}}=2.25$ eV ($E_{\mathrm{g}}^{\mathrm{PBE}}=1.33$ eV), respectively. Again, these compare with $a_{0}^{\mathrm{exp}}=4.360$ Å and $E_{\mathrm{g}}^{\mathrm{exp}}=2.42$ eV from experiments, as expected [36]. We investigated two structures of the Ci defect reported in the literature, namely a tilted-$\langle 001\rangle$ split-interstitial with $C_{1h}$ symmetry [37] and an upright-$\langle 001\rangle$ split-interstitial with $D_{2d}$ symmetry (or simply $\langle 001\rangle$ split-interstitial) [38]. The latter comprises a $\langle 001\rangle$-aligned C-C dimer at the carbon site (both C atoms are three-fold coordinated and symmetrically equivalent), while in the monoclinic structure the C-C bond makes an angle with the $\langle 100\rangle$ axis, ending up in a structure where one of the C atoms has four-fold coordination, while the other keeps the three-fold coordination.

3 Results and discussion

3.1 Oxygen vacancy in MgO

The oxygen vacancy in MgO has the same symmetry ($O_{h}$) in all stable charge states ($q=0,\,+1$ and $+2$). However, the displacement of its neighbours (from their crystallographic positions) is variable. Accordingly, whereas Mg2+ first neighbours are pushed away from the vacant site by $d\approx 0.085$ Å per ionised electron, O2- next neighbours are attracted to the centre by $d\approx 0.035$ Å for each ionisation. These results are also close to those reported by Rinke and co-workers [39] using the local density approximation, where outward and inward displacements of about 0.09 Å and 0.03 Å per ionisation were reported for Mg and O atoms, respectively.

In the top-most section of Table 1, we report the displacement ($d$) of Mg and O nearest neighbours to VO from fully-relaxed geometries with respect to pre-relaxed ones. Results are shown for neutral, positively and double positively charged defects. Pre-relaxed and fully-relaxed defect geometries are very similar. The largest atomic displacement with respect to the pre-relaxed geometry was 0.015 Å, and that was observed for the Mg first neighbours. Still, after relaxation, total energies differed by at most $\Delta E\sim-40$ meV from those of pre-relaxed structures.

In the above, we inspected the quality of pre-relaxed structures for the calculation of HSE total energies. The same analysis could be done for the calculation of defect formation energies,

[TABLE]

This is the amount of energy required to combine $n_{i}$ elements of species $i$ and form a defective crystal (with computed energy $E$). Elements are assumed to be available in standard phases with respective chemical potentials $\mu_{i}$ (see Ref. [40] and references therein for further details). In Eq. 1, $q$ is the charge state of the defect, obtained by exchanging electrons between an electronic reservoir with chemical potential $E_{\mathrm{v}}+E_{\mathrm{F}}$ (where $E_{\mathrm{F}}$ is the Fermi level with respect to the valence band top energy, $E_{\mathrm{v}}$) and the highest occupied or unoccupied states. It is assumed that the calculation of $\mu_{i}$ and $E_{\mathrm{v}}$ in Eq. 1 does not involve a pre-relaxation step. Hence, the error of the calculated formation energies using pre-relaxed structures is the same as that of total energies, $\Delta E_{\mathrm{f}}\equiv\Delta E$.

We emphasise that $\Delta E\leq 0$ , so that the error involving the energy difference between two pre-relaxed structures benefits from cancelation effects. An example of such a quantity is a defect transition level, which by definition, is given by the location of $E_{\mathrm{F}}$ where $E_{\mathrm{f}}(q)=E_{\mathrm{f}}(q+1)$. In the case of a double donor such as MgO:VO, the $(q+1)$-th donor level (with respect to the conduction band minimum, $E_{\mathrm{c}}$) is given by $E_{\mathrm{c}}-E(q/q+1)=[\epsilon_{\mathrm{c}}+E(q+1)]-E(q)$, where $\epsilon_{\mathrm{c}}$ is the lowest unoccupied state in a bulk calculation. Hence, for MgO:VO we obtain $E_{\mathrm{c}}-E(0/+)=2.86$ eV and $E_{\mathrm{c}}-E(+/\!+\!+)=4.16$ eV using fully-relaxed HSE energies. These figures would be closer to those obtained from $G_{0}W_{0}$ quasi-particle energies [39] should we have accounted for Frank-Condon relaxation effects and for a better correlation treatment to improve the gap width. However, these issues are not relevant for the present analysis. We are interested in assessing the quality of $E_{\mathrm{c}}-E_{\mathrm{pre}}(q/q+1)$ obtained from pre-relaxed energies, with respect to the analogous calculation using fully-relaxed structures, $E_{\mathrm{c}}-E_{\mathrm{full}}(q/q+1)$. The results are shown in Table 1. The error bar of the pre-relaxed results is of the order of few tens of meV, which is quite acceptable for semiconductors and insulators with a gap width in the eV range.

3.2 Interstitial oxygen in silicon

Our second case is a well-established defect in crystalline silicon, namely Si:Oi. In the ground state, Oi adopts a puckered bond-centre configuration [41]. The potential energy surface for rotation of the O atom around the Si-Si bond has the shape of a flat ‘Mexican hat’ (with small bumps in the meV range). Accordingly, the O atom is slightly displaced (0.3 Å) away from the centre of a $[111]$-aligned Si-Si bond, either along $[1\bar{1}0]$ or along $[11\bar{2}]$ (resulting in defects with $C_{2}$ or $C_{1h}$ symmetry, respectively). The perfect bond-centred structure has $D_{3d}$ symmetry and it is a maximum in the potential energy landscape.

Long-range diffusion of Oi in Si occurs via sequential jumps between neighbouring puckered configurations. At the saddle-point, the structure passes close to a $C_{2v}$-symmetric configuration (often referred to as ‘Y-lid’), consisting of a $\langle 100\rangle$-aligned O-Si dimer sharing a Si site, where both O and Si atoms are three-fold coordinated [41]. The relaxation of this structure was achieved with help of force symmetrisation.

The results are shown in Table 1. Concerning HSE-relaxed energies, we obtained the following relative energies with respect to the $C_{1h}$ ground state: $E_{\mathrm{full}}-E_{\mathrm{full},\mathrm{GS}}=2$ meV, 10 meV and 2.76 eV for the $C_{2}$, $D_{3d}$ and $C_{2v}$ structures, respectively, confirming the flatness of the potential around the bond centre. The energy of the $C_{2v}$ structure is consistent with analogous results obtained by Binder and Pasquarello [42], where it was demonstrated that hybrid-DFT was able account for the observed 2.53 eV migration barrier of Oi in silicon [43]. This contrasts with local and semi-local functionals which predict a barrier of about 2 eV [41].

Unlike MgO:VO, pre-relaxed structures of Oi in Si hold a lingering strain when plugged into a HSE calculation. This effect is unavoidable and stems from the fact that PBE over-estimates Si-O bond lengths with respect to those from HSE. The result is the straightening of the Si-O-Si unit after HSE full relaxation, with the O atom approaching the bond centre site.

Due to the soft nature of the Si-O-Si bending potential, the above displacements lead to tolerable energy changes. Values of $\Delta E$ in Table 1 show that fully relaxed HSE total energies are at most a few tens of meV below those obtained from single-point calculations performed on pre-relaxed structures. Interestingly, and again due to cancelation effects, the difference between relative energies whether using fully-relaxed ($E_{\mathrm{full}}-E_{\mathrm{full},\mathrm{GS}}$) or pre-relaxed ($E_{\mathrm{pre}}-E_{\mathrm{pre},\mathrm{GS}}$) structures is negligible (<10 meV).

3.3 Si(001)-$\mathit{(2\times 1)}$ surface

The Si(001) surface reconstruction is known to be made of dimers, each possessing two unsaturated radicals. Experiments show that charge transfer between these radicals leads to dimer buckling [35], where one of the atoms becomes protruded, while the other drops into the surface. The smallest surface unit cell which is able to capture this effect contains a single dimer. We can either calculate a symmetric Si(001)-$(2\times 1)$ reconstruction, where the dimerised Si atoms have the same height (no buckling), or the buckled Si(001)-b$(2\times 1)$.

Table 1 shows that, for both reconstructions, atom displacements of pre-relaxed slabs are large. Large displacements were also found for subsurface Si atoms.

The surface formation energy, $\sigma$, was obtained as $2A\sigma=E-\sum n_{i}\mu_{i}$, where $A$ is the area of the surface unit cell and the factor of two accounts for the identical and opposite-facing surfaces on the slab with total energy $E$. Fully-relaxed surface formation energies are about $|\Delta\sigma|\sim 1$ meV/Å2 lower than pre-relaxed analogues, a value which is regarded too large to make pre-relaxed energies reliable for absolute formation energy calculations.

Such strong relaxation energies and displacements arise from a mismatch between PBE- and HSE-level lattice parameters of about 0.7%, which leads to a corresponding contraction of the whole slab thickness upon full HSE-relaxation.

Despite the above, relative energies benefit from cancelation effects. Table 1 shows that Si(001)-b$(2\times 1)$ is more stable than Si(001)-$(2\times 1)$ by 1.84 meV/Å and 1.65 meV/Å when using pre-relaxed and fully-relaxed energies ($\sigma_{\mathrm{pre}}$ and $\sigma_{\mathrm{full}}$), respectively.

3.4 Carbon self-interstitial in in 3C-SiC

The stable configuration of the carbon self-interstitial in $3C$-SiC has been investigated concurrently by different groups, with different structures being proposed for the neutral charge state. Gali et al. [38] reported a spin-1 split-interstitial with $D_{2d}$ symmetry, being made of a $\langle 001\rangle$-aligned C-C dimer on a carbon site. The electronic structure consisted on a semi-occupied doublet state arising from orthogonal $\pi$-like orbitals on both three-fold coordinated carbon atoms.

Conversely, Bockstedte et al. [37] found a diamagnetic monoclinic structure ($C_{1h}$ symmetry), where the C1-C2 dimer was tilted towards the $\langle 110\rangle$ direction, resulting in three-fold and four-fold coordinated atoms, respectively. We note that, while the paramagnetic state was obtained from a post-corrected hybrid-DFT calculation, the diamagnetic and low-symmetry state was found within PBE.

Our calculations confirm the above conflicting results. PBE- and HSE-relaxations lead to C${}_{\mathrm{i}}(C_{1h})$ and C${}_{\mathrm{i}}(D_{2d})$ ground state structures with total spin $S=0$ and $S=1$, respectively. Therefore, after obtaining the pre-relaxed structure ($C_{1h}$), a subsequent HSE calculation displaces the four-fold coordinated C2 atom by more than 0.5 Å, raising the symmetry to $D_{2d}$ and flipping the spin to $S=1$ (see Table 1). The difference between pre-relaxed and fully-relaxed HSE energies is as much as $\Delta E=-0.1$ eV. Obviously this poses a major problem to the use of the pre-relaxed structures. We will return to this issue below.

In the positive charge state, the Ci defect suffers a weak Jahn-Teller distortion involving a dynamic overlap of the $\pi$-orbitals [37, 38], and the symmetry is lowered from $D_{2d}$ to $D_{2}$. Pre-relaxed and fully-relaxed structures show essentially the same geometry (as shown by the small atomic displacements in the bottom data set of Table 1), and the relaxation energy is only $\Delta E=-5$ meV.

Due to the large relaxation energy obtained for the neutral charge state, the calculation of the donor level using pre-relaxed energies does not profit from error cancelation effects. The result differs from the fully-relaxed donor level by about 0.1 eV (ten times larger than MgO:VO).

We investigated the origin of the structure/spin disparity between the above PBE and hybrid-DFT calculations. In order to rule out dispersion effects due to the finite size of the supercell, we also carried out 216-atom and 512-atom supercell calculations (plus one C atom), using $2\times 2\times 2$ and $1\times 1\times 1$ ($\Gamma$-point) $\mathbf{k}$-point grids for BZ sampling. Full HSE relaxation was not possible for these cells. Instead, the $D_{2d}$ structure was subjected to a symmetry-constrained PBE relaxation, followed by an HSE single-point calculation. Irrespectively of the supercell size, the HSE energy of the pre-relaxed C${}_{\mathrm{i}}(D_{2d},S=1)$ state was lower than that of C${}_{\mathrm{i}}(C_{1h},S=0)$ by 0.1 eV.

We went on and inspected the Kohn-Sham levels of seven neutral diamagnetic structures, obtained after linear interpolation between C${}_{\mathrm{i}}(C_{1h})$ and C${}_{\mathrm{i}}(D_{2d})$. Analysis of defect levels is more conveniently done at the $\mathbf{k}=\Gamma$ point, where the wave-functions are real. However, in order to preserve the BZ-sampling quality, large cubic supercells with 512 atoms (plus one carbon atom) were employed for PBE-relaxations and respective HSE single-point calculations.

The results are shown in Figure 2 by data points connected by Bézier curves for better perception. Circles and squares apply to $C_{1h}$ and $D_{2d}$ defect symmetries, respectively. Also for the sake of clarity, only the highest occupied (closed symbols) and the next four unoccupied levels (open symbols) are shown. Left- and right-hand sides of the figure refer to results obtained with PBE and HSE exchange-correlation functionals, respectively. The valence band top was aligned on both sides at the origin of the energy scale.

The paramagnetic state C${}_{\mathrm{i}}(D_{2d},S=1)$ was also investigated. Its electronic structure is shown in the middle region of the figure (for both PBE and HSE functionals).

If we ignore the difference in the gap width, the electronic structure of C${}_{\mathrm{i}}(C_{1h},S=0)$ is rather similar whether it is calculated using PBE or HSE (left and right edges of Figure 2, respectively). The defect is a symmetric singlet state ($A$ within $C_{1h}$ point group), with a fully occupied p-like deep level localised on C1 and C2 atoms (see inset of 2). As the structure evolves to $D_{2d}$, the $A$-state mixes with the conduction band states and the defect becomes a shallow donor. Note that C${}_{\mathrm{i}}(D_{2d},S=0)$ is not stable – within PBE-level, atomic relaxation drives the geometry back to the $C_{1h}$ structure, while within HSE the ground state is C${}_{\mathrm{i}}(D_{2d},S=1)$.

Hence, comparing PBE and HSE results in Figure 2, we readily conclude that within PBE the exchange interactions are underestimated due to the strong resonance between the doubly degenerate $E$-level and the low-lying conduction band states. For that reason, the PBE functional fails to describe the ground state structure of neutral Ci in $3C$-SiC, and that undermines any pre-relaxed HSE calculation.

4 Conclusions

We investigated the suitability of atomistic geometries, particularly of defects in semiconductors and insulators, obtained within a semi-local DFT method (referred to as pre-relaxed structures), to be used on single-point hybrid-DFT calculations. To that end, four distinct case studies were investigated, namely the oxygen vacancy in magnesium oxide, the oxygen interstitial in silicon, the Si(001) surface and the carbon self-interstitial in cubic silicon carbide.

We found at least two important sources of error that should be monitored. The first is the presence of lingering strain within the pre-relaxed structure, which will be released should a full hybrid-DFT calculation be performed. The relaxation energy, $\Delta E$, arises from slight differences in equilibrium bond lengths as obtained from semi-local and hybrid-DFT methods. It is interpreted as the error bar of single-point HSE energies based on pre-relaxed structures, including of formation energies.

The magnitude of $\Delta E$ is system-dependent. For localised point defects, the effect was estimated to be in the range of a few tens of meV. This is acceptable for the calculation of most defect-related observables, including formation energies. Extended defects, on the other hand, are expected to be affected by larger errors. The surface formation energy of two Si(001) pre-relaxed reconstructions were calculated with a discrepancy of $\Delta\sigma\leavevmode\nobreak\ \sim-1$ meV/Å. This is close to the difference between the surface formation energies of symmetric and buckled dimerised Si(001)-$(2\times 1)$ reconstructions, and larger than the usual error bar needed for this type of calculation. The problem arises from the lattice mismatch generated by the different functionals, that leads to the accumulation of lingering strain across a large volume of the pre-relaxed geometry (particularly in bulk-like regions). The presence of vacuum in the slab allows the strain to relax, releasing relatively large amounts of energy during a full hybrid-DFT relaxation.

We also found that calculations based on energy differences of pre-relaxed structures benefit from error cancelation effects. This is because $\Delta E$ is always negative. The error bar in this case was about one order of magnitude lower than that affecting total energies. This feature applies for instance to transition levels, binding energies, migration barriers, but also to surface formation energy differences.

The second source of error results from an over-mixing of defect states with the host bands during the pre-relaxation stage. This effect also depends strongly on the problem at hand. Local and semi-local functionals are known to underestimate the width of band gaps of insulators and semiconductors. This favours resonances involving defect levels edging conduction band and valence band states. The result is rather similar to the pseudo-Jahn-Teller effect and, as such, leads to artificial bond formation and breaking. Obviously, the spurious pre-relaxed geometries will lead to misleading single-point hybrid-DFT energies. We discuss the effect in the light of a detailed analysis of the carbon interstitial in $3C$-SiC. Many examples which are expected to show analogous resonances have been reported in the literature. These include the negative carbon vacancy in $4H$-SiC [16], or the cadmium vacancy in CdTe [23].

This work is supported by the NATO SPS programme, project number 985215. JC thanks the Fundação para a Ciência e a Tecnologia (FCT) for support under project UID/CTM/50025/2013, co-funded by FEDER funds through the COMPETE 2020 Program. JDG acknowledges the financial support from the I3N through grant BPD-11(5017/2018).