4.11 Formation energies of charged defects

The Gibbs free energy of formation of a defect is given by

$\mathrm{\Delta }{\text{G}}_{\text{f}}^{\text{D}}={\text{E}}_{\text{tot}}^{\text{D}}-{\text{E}}_{\text{tot}}^{\text{perf}}-{\displaystyle \sum _{\text{i}}}{\text{n}}_{\text{i}}{\mu }_{\text{i}}^{\text{ref}}+\text{q}{\epsilon }_{\text{F}}^{\text{ref}}-{\displaystyle \sum _{\text{i}}}{\text{n}}_{\text{i}}\mathrm{\Delta }{\mu }_{\text{i}}+\text{q}\mathrm{\Delta }{\epsilon }_{\text{F}},$

(4.15)

where ${\text{E}}_{\text{tot}}^{\text{D}}$ and ${\text{E}}_{\text{tot}}^{\text{perf}}$ are the total energies of the defected and the perfect sytem, n${}_{\text{i}}$ is the number of atoms of type i added ($>0$) or removed (), $\mathrm{\Delta }{\mu }_{\text{i}}$ are the corresponding atomic chemical potentials referenced to ${\mu }_{\text{i}}^{\text{ref}}$, $\mathrm{\Delta }{\epsilon }_{\text{F}}$ is the Fermi level referenced to ${\epsilon }_{\text{F}}^{\text{ref}}$ and q is the charge of the system.
A common choice as a reference for the electron chemical potential ${\epsilon }_{\text{F}}^{\text{ref}}$ is the valence band maximum (VBM), so that the Fermi level can be assumed in the range between the VBM and the conduction band minimum (CBM), but in principle the choice of references for the chemical potentials is arbitrary.
FHI-aims uses the Ewald summation technique to calculate the electrostatic Hartree potential for a periodic system. For a charged periodic system (specified by the keyword charge in control.in) a neutralizing homogeneous background charge density is introduced to remove the divergent $𝐆=0$ component of the long-range part of the electrostatic potential.
This scheme is not suitable for periodic surface models, because the background charge density would be spread over the whole unit cell including the vacuum region. Instead charged surface defects can be treated within a virtual crystal approach (VCA), which corresponds to distributed doping of the material. The following scheme can be used for an insulating system with a localized defect level in the bandgap. By modifying the charge of the atomic nuclei (using the keyword nucleus in control.in), while keeping the system neutral, additional delocalized states can be introduced at the top of the valence band or at the bottom of the conduction band. The occupation of the defect levels can thus be tuned by the amount of charge distributed on the cations or anions in the system. To ensure that the defect has the desired charge q, the sum of the modified nuclear charges ${\text{Z}}_i^{\prime }$ should differ from the sum of the original nuclear charges ${\text{Z}}_i$ by the value of q:

${\displaystyle \sum _i^{{\text{N}}_{\text{atoms}}}}{\text{Z}}_i^{\prime }={\displaystyle \sum _i^{{\text{N}}_{\text{atoms}}}}{\text{Z}}_i-q.$

Note that for calculating the formation energy of a charged defect within the VCA the reference system should be the doped undefected system, not the perfect undoped system. Since doping pins the Fermi level $\mathrm{\Delta }{\epsilon }_{\text{F}}$ vanishes for this method.
For example, a way to model a positively charged oxygen vacancy at a metal oxide Me_xO_y surface is to distribute the charge uniformly on the metal atoms Me by changing their nuclear charge from Z(Me) to

$\text{Z}({\text{Me}}^{\text{VCA}})=\text{Z}(\text{Me})-{\displaystyle \frac{\text{q}}{\text{N(Me)}}},$

where N(Me) is the number of metal atoms in the system. This introduces vacant states at the VBM which in the defected system will be occupied by electrons from the defect level.
When using the neutralizing background method for bulk systems the additional term may introduce an arbitrary shift, so that it is necessary to find a common energy reference for the charged and the neutral periodic system to which the respective potentials can be aligned. For example alignment of the core levels of an atom far away from the defect can be done according to

$\mathrm{\Delta }{\epsilon }_{\text{F}}=({\epsilon }_{\text{F}}-{\epsilon }_{\text{core}}^{\text{D}})-({\epsilon }_{\text{VBM}}^{\text{perf}}-{\epsilon }_{\text{core}}^{\text{perf}}).$

Plot the atom projected density of states (output option output atom_proj_dos in control.in) for this atom for the charged defected and the neutral perfect system to visualize changes in the core states. (Be aware that in an all-electron approach the deeply lying core states are sensitive to local changes in electron density due to relaxation and charge redistribution, so that their shift in a defected system with respect to the perfect host system may not include only the average potential shift.)
Due to spurious electrostatic interaction as a result of the employed periodic boundary conditions the formation energy of a charged defect depends on the dimensions of the supercell. The formation energy scales as $\mathrm{\Delta }{\text{G}}_{\text{f}}^{\text{D}}(L)\approx \text{a}\frac{1}{L}+\text{c}$ for sufficiently large supercells [206]. For a simple cubic unit cell $L$ corresponds to the supercell lattice constant and can take up integer multiples of the unit cell lattice constant L⁽⁰⁾. For differently shaped unit cells with lattice constants L${}_{\text{1}}{}^{}{}_{}{}^{(0)}$, L${}_{\text{2}}{}^{}{}_{}{}^{(0)}$, L${}_{\text{3}}{}^{}{}_{}{}^{(0)}$ set for example L:=L${}_{\text{1}}$ and build supercells $L_{\text{1}}=\text{n}\cdot {\text{L}}_{\text{1}}^{(0)}$, $L_{\text{2}}=\text{n}\cdot {\text{L}}_{\text{2}}^{(0)}$, $L_{\text{3}}=\text{n}\cdot {\text{L}}_{\text{3}}^{(0)}$ with integer n. The desired formation energy of a single defect in an infinite supercell $\mathrm{\Delta }{\text{G}}_{\text{f}}^{\text{D}}(L\to \mathrm{\infty })$ can then be obtained by extrapolation. Note, that the convergence of the extrapolated energy with respect to the supercell size should be tested carefully. Taking into account geometric relaxation can improve the convergence significantly. Alternatively postprocessing correction schemes that allow to remove the spurious interaction terms have been suggested in literature [206, 96].