To calculate the formation energy of a material AB with 2 atoms, I created a 2x2x2 supercell, resulting in a total of 16 atoms. Then, I removed one atom of type A from the supercell. From the calculations , I obtained the following energies (extracted from the output):

-Total energy of the defect-containing supercell: Edefect =-5193.89497567 Ry (for 15 atoms)
-Total energy of the perfect crystal unit cell: Eperfect =-697.76209902 Ry (for 2 atoms)
-Total energy of the isolated A atom: Evacancy =-462.13675745 Ry (for 1 atom)

To determine the formation energy, I initially used the formula:

Eform=Edefect (supercell)−Eperfect (unit cell)+Evacancy

Eform= (-5193.89497567/15) 
- (-697.76209902/2) 
+ (-462.13675745/1) 
= -459.515372984667 Ry

I have doubts about the correctness of the formula I used. Could someone please review and correct my approach?

PS: For my simulations I use Quantum ESPRESSO.

  • 1
    $\begingroup$ Since someone asked you this in Discord, it would help to have the question visible here too: are you still in active need of an answer to this question? $\endgroup$ Commented Mar 6 at 11:45

1 Answer 1


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_{q/b}).$$


  • $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)

enter image description here


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .