Defect Formation Energy (DFE) of a point defect
The case you described is an example of a point defect (vacancy, substitution/impurity, interstitial) in bulk material. Such point defects can be neutral as well as electrically charged. In the context of periodic Density Functional Theory (DFT) based calculation, for a point defect with charged state $q$, the defect formation energy, $E_f^D(q,E_F)$, is defined as a function of the Fermi energy $E_F:$
$$ E_f^D(q,E_F) = \left[E^D(q) + E^D_{corr}(q)\right]-E^{Pristine}-\sum n_i\mu_i + q(E_F + E_{VBM}^{Pristine}-\Delta V_{0/b}).$$
Here,
- $E^D(q)$ is the total energy of the defective supercell,
- $E^D_{corr}(q)$ is the correction term due to the finite size effects (this also accounts for the compensating homogeneous electron gas in the background that is added by the code to prevent the divergence in case of charged defect calculation)
- $E^{Pristine}$ is the total energy of the pristine supercell,
- $n_i$ is the number of the $i$-th atoms being added to the pristine system to create the defective supercell (so $n=-1$ for a single vacancy),
- $\mu_i$ is the chemical potential of the $i$-th atoms calculated using the same pseudopotential,
- $E_{VBM}^{Pristine}$ is the valence band maximum of the pristine supercell. You can set the
verbosity='high' in the &CONTROL namelist and then you can check the occupations in the output and note the VBM.
- $\Delta V_{0/b}$ is the potential alignment term between pristine bulk supercell and neutral $(q=0)$ defective supercell.
For neutral defects, this gets simplified as:
Calculating DFE using Quantum ESPRESSO (QE)
In all of the following steps, I am assuming you did proper convergence tests before proceeding to the actual calculation. Now, as the first step, you have to build a supercell of the pristine material and run an SCF calculation (calculation='scf' using pw.x). This energy is the $E^{Pristine}$. Next, you have to introduce the defect in the supercell and run another SCF calculation. This will give you $E^D(q=0)$. You can then change the $q$ by using the tot_charge keyword in the &SYSTEM namelist and similarly obtain the energy for the corresponding charged states.
The $E^D_{corr}(q)$ term can be calculated in several ways. For example, here (arxiv) they did it in 3 ways:
- Manually applying 1st order Makov-Payne correction: $q^2\alpha/2\epsilon L$ (in atomic units) or $14.39952q^2\alpha/2\epsilon L$ (in $eV$ and $L$ in Angstrom)
- Using the
assume_isolated='makov-payne' in QE calculation which gives you higher order corrections
- Using tools such as SXDEFECTALIGN - see section 4.5 of this document about how to use this software with QE. Note that the potential alignment term can also be obtained using this.
Now, the only term left is the $-\sum n_i\mu_i$, the chemical potential terms. You can calculate them from their elemental configuration or from their stable configuration. The $\mu_i$ is basically the total energy per atom of the $i$-th species. This way you can calculate the defect formation energy as a function of the Fermi level.
Usually, people vary the Fermi level in the range of the band gap. For example, if your calculated band gap is 2 eV, then you will vary the Fermi level from 0 to 2 and plot the defect formation energy. Since the formula is of form: $ y = mx + c$ where $m$ corresponds to $q$, you will get straight lines plot like the following: (this is one of my calculation for substitutional Cu impurity in ZnS)
