Defect formation energy
***********************

BAND Quantum ESPRESSO Semiconductors Materials Science Defect
formation energy Charged defects

Point defects are omnipresent in materials and influence their
electrical and optical properties. Understanding the formation of
defects is critical for many industries in the area of physics and
materials science. For example, this knowledge can be used to:

* design materials with specific and improved properties

* understand quality and ensure the reliability of materials

* gain insights into the manufacturing processes

Therefore, first principles modeling of point defects has become an
invaluable tool for understanding materials properties. Point defects
range from simple vacancies to intrinsic or extrinsic substitutions
and interstitials.

   [image]

Because defects are usually present in very low amounts in crystalline
materials (ppm), it is often of interest to examine the properties of
a single defect in a pristine crystal. However, because of periodic
boundary conditions, the size of the supercell must be sufficiently
large to neglect the interactions between the periodic replicas of the
defects. In the situation of charged defects, the strong and long-
range nature of the Coulomb interaction between the localized charge
distributions converges slowly. Consequently, the sizes of the
supercells needed to produce accurate energies become excessively
large. Thankfully, an ensemble of methods [2] [3] have been proposed
to correct for the spurious periodic interaction in small supercells
to enable accurate computation of both neutral and charged defect
formation energies, from first principles, and with reasonable
computing resources.

In this tutorial we will learn how to compute defect formation energy
in solids based on density functional theory. We will present the
procedure to compute two types of defects:

* neutral defects

* charged defects

Note:

  In the following, we will use the indices *p*, *0*, and *q* for
  perfect, neutral defect, and charged defect quantities. Moreover, to
  accelerate the calculations we will only consider unrelaxed defect
  cells, unless specified.

Note:

  To perform perform proper calculations on periodic systems, the
  number of kpoints should be wisely selected based on convergence
  tests, for each supercell size.


Neutral defects
===============

Considering neutral defects, the defect formation energy can be
expressed as:

   E^f_0 = E_0 - E_p - \sum_i n_i\mu_i

with E^f_0, E_0, and E_p the defect formation energy, the total energy
of the neutral defective structure, and the total energy of the
perfect structure, respectively. The coefficients n_i represent the
number of atoms of element i forming the defect (negative if removed,
and positive if added) with chemical potentials \mu_i, being reference
energies for the added/removed atoms. Therefore, computing the neutral
defect formation energy corresponds to finding the energies of the
neutral defective structure, the perfect structure, and the chemical
potentials of the atoms participating in the defect. Let's have a look
at some typical defects including vacancy, substitution, and
interstitial. For example we will compute the defect formation energy
corresponding to:

* C vacancy in diamond

* O substituted with N in MgO

* H interstitial in Si


Vacancy formation energy
------------------------

We will demonstrate how to compute the vacancy formation energy in
diamond with BAND.

   Start **AMSjobs**
   **1.** from the **menu bar** select **SCM → New input**
   **2.** switch to the BAND engine [image: ADFPanel][image] **→** [image: BANDPanel][image]
   **3.** from the **menu bar** select **Builders → Crystal Structure → Cubic → Diamond**
   **4.** in the builder panel, select **C** from the **Presets** dropdown menu and click **OK**
   **5.** from the **menu bar** select **Edit → Crystal → Convert To Conventional Cell**
   **6.** from the **menu bar** select **Edit → Crystal → Generate super cell..**
   **7.** in the **Super Cell Tool**, set the diagonal values to **3** to generate a 3x3x3 supercell then click **OK**
   **8.** click on [image: MoreBtn][image] next to **Numerical quality**, and in **K-space** and select **Good** from the dropdown menu

Note:

  There is a shortcut to **Edit → Crystal** (in the **menu bar**) by
  clicking the [image: CrystalTool][image] icon at the bottom of the
  molecule panel in AMSinput.

Note:

  It can also be a good idea to enforce symmetry of the crystal in
  order to optimize the performances of the DFT calculation.

Furthermore, we will mark the position of the future vacancy by
creating a region containing the C atoms that we will delete in the
next calculation. This will ensure the exact same number of symmetry
operators in all the calculations.

   **1.** select a C atom
   **2.** go to **Model → Regions**
   **3.** click the **+** symbol left of **Regions**
   **4.** rename the region as vacancy

   [image]

Later on we want to align the potential of a charged defect, and that
is why two extra steps are needed

   **1.** Set the origin to the C atom in the vacancy region (because we want to align the potential far away from the vacancy, at the edge of the unit cell): with the C atom selected, click on **Edit → Set Origin → X,Y and Z**,
   **2.** Using the looking glass icon ([image: Search][image]), search for the option **UpdateStdVec**, and disable it ("Programmer%UpdateStdVec") in the **Details → Expert BAND** panel. Otherwise the systems may be shifted with respect to each other, spoiling the alignment, see for instance Fig. %s.

You can keep all defaults parameters, save the job as
**C_diamond_3x3x3_p** and run the calculation. This corresponds to
running a single point DFT calculation with the LDA exchange-
correlation potential and the DZ basis set. The number of k-points was
set to Good, which will result in 3 kpoints along each lattice
direction.

Once the calculation is over, we can compute the energy of the
defective cell.

   **1.** from **AMSjobs** click the BAND icon next to the job you just executed to open **AMSinput**
   **2.** select the **C** atom that we previously tagged by the Regions by left-clicking on it and delete the atom with the **delete** key
   **3.** you can also delete the region from **Model → Regions** and then click the **-** symbol left the vacancy region

   [image]

Save the job as **C_diamond_3x3x3_0** and run the calculation.

To create the vacancy, we have simply removed a carbon atom,
therefore, -\sum_i n_i\mu_i = +\mu_C. In general, the chemical
potential of the atoms that participate in the defects are taken as
the energy of that atom in its ground state structure. In the case of
carbon, we can use the total energy of the diamond structure as
reference.

We can summarize the results in the table below:


Summary of the total energies computed with LDA for C.
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

+-----------------------------------+-----------------------------------+-----------------------------------+
| Calculations                      | Perfect diamond (p)               | Diamond vacancy (0)               |
|===================================|===================================|===================================|
| Number of atoms                   | 216                               | 215                               |
+-----------------------------------+-----------------------------------+-----------------------------------+
| Energy (eV)                       | -2118.32                          | -2099.94                          |
+-----------------------------------+-----------------------------------+-----------------------------------+

The neutral vacancy formation energy in C diamond can be calculated
as:

   E^f_0 = E_0 - E_p +\mu_C = -2099.94 - -2118.32 + -2118.32/216 =
   8.57 \text{ eV}

It is interesting to appreciate how increasing the supercell affects
the vacancy formation energy. You can repeat the previous calculations
in a smaller (2x2x2) and larger (4x4x4 and 5x5x5) supercells.

   [image]

We can conclude that from the supercell size of 4x4x4, the vacancy
interaction becomes negligible.


Some considerations
-------------------

Here follow some thoughts on the choice of the chemical potentials and
the k-grid.


Chemical potential
~~~~~~~~~~~~~~~~~~

The chemical potential represents the energy of a single atom being
taken from or added to some reservoir. As such there is some
arbitrariness, and it is important to specify how the reference
systems are chosen. Above we have chosen for the Carbon atom the
energy of the Diamond crystal per atom. Alternative choices could have
been the energy per atom in Graphite, or a Bucky ball. Sometimes
combinations are taken: for instance in the work of Komsa [2] the
chemical potential of Ga is defined as the energy of the GaAs crystal
minus the one from As

   \mu_\text{Ga} = \mu_\text{GaAs} - \mu_\text{As}

In an experimental setup the chemical potential may depend on the
temperature and pressure.


K-points
~~~~~~~~

As neutral defects formation energies converge quite quickly with
super cell size no super cell extrapolation is needed, and it may be
beneficial to use the "GammaOnly" k-grid (with only one k-point). For
small super cells the "GammaOnly" calculations are far from k-space
converged, but you can go to bigger super cells, that may be cheaper
to calculate than a smaller unit cell with more k-points. With
increasing super cell size, the "GammaOnly" grid becomes progressively
more accurate.

   [image]

In the figure we see the vacancy formation energy for several choices
of the "KSpace%Quality" (obtained with DFTB). While for the 4x4x4
super cell (abbreviated as super cell 4) the energy may be k-space
converged with "KSpace%Quality=Good", a better value is obtained with
"GammaOnly" at super cell=6. The latter may be faster.


Substitution
------------

In this second example, we will explore the case of a nitrogen
substitution in the oxide MgO. We start with the perfect supercell.

   **1.** in the **menu bar** of **AMSjobs** select **SCM → New input**
   **2.** switch to the QE engine [image: ADFPanel][image] **→** [image: QuantumESPRESSOPanel][image]
   **3.** click the magnifier icon and search for MgO
   **4.** **Edit → Crystal → Convert To Conventional Cell**
   **5.** **Edit → Crystal → Generate super cell..**
   **6.** click **OK** in the Panel to generate a 2x2x2 supercell
   **7.** set the **Wavefunction energy cutoff** to 60.0 Ry
   **8.** set 2 **K-points** in each dimension

   [image]

Save the job as **MgO_2x2x2_p** and run the calculation. With these
settings, we are going to perform a single point calculation with QE
with a 60 Ry kinetic energy cutoff, 2x2x2 k-grid, at the PBE level.
Once the calculation is finished, you can reopen the calculation to
prepare the defective cell.

   **1.** in **AMSjobs** click the QE icon next to **MgO_2x2x2_p** to open the previous calculation
   **2.** select an O atom and click **Atoms → Change Atom Type** and select N (ignore the warnings)
   **3.** Go to **Model → Magnetization** and set the **Magnetization** to **Collinear**
   **4.** Fix the total magnetization to 1

Note:

  Here we need to set the magnetization because we now have an odd
  number of electrons. It is also possible to perform a charged state
  calculation with an even number of electron (see below).

   [image]

Save the job as **MgO_2x2x2_OsubN_0** and run the calculation.

Now we need 2 more calculations to compute the chemical potential of O
and N. The ground state of oxygen is taken as the energy per atom of
the \text{O}_2 molecule and that of nitrogen is taken as the energy
per atom of the \text{N}_2 molecule. We can start from the
**MgO_2x2x2_OsubN_0** calculation in order to keep settings from the
previous job.

   **1.** in **AMSjobs** click the QE icon next to **MgO_2x2x2_OsubN_0**
   **2.** select all the atoms and delete them (hold shift and drag the mouse over the crystal and tap delete)
   **3.** go to **Model → Lattice**, set the diagonal of the **Lattice vectors** to 12 Å
   **4.** select the O icon and click twice in the molecule area followed by the **esc** key
   **5.** select the two O atoms and set the distance to 121 pm (you can either input the distance in the input cell at the bottom-right of the molecule area or adjust the distance with the slide bar at the bottom of the molecule area)
   **6.** set the **Task** to **Geometry Optimization**
   **7.** from the **Magnetization** tab fix the total magnetization to 2

   [image]

Save the file as **O2** (**File → Save As..**) and run the
calculation. Once finished, open the calculation in AMSinput to create
the \text{N}_2 run.

   **1.** in **AMSjobs** click the QE icon next to **O2**
   **2.** select all atoms and click **Atoms → Change Atom Type** to N
   **3.** select the two N atoms and set the distance to 110 pm
   **4.** remove the 2 in **Fix total magnetization** and set the **Magnetization** to **None**

   [image]

Save the job to **N2** and run the calculation. We can summarize the
results in the table below:


Summary of the total energies computed with LDA for MgO.
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

