First-principles calculations for point defects in solids

Point defects and impurities strongly affect the physical properties of materials and have a decisive impact on their performance in applications. First-principles calculations have emerged as a powerful approach that complements experiments and can serve as a predictive tool in the identification and characterization of defects. The theoretical modeling of point defects in crystalline materials by means of electronic-structure calculations, with an emphasis on approaches based on density functional theory (DFT), is reviewed. A general thermodynamic formalism is laid down to investigate the physical properties of point defects independent of the materials class (semiconductors, insulators, and metals), indicating how the relevant thermodynamic quantities, such as formation energy, entropy, and excess volume, can be obtained from electronic structure calculations. Practical aspects such as the supercell approach and efficient strategies to extrapolate to the isolated-defect or dilute limit are discussed. Recent advances in tractable approximations to the exchange-correlation functional ($\mathrm{DFT}+U$, hybrid functionals) and approaches beyond DFT are highlighted. These advances have largely removed the long-standing uncertainty of defect formation energies in semiconductors and insulators due to the failure of standard DFT to reproduce band gaps. Two case studies illustrate how such calculations provide new insight into the physics and role of point defects in real materials.

Figure 1

Schematic illustration of formation energy $E^f$ vs Fermi level $E_F$ for an amphoteric defect that can occur in three charge states $q$: $+1$, 0, and $-1$. Solid lines correspond to the formation energy as defined by Eq. (1). The defect exhibits two charge-state transition levels (see Sec. 2d): a deep donor level $ϵ(+/0)$ and a deep acceptor level $ϵ(0/-)$. The thick solid lines indicate the energetically most favorable charge state for a given Fermi level.

Figure 2

Quasiharmonic formation free energy of the vacancy and the $⟨110⟩$ and $⟨111⟩$ split-interstitial dumbbell in bcc iron. Adapted from [211].

Figure 3

Equilibrium concentration of vacancies in aluminum at zero pressure as a function of $T^m/T$, where $T^m$ is the melting temperature. The first-principles results have been obtained according to Eq. (41) with (solid lines) and without (dashed lines) the inclusion of anharmonic (ah) contributions to the free energy of formation. Both LDA and GGA functionals have been applied. Experimental results are included for comparison. Adapted from [115].

Figure 4

DFT-computed temperature dependence of the Gibbs energy of vacancy formation in Al. The horizontal dashed line indicates the formation enthalpy at 0 K. The solid lines provide temperature-dependent values excluding and including anharmonic contributions. The dotted lines show the linear (Arrhenius) extrapolation of high-temperature formation enthalpies down to $T=0\mathrm{K}$ as discussed in the text.

Figure 5

Schematic illustration of formation energy vs Fermi level for an acceptor-type defect that can occur in two charge states: 0 and $-1$. Solid lines correspond to formation energies for the fully relaxed atomic geometry for each charge state. Dashed lines correspond to formation energies of configurations where the defect geometry of the other charge state has been frozen in. Both thermodynamic [$ϵ(0/-)$] and optical transition levels ($E_{\text{opt}}^{-/0}$ and $E_{\text{opt}}^{0/-}$) are shown.

Figure 6

Schematic illustration of pressure dependence of band edges and defect levels.

Figure 7

Formation energy vs Fermi level for a defect that can occur in $q=+1$, 0, or $-1$ charge states, illustrating the definition of the $U$ parameter. Solid lines correspond to formation energies for the fully relaxed atomic geometry for each charge state. Dashed lines correspond to formation energies where the defect geometry of the neutral charge state has been frozen in. Charge-state transition levels derived from the fully relaxed geometries (i.e., the thermodynamic transition levels) are marked by $ϵ(q/q^{\prime })$, and the levels where the geometry has been frozen in by $E_{\text{opt}}^{q/q^{\prime }}$. The latter also correspond to optical transition levels (see Sec. 2e). The difference between two subsequent transition levels is called the $U$ parameter. Note that the purely electronic $U$ parameter ($U^{\text{el}}$) is positive whereas the effective $U^{\text{eff}}$ parameter in this example is negative.