+----------------------+----------------------+----------------------+----------------------+----------------------+
| Calculations         | MgO (p)              | MgO with N           | O_2 molecule         | N_2 molecule         |
|                      |                      | substituted (0)      |                      |                      |
|======================|======================|======================|======================|======================|
| Number of atoms      | 64                   | 63                   | 2                    | 2                    |
+----------------------+----------------------+----------------------+----------------------+----------------------+
| Energy (Ry)          | -1208.594            | -1197.701            | -41.523              | -20.192              |
+----------------------+----------------------+----------------------+----------------------+----------------------+

To create the substitution, we have removed a O atom and added a N
atom, therefore, -\sum_i n_i\mu_i = +\mu_O - \mu_N. The substitution
formation energy in MgO can be calculated as:

   E^f_0 = E_0 - E_p +\mu_O - \mu_N = -1197.701 - -1208.594 +-41.523/2
   - -20.192/2 = 0.228 \text{ Ry}

We found that the formation energy corresponding to an oxygen
substituted with nitrogen in MgO is 6.20 eV.

Note:

  The supercell size for this example is probably too small to neglet
  the defect interaction between periodic images, and one should check
  convergence of the defect formation energy as a function of
  supercell size.

You can also verify that the numbers of kpoints used for the
calculations of the defect formation energy were sufficient by
repeating the calculation with various number of kpoints. This can be
achieved with this "PLAMS script". Remember you can always create a
PLAMS script from an AMSinput with **File → Export PLAMS script...**.
From there you you get all the PLAMS settings for the calculation and
you can make yourself a custom workflow!

   [image]

Based on the above convergence test, we can conclude that using a
2x2x2 k-grid to described the defect formation energy of N substituted
MgO leads an approximate error of ~0.05 eV as compared to a fine 5x5x5
k-grid. Therefore, a 3x3x3 k-grid is recommended.


Band and DFTB calculations
~~~~~~~~~~~~~~~~~~~~~~~~~~

Remember that we are interested in the super cell limit. If you use
the "GammaOnly" k-grid the super cell dependence of the substitution
energy is

   [image]

We can conclude that BAND predicts a super cell converged substitution
energy of about 6 eV with a DZ basis. Using a TZP basis set this
becomes about 6.35 eV.

The band procedure can be summarized as.

* Start from the minimal cell of MgO

* Make a larger super cell

* Make a region for the atom that you want to substitute

* Use the GammaOnly k-grid.

* Use the "Unrestricted" keyword

* Run the pristine and substituted system

* The N2 and O2 molecules can be calculated without an artificial
  lattice, i.e. as true molecules.

* There is no need to set the spin explicitly. The "Unrestricted" is
  only really needed for the substituted systems and for the O2
  molecule.

If you try this with DFTB a very different value is obtained of -30
eV. DFTB is not able to describe substitution energies.


Extrinsic interstitial
----------------------

For the last example, we will compute the formation energy of a
neutral H interstitial in Silicon. Neutral H interstitial will seek
out the regions of highest electronic charge density, which
corresponds in Silicon crystals at the midway between two Si atoms.
Let's start by computing the energy of the perfect supercell.

   Start **AMSjobs**
   **1.** from the **menu bar** select **SCM → New input**
   **2.** switch to the BAND engine [image: ADFPanel][image] **→** [image: BANDPanel][image]
   **3.** from the **menu bar** select  **Builders → Crystal Structure → Cubic → Diamond**
   **4.** in the builder panel, select **Si** from the **Presets** dropdown menu and click **OK**
   **5.** from the **menu bar** select **Edit → Crystal → Convert To Conventional Cell**
   **6.** from the **menu bar** select **Edit → Crystal → Generate super cell..**
   **7.** click **OK** in the panel to generate a 2x2x2 supercell

   [image]

Save the job to **Silicon_2x2x2_p** and run the calculation. We can
prepare the defective cell from the previous settings.

   **1.** from **AMSjobs** click the BAND icon next to the job you just executed to open **AMSinput**
   **2.** select two bonded **Si** atoms by left-clicking on them with the shift key pressed in
   **3.** from the **Atoms** menu select **Add Dummy At Center Of Selection**
   **4.** left-click the dummy atom and select **Atoms → Change Atom Type**, and choose H
   **5.** select the **Task** **Geometry Optimization**

Before running the geometry optimization, we need to make some room
for the H interstitial. This can be achieved as follow.

   **1.** select one of the bonded Si and move it slightly away from the H interstitial by right-clicking on it and moving the mouse in the desired direction
   **2.** select the second bonded Si and move it slightly away from the H interstitial

Ideally, you would like the distance between the H and neighbors Si to
be 1.6 Å.

   [image]

Save the job to **Silicon_2x2x2_0** and run the calculation. Now we
can calculate the chemical potential of hydrogen, taken as the total
energy per atom for \text{H}_2.

   **1.** from **AMSjobs** click the BAND icon next to the job you just executed to open **AMSinput**
   **2.** delete all the atoms in the supercell
   **3.** create two hydrogen atoms and set their distance to 0.77 Å
   **4.** set the **Periodicity** to **None**

   [image]

Save the job to **H2** and run the calculation. We can summarize the
results in the table below:


Summary of the total energies computed with LDA for Silicon.
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

+---------------------------+---------------------------+---------------------------+---------------------------+
| Calculations              | Silicon (p)               | Silicon with H            | H_2 molecule              |
|                           |                           | interstitial (0)          |                           |
|===========================|===========================|===========================|===========================|
| Number of atoms           | 64                        | 65                        | 2                         |
+---------------------------+---------------------------+---------------------------+---------------------------+
| Energy (eV)               | -338.7376                 | -340.3065                 | -6.6213                   |
+---------------------------+---------------------------+---------------------------+---------------------------+

To create the interstitial, we have added a H atom, therefore, -\sum_i
n_i\mu_i = - \mu_H. The interstitial formation energy in Silicon can
be calculated as:

   E^f_0 = E_0 - E_p - \mu_H = -340.3065 - -338.7376 - -6.6213/2 =
   1.74 \text{ eV}


Charged defects
===============

Considering now charged defects, the defect formation energy can be
expressed as:[2]

   E^f_q = E_q - E_p + E_\text{corr} - \sum_i n_i\mu_i + q[\epsilon_F
   + \Delta V_{0/p}]

with E^f_q and E_q, the charged defect formation energy and the total
energy of the charged defective structure, respectively. E_\text{corr}
is the finite-size supercell correction that we will ignore first and
discuss later. q is the charge of the defect, \epsilon_F the fermi
energy, and \Delta V_{0/p} is a term used to align the electrostatic
potential of the perfect and neutral defective supercells. Indeed, the
electrostatic potential between the perfect and neutral defective
supercells can be slightly different and is estimated by the
difference of the electrostatic potential away from the defect [2].
The term q\epsilon_F here is equivalent to the chemical potential of
the electron which is added if the charge is positive or subtracted if
the charge is negative, similar to the \sum_i n_i\mu_i term for the
atoms.

Note:

  Note that because of the periodic boundary conditions it is not
  possible to compute the energy of a charged supercell and most
  engines will introduce a compensating homogenous background charge
  to neutralize the total charge.

Note:

  Because of the computational limitations of density functional
  theory, it is not possible to perform a supercell calculation with
  similar concentration as in the real crystal. This limitation leads
  to a highly concentrated defect in the supercell. As we have seen
  above, in the case of neutral defects, the choice of a reasonably
  large supercell is enough to limit the periodic interactions between
  defects. However, because of the long-range effect of Coulomb
  interaction, the proper calculation of a charged defect must be
  extrapolated to an infinitely large supercell.

We provide below two examples based on diamond to compute the
formation energy of charged C vacancy using:

* the supercell method

* finite-size supercell correction


Comparing potentials of different geometries
--------------------------------------------