Figure 8

Configuration coordinate diagram for a ${\mathrm{N}}_{\mathrm{O}}$ substitutional impurity in ZnO, illustrating the difference between thermodynamic and optical transition levels. The negative charge state is considered the ground state, and the variation of the energy as a function of atomic displacements around the stable configuration is shown. The curve for ${\mathrm{N}}_{\mathrm{O}}^0$ is vertically displaced from that for ${\mathrm{N}}_{\mathrm{O}}^{-}$ assuming the presence of an electron at the CBM. $E_{\text{rel}}$ is the relaxation energy that can be gained, in the negative charge state, by relaxing from configuration $\{R{\}}_{q=0}$ (the equilibrium configuration for the neutral charge state) to configuration $\{R{\}}_{q=-1}$ (the equilibrium configuration for the negative charge state). This relaxation energy is sometimes called the Franck-Condon shift. The peak energies that would be observed in optical absorption or emission experiments are indicated. Adapted from [213].

Figure 9

Schematic illustration of defect-band dispersion for a $p$-like defect state ($s$-like states would have a minimum at $\mathrm{\Gamma }$). The shaded area indicates the effect of occupation according to a Fermi-Dirac distribution: electrons accumulate in the lower-lying, bonding parts of the defect band, shifting the occupied average away from the desired value for the isolated defect.

Figure 10

Calculated formation energy of the unrelaxed, neutral vacancy in diamond as a function of $\mathbf{k}$-point sampling [Monkhorst-Pack mesh ([238]) with offset $[\frac{1}{2},\frac{1}{2},\frac{1}{2}]$] for different supercells. Solid lines: constant defect-band occupations. Dashed lines: Fermi distribution ($k_{B}T=0.05\mathrm{eV}$).

Figure 11

Schematic illustration of origin and supercell dispersion of a shallow donor state. The defect gives rise to a localized state (resonance) within the conduction band. The electrically active level is not directly associated with this localized state, but rather arises from a perturbed host state below the CBM, offset in energy by an approximately constant amount.

Figure 12

Alignment problem for charged defects: Isolated defects (dashed lines) have a well-defined asymptotic limit for the $1/r$ Coulomb potential (thin solid line), where the potential can be aligned to the bulk. From the periodic array (thick solid line)—here aligned at the defect centers—the bulk limit cannot be determined.

Figure 13

Energy vs number of electrons for the LDA, Hartree-Fock theory, and a hybrid functional (HSE) for a ${\mathrm{Si}}_4{\mathrm{H}}_4$ cluster. The LDA underestimates the discontinuities at integer numbers, yielding a convex behavior below the straight line, whereas the Hartree-Fock calculation yields a concave behavior lying above the ideal straight line (dotted lines). The KS (HF) eigenvalues of the LUMO and HOMO correspond to the derivatives for $20\pm \delta$ electrons.

Figure 14

Calculated electronic band structures of ZnO using the LDA (left) and the $\mathrm{LDA}+U$ (right). The band alignment between LDA and $\mathrm{LDA}+U$ was taken into account (see text and Fig. 15). The lowest-energy bands between $-6.5$ and $-4.5\mathrm{eV}$ in the LDA ($-8$ and $-7\mathrm{eV}$ in the $\mathrm{LDA}+U$) are derived from Zn $3d$ states; the bands between $-4.5$ and 0 eV in the LDA ($-6.0$ and $-0.7\mathrm{eV}$ in the $\mathrm{LDA}+U$) are derived mostly from O $2p$ states; the bands above 0.8 eV in the LDA (1.5 eV in the $\mathrm{LDA}+U$) are conduction-band states, with the lowest conduction-band states derived mostly from Zn $4s$ states. The character of the bands was determined by projecting a given band on atomic-orbital states centered on the Zn and O atoms.

Figure 15

Calculated band alignment at a hypothetical ${\mathrm{ZnO}}^{\mathrm{LDA}}/{\mathrm{ZnO}}^{\mathrm{LDA}+U}$ interface, showing the effects on individual band edges of including the on-site Coulomb interaction $U$ (with a value $U=-4.7\mathrm{eV}$) for Zn $3d$ states. All values are in electron volts.