For the calculation of \Delta V_{0/p} we determine the difference in
the coulomb potential of the defect and the pristine system. By
default band may shift a system, to maximize the symmetry. If this
shift is different for the two systems we get a wrong result. To avoid
this we have to use "Programmer%UpdateStdVec=No".

   [image]The plane averaged potential difference requires that both
   systems have the same origin, otherwise a wrong result may be
   obtained, as shown with the blue curve. To avoid letting band
   change the origin use "Programmer%UpdateStdVec=No".


Supercell method
----------------

The supercell method consists of computing the formation energy of
charged defects for various supercell sizes, and to extrapolate the
defect formation energy of an isolated defect in the limit of infinite
supercell size. One of the major drawback of this method is that it
requires the calculations of at least 3 supercell sizes.

To illustrate the procedure, we will compute the -2 charged vacancy
formation energy in diamond (V^{-2}). We need to perform the
calculations corresponding to neutral perfect, neutral defect, and
charged defect supercells. We already performed calculations for
neutral perfect and neutral defect cells. We can just reopen the
neutral defect and set the charge to -2.

   **1.** from **AMSjobs** click the BAND icon next to the job we saved as **C_diamond_3x3x3_0**
   **2.** set the **Total charge** to -2

Save the job to **C_diamond_3x3x3_q** and run it. You can approve when
the pop-up will tell you that a neutralizing density will be added to
the charged density.

   [image]

Looking at the charged defect formation energy equation, the energies
can be extracted as usual, and the Fermi energy will be approximated
to the energy of the top of the valence band (~ -7 eV). We just need
to extract the potential alignment term V_{0/p} between the neutral
perfect and neutral defective cells. We will compute and save the
planar averaged of the Coulomb potential and execute a Python script
to extract the value of V_{0/p}.

   **1.** in **AMSjobs** select the neutral perfect job and click **SCM → View**
   **2.** select **Fields → Grid → Medium**
   **3.** select **Add → Graph Along Line**
   **4.** next the **Use full cell** select **Line along b** in the dropdown menu. (Line along a and along c are equally good choices, as the cell is a cube)
   **5.** select **Averaged** next to **Planes S are to be**
   **6.** check the button next to \int_S F(r) dS and uncheck F(r)
   **7.** click **Field → Potential → Coulomb Potential**

   [image]

   [image]

Wait for the planar averaged to be completed. In the **Coulomb
Potential Line Plot** window, click the **Save XY Values** and save
the plot as **C_diamond_3x3x3_p_coul.xy**.

   [image]

Reproduce the previous steps to extract the Coulomb potential of the
neutral defective cell and save the graph as
**C_diamond_3x3x3_0_coul.xy**. Plotting both potentials (using for
instance a spreadsheet program) we see that the pristine and defect
potentials are very similar with some difference near the center of
the cell, exactly as hoped for

   [image]

Download "this" Python script and execute it in the same directory of
the two XY files saved previously.

   [image]

This is the Coulomb potential corresponding to the neutral vacancy.
The value of the potential difference away from the vacancy gives us
V_{0/p} which is also outputted by the Python script. We can now
summarize the calculations.

   E^f_q = E_q - E_p - \sum_i n_i\mu_i + q[\epsilon_F + \Delta
   V_{0/p}]


Summary of the total energies computed with LDA for C (all in eV).
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

+----------------+----------------+----------------+----------------+----------------+----------------+----------------+
| Calculations   | C (p)          | C (0)          | C (q)          | E VB (p)       | ΔV_{q/0}       | E^f_q          |
|================|================|================|================|================|================|================|
| 2x2x2          | -627.68        | -609.58        | -620.57        | -7             | -0.097         | 11.49          |
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+
| 3x3x3          | -2118.32       | -2099.94       | -2110.44       | -7             | -0.037         | 12.15          |
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+
| 4x4x4          | -5021.28       | -5002.84       | -5013.03       | -7             | -0.005         | 12.45          |
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+
| 5x5x5          | -9807.05       | -9788.61       | -9798.59       | -7             | -0.011         | 12.68          |
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+

The supercell method consists in reproducing the previous calculations
within various supercell sizes. We performed  calculations in 2x2x2
and 4x4x4 supercells and we provide below the plot of the charged
defect formation energy as a function of lattice parameter. By
extrapolation, we find the isolated defect formation energy to be of
13.5 eV.

   [image]

Note:

  Comparing the previous charged defect formation energy and the
  neutral defect formation energy scalings we can appreciate the
  strong long-range effect of the charge Coulomb potential compared to
  that of the neutral defect. Indeed, while a 4x4x4 supercell was
  large enough to compute an accurate neutral defect formation energy
  in diamond, here we see that the same supercell is still off by
  almost 1 eV for charged defect!


Finite-size supercell correction
--------------------------------

Several flavors of finite-size supercell corrections have been
proposed and we will describe the scheme of Freysoldt, Neugebauer, and
Van de Walle (FNV) [1]. The idea of the finite-size supercell
correction is that since direct first principle calculations of the
charged defect formation energy with various supercell sizes is
computer intensive, one would only perform calculations in a
relatively large supercell and an approximated model will be used to
correct for the finite-size supercell interaction including the proper
isolated defect interactions. In the FNV scheme, the defect density is
approximated with a Gaussian function and the corresponding potential
and energy are calculated by solving the Poisson equation. The
correction term is defined as

   E_\text{corr} = E_{lat} - q\Delta V_{q/0}

where E_{lat} corresponds to the lattice finite-size correction and
\Delta V_{q/0} the potential alignment between the model and the DFT
calculation. The model supercell calculation for E_{lat} involves
solving the Poisson equation for the periodic model potential,
V_{lat}(r).

   \nabla\cdot [\epsilon (r)\nabla V_{lat}(r)] = -4\pi\rho(r)

with \epsilon (r) the dielectric tensor profile of the materials and
\rho(r) the model charge distribution. The dielectric tensor can be
computed from density functional theory, as described in **this
tutorial**. In a 3D bulk system with an isotropic dielectric,
V_{lat}(r) can be obtained as:

   V_{lat}(G) = \frac{4\pi \rho(G)}{\epsilon |G|^2}

The number of G vectors used in the calculation is determined by an
energy cutoff. To best define the charge density of the model Gaussian
distribution, one needs the width of the distribution. This can be
obtained by fitting the DFT defect charge density. First, you need to
identify the molecular orbital of the defect. This can be achieved by
looking at the DOS (or even better at the band structure).

   **1.** in **AMSjobs** select the charged defect job and click **SCM → DOS**

   [image]

We can clearly see a state in the bandgap corresponding to the defect
around -0.13 Hartree. Let's have a look at the orbitals.

   **1.** in **AMSjobs** select the charged defect job and click **SCM → View**
   **2.** from the **menu bar** select **Add → Graph Along Line**
   **3.** next the **Use full cell** select **Line along b** in the dropdown menu
   **4.** select **Averaged** next to **Planes S are to be**
   **5.** check the button next to \int_S F(r) dS and uncheck F(r)
   **6.** click **Field → Orbitals (occupied)**

   [image]

   [image]

We can see that the highest occupied orbital is indeed around -0.13
Hartree and is likely corresponding to the defect. Click on it and
AMSview will calculate the planar averaged.

   [image]

Save the graph as XY Values and name it **C_diamond_3x3x3_q_HOMO.xy**.
Download "this" Python script and execute it in the same directory of
the XY file you just saved.

   [image]

Now we know the shape of the defect, we can solve Poisson to get the
potential and the energy correction. Download "this" Python script
that will solve Poisson. You can edit the file for your specific case.
Parameters are the lattice, an energy cutoff for the basis set used in
the model, the dielectric profile (here simply a scalar), the width of
the Gaussian, the charge of the defect, and its location. The script
will solve Poisson in 2x2x2 up to 5x5x5 supercells by defaults.

   [image]

The script will output at the end the lattice correction defined above
E_{lat} = E_{iso} - E_{per}. Indeed, the lattice correction removes
for the spurious periodic interaction in the supercell E_{per} as well
as adds the isolated defect energy E_{iso} (extrapolated).