Figure 16

Formation energy as a function of Fermi level for an oxygen vacancy ($V_{\mathrm{O}}$) in ZnO. Energies according to the (a) LDA and (b) $\mathrm{LDA}+U$ calculations for Zn-rich conditions. The 0, $1+$, and $2+$ charge states are shown and the calculated band gaps are indicated. (c) Energies according to the LDA and $\mathrm{LDA}+U$ extrapolation scheme described in Sec. 4d4 ([148]). (d) Energies obtained with the HSE hybrid functional described in Sec. 4g1 ([259]). (c), (d) The formation energies for Zn-rich (lower curve) and O-rich conditions (upper curve) are shown for the 0 and $2+$ charge states. The position of the $ϵ(2+/0)$ transition level is indicated in (c) and (d).

Figure 17

Feynman diagrams for the screened interaction $W$ in the random phase approximation. In the random phase approximation, an “incoming” field $v$ can create a particle-hole pair which annihilates ($P$) creating a new field $v$. As a response to this induced field, another particle-hole pair can be created corresponding to the third term, and this process continues ad infinitum. The closed particle-hole diagrams are usually termed “bubble” diagrams; they are summed up to infinity.

Figure 18

Selected Feynman diagrams for the vertex connecting a Coulomb line with two propagators (vertex $\mathrm{\Gamma }$). The diagrams shown are the particle-hole ladder diagrams and correspond to the diagrams typically used in the Bethe-Salpeter equation for the calculation of optical properties ([7]; [302]). This is the simplest possible vertex correction and it describes the electrostatic interactions between particles and holes.

Figure 19

Results for $G_0{W}_0$ and $G{W}_0$ and self-consistent quasiparticle $GW$ ($\mathrm{scQP}GW$) band gaps, along with results using the semilocal Perdew-Burke-Ernzerhof (PBE) functional ([318]). Lattice constants are from low-temperature experiments, where available. Note the double logarithmic scale.

Figure 20

Diagonal part of the electronic contribution to the inverse of the macroscopic dielectric function vs wave vector $g$. DFT RPA screening results for magnesia (MgO) (circles) and GaAs (squares) are shown. The lines are fits to the calculated data. Thomas-Fermi screening for GaAs (broken line) and the hybrid functional HSE (full line) are shown as well.

Figure 21

Band gaps for the HSE03 ([131]) and for a modified HSE functional with $\mu =0.5{\AA }^{-1}$ and $\alpha =0.6$. The latter yields consistently improved band gaps. Also shown are the band gaps for $G_0{W}_0$ calculations using HSE03 orbitals and one-electron energies. Note the double logarithmic scale.

Figure 22

Formation energy as a function of Fermi level $E_F$ for ${\mathrm{N}}_{\mathrm{O}}$ in ZnO, for Zn-rich conditions. $E_F$ is referenced to the VBM, and the position of the transition level $(0/-)$ is indicated. Adapted from [213].

Figure 23

(a) Thermal expansion coefficient and (b) isobaric heat capacity of aluminum including electronic, quasiharmonic, anharmonic, and vacancy contributions. Experimental values (divided into pre-1950 and post-1950 data) are included for comparison. The melting temperature $T^m$ (933 K) is indicated by the vertical dashed line. At $T^m$, the crosses indicate the sum of all numerical errors (e.g., pseudopotential error, statistical inaccuracy, etc.) in all contributions for the GGA. Adapted from [114].

Figure 24

Calculated contributions from vacancies (vac), electronic (el), and anharmonic (ah) excitation mechanisms to the Gibbs energy $G$, the expansion coefficient $\alpha$, and the isobaric heat capacity $C_P$ of aluminum (all quantities at zero pressure). Thick (thin) lines show GGA (LDA) results. The black solid lines show the sum of all excitations. The right axes are scaled with respect to the “full” GGA values at the melting temperature $T^m$ (indicated by the vertical dashed line). The inset enlarges the vacancy contribution to $G$ at high temperatures. Adapted from [115].