>>> 2x2x2 E_latt = E_iso - E_per =  2.821 -  0.985 =  1.836 eV
>>> 3x3x3 E_latt = E_iso - E_per =  2.821 -  1.530 =  1.291 eV
>>> 4x4x4 E_latt = E_iso - E_per =  2.821 -  1.835 =  0.986 eV
>>> 5x5x5 E_latt = E_iso - E_per =  2.821 -  2.026 =  0.795 eV

The script will also generate the potentials corresponding to each
model defective cells as **V_X_r.txt** files. We can therefore compare
the potential of the model defect to that of the DFT calculation
(taken as the difference between the charged and neutral defect). This
"this" can be used to get the validation plot below.

   [image]

The potential difference between the model and DFT gives us access to
the last term in the charged defect formation energy \Delta V_{q/0}.
To summarize, we have the following.

   E^f_q = E_q - E_p + E_\text{corr} - \sum_i n_i\mu_i + q[\epsilon_F
   + \Delta V_{0/p}]


Summary of the total energies computed with LDA for C (all in eV).
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

+--------------+--------------+--------------+--------------+--------------+--------------+--------------+--------------+
| Cells        | C (p)        | C (0)        | C (q)        | ΔV_{0/p}     | E_{latt}     | ΔV_{q/0}     | E^f_q corr.  |
|==============|==============|==============|==============|==============|==============|==============|==============|
| 2x2x2 (eV)   | -627.68      | -609.58      | -620.57      | -0.097       | 1.84         | 0.0          | 13.33        |
+--------------+--------------+--------------+--------------+--------------+--------------+--------------+--------------+
| 3x3x3 (eV)   | -2118.32     | -2099.94     | -2110.44     | -0.037       | 1.29         | 0.0          | 13.45        |
+--------------+--------------+--------------+--------------+--------------+--------------+--------------+--------------+
| 4x4x4 (eV)   | -5021.28     | -5002.84     | -5013.03     | -0.005       | 0.99         | 0.0          | 13.44        |
+--------------+--------------+--------------+--------------+--------------+--------------+--------------+--------------+
| 5x5x5 (eV)   | -9807.05     | -9788.61     | -9798.59     | -0.011       | 0.80         | 0.0          | 13.47        |
+--------------+--------------+--------------+--------------+--------------+--------------+--------------+--------------+

   [image]

With the smallest supercell 2x2x2 we only do an error of 0.1 eV while
the 3x3x3 gives an error of 0.01 eV with respect to the large 5x5x5
supercell.

The formation energies of charged defects depend on the Fermi level,
which is the reservoir with which the impurity exchanges electrons.
The Fermi-level position (traditionally taken with respect to the
valence band maximum) for which +1 and neutral states have equal
energies defines the +/0 transition level, that is, the donor level
\epsilon_D; similarly, the acceptor level \epsilon_A (i.e., the 0/-
transition) is defined as the point where neutral and -1 states have
equal energy [3] .

   [image]

Finally, 1D and 2D charged defect systems can also be corrected, as
described in these papers [4] [5] .


References
==========

[1] Freysoldt C, Neugebauer J, Van de Walle CG. Fully ab initio
    finite-size corrections for charged-defect supercell calculations.
    Physical review letters 102.1 (2009): 016402

[2] Komsa HP, Rantala TT, Pasquarello A. Finite-size supercell
    correction schemes for charged defect calculations. Physical
    Review B 86.4 (2012): 045112

[3] Herring C, Johnson NM, Van de Walle CG. Energy levels of isolated
    interstitial hydrogen in silicon. Physical Review B 64.12 (2001):
    125209

[4] Naik MH, Jain M. CoFFEE: corrections for formation energy and
    eigenvalues for charged defect simulations. Computer Physics
    Communications 226 (2018): 114-126

[5] Github CoFFEE package


Challenges in computing charged defects
=======================================
