# 7. Diamond vacancy defect tutorial

## 7.1. Introduction

In this tutorial, I explain how to calculate the defect level diagram of a carbon vacancy in diamond. We will calculate the formation energy as a function of chemical potential of several charge states of the defect. The formation energy of a charged defect is given by

$$E_f[D^q] = E[D^q] - E[Host] - \sum n_s \mu_s + q E^{Fermi} + E_{corr}(q,\epsilon)$$

where $$E[D^q]$$ is the total energy of the defect system with charge q, $$E[Host]$$ is the total energy of the perfect system, $$n_s$$ is the number of atoms of species $$s$$ added (if positive) or removed (if negative). For instance, in the diamond vacancy, we have $$n_C = -1$$. $$\mu_s$$ is the chemical potential of species $$s$$, $$E^{Fermi}$$ is the chemical potential of the electrons and $$E_{corr}(q,\epsilon)$$ is a correction that account for finite size effects. I will use the correction scheme of Kumagai and Oba, which requires determining the dielectric constant of diamond (and more generally the dielectric tensor).

After obtaining the defect formation energies of several charge states, we will compute the following transition level diagram

I shall go through the following steps:

In this tutorial, I will use 64-atom supercells and neglect structure relaxation to minimize the computational intensity.

## 7.2. Vacuum atom

Since we will investigate a vacancy defect with atomic orbitals, we will introduce a so-called vacuum atom in place of the removed carbon atom. Simply put, we remove the pseudopotential but keep the orbitals. I create a vacuum atom from the pseudopotential file C_SZDP.mat as follows

load('C_SZDP.mat')
data = makeVacuum(data);
save Va_SZDP data


In what follows, a carbon vacancy will draw its atomic orbital basis from the file Va_SZDP.mat, and I shall use Va to label the element.

## 7.3. Neutral defect band structure

In this section, I compute the band structure of the neutral vacancy in diamond. This is to give myself an idea of the electronic structure of the defect, which will determine how we do things next. As we shall see, the diamond vacancy has a six-fold degenerate defect state (including spin degeneracy). But the finite supercell size yields an artificial dispersion of the defect levels across the Brillouin zone. This so-called band dispersion error is mitigated using a uniform occupation scheme instead of, say, Fermi-Dirac occupancies.

So we need to know the band structure to know how to populate it. We thus perform the ground state calculation scf.input and get a band structure from bs.input. The band structure looks like

The Fermi level clearly cuts through three defect bands around the center of the band gap. In the limit of an infinite system size, the three bands become flat and degenerate (as they are at the $$\Gamma$$-point). The dispersion is a finite size effect. Further inspecting the eigenvalues and occupancies reveals that the neutral defect level is 1/3 occupied. To mitigate the (incorrect) dispersion effects, we will adopt is fixed population scheme to force the occupancy to 1/3 across the entire Brillouin zone, as it would in the infinite size system. This not only improve convergence with respect to supercell size, but also with respect to k-sampling, as shown in this figure

Here, the subscript vac is for the regular occupancy scheme and fix for the fixed occupancy scheme.

## 7.4. Neutral formation energy

Now that I know the electronic structure, I’m ready to perform the supercell calculations. First, we’ll need the pristine system potential for reference. I thus set up a bulk supercell as bulk/scf.input. Note the following parameter is set as follows

option.savePotential = 1


to make sure to save the potential.

Next, let’s see what’s in the neutral directory. The script scf.input defines the parameters for the ground state calculation of the neutral defect. The important keywords are at the bottom of the file

nv = 64 * 4 - 4
atom.nvalence = nv - 0
kpoint.sampling = 'fixed'
kpoint.nocc = [ones(floor(nv/2)-1,1);2/6;2/6;2/6]
element(2).species = 'Va'
element(2).path = '../Va_SZDP.mat'


First nv is not a keyword, but simply a useful variable to write the rest of the file. It corresponds to the total valence of the neutral defect system, i.e. 64 atoms times 4 electrons, minus 1 atom times 4 electrons. The keyword atom.nvalence is set to nv - 0. We could simply set it to nv but this is to make it obvious that the system is neutral; later we’ll change 0 for something else when we study charged defects. The fixed occupancy is defined using the keywords kpoint.sampling and kpoint.nocc. The value of kpoint.sampling is self-explanatory. For kpoint.nocc, one needs to specify the occupancy of each band. In a vacancy defect system, there are floor(nv/2)-1 fully occupied bands, three partially occupied bands and then unoccupied bands. We only need to specify the occupancy of the (partially) occupied bands. I thus use the MATLAB function ones to generate a floor(nv/2)-1-entry vector of ones (the spin factor should not be included here). The following three defect levels should then be populated with 1/3 electron each, for a total of nv/2 electrons (again neglecting the spin factor). Finally, we add the vacuum atom definition with label Va. Note that this is all done in a programmatic way, but you could input all numbers in a more manual way if that makes you feel more comfortable.

Once the scf.input calculation is done, we can proceed to compute the formation energy of the neutral defect. This is done by the baseFormEnergy function. We must create a structure to store the following parameters:

This is done as follows

param.n = [-1];
param.mu = [-163.887300];
param.efermi = 10.126067;
ef = bareFormEnergy('bulk/results/scf','neutral/results/scf',param);
fprintf('ef = %f eV\n',ef);


which prints

ef = 7.808313 eV


Note the data in param are in units of eV. Usually, one includes ionic relaxation in the formation energy contribution by using the total energy of the relaxed defect supercell. If you have the time and resources, you may do the computation rlx.input instead of scf.input.

## 7.5. Charged defect formation energy

We now proceed to a charged defect formation energy. Charged defect supercell calculations are quite similar to the neutral defect one. The main difference is that we will use a correction formula to mitigate finite-size electrostatic errors. Let’s take the +2 charged defect for example (i.e. change directory to plus2). There is an input file scf.input quite similar to that of the neutral defect. It ends as follows

nv = 64 * 4 - 4
atom.nvalence = nv - 2
kpoint.sampling = 'fixed'
kpoint.nocc = [ones(floor(nv/2)-1,1);0/6;0/6;0/6]
element(2).species = 'Va'
element(2).path = '../Va_SZDP.mat'


As before nv is not a keyword but a variable which corresponds to the total valence of the neutral defect system. The keyword atom.nvalence is set to nv - 2; we remove two electrons making the total charge +2. The occupancy is defined using the keywords kpoint.sampling and kpoint.nocc. Remember that the three defect levels of the neutral defect were populated with 1/3 electron each. Here we remove two electrons (one in each spin channel) from those levels leaving zero electron. We thus set the occupancy of the last three bands to 0/6 (to keep the template file intact). Finally, we add the vacuum atom definition with label Va.

As before, the (bare) formation energy may be obtained with bareFormEnergy

param.n = [-1];
param.mu = [-163.887300];
param.efermi = 10.126067;
ef = bareFormEnergy('bulk/results/scf','neutral/results/scf',param);
fprintf('ef = %f eV\n',ef);


which yields

ef = 6.875340 eV


Now, this formation energy includes defect-defect electrostatic interactions which shouldn’t be present in the infinite-size, isolated defect system. Here, I’ll use the the correction scheme of Kumagai and Oba, which requires determining the dielectric constant of diamond. This may be done following the Dielectric constant tutorial. In the end, one gets a value close to 5.41. The electrostatic correction is performed by the function corrKO as follows

param.dielectric = 5.41*[1,0,0;0,1,0;0,0,1];
param.site = 1;
ec = corrKO('bulk/results/scf','plus2/results/scf',param);
fprintf('ec = %f eV\n',ec);


which yields

ec = 0.765893 eV


It is quite significant indeed, more than 10%. The interface of corrKO is similar to that of bareFormEnergy, but it requires different parameters. First, the dielectric tensor is required in param.dielectric. Here it is set to 5.41 times the identity matrix. Next, param.site is the atom index corresponding to the center of the defect (here the vacancy). The total formation energy is finally ef + ec = 7.641232 eV.

Again, one should usually include ionic relaxation in the formation energy contribution by using the total energy of the relaxed defect supercell.

## 7.6. Transition level diagram

In this section, we derive the transition level diagram from the defect formation energy. Out starting point is the corrected formation energy for all charged states. If you repeat the steps from the previous section for all charged states, you should find values close to

ef = [5.746588,6.592742,7.808313,9.448148,11.544052,14.094930,17.067092];


for the bare formation energy and

ec = [0.681561,-0.010808,0.000000,0.760446,2.300168,4.624176,7.712705];


for the corrections, which gives finally

ef = [6.428149,6.581935,7.808313,10.208593,13.844220,18.719106,24.779797];


You can read everything using the script readall.m. Recall the formula for the formation energy of a defect is

$$E_f[D^q] = E[D^q] - E[Host] - \sum n_s \mu_s + q E^{Fermi} + E_{corr}(q,\epsilon)$$

Therefore, as we vary $$E^{Fermi}$$ the most energetically favorable charge state will change. The crossing point between the formation energy of a charged state and another is a transition level. Those levels can be computed by the function transitionLevels.

ef = [6.428149,6.581935,7.808313,10.208593,13.844220,18.719106,24.779797];
q = [2,1,0,-1,-2,-3,-4];
trlvl = transitionLevels(q,ef);


Note that if the data is not in decreasing q-order, transitionLevels will reorder the data internally. The first entry in trlvl is thus the transition level between the highest charge state and the second highest, and so on.

I get

trlvl = [0.153786,1.226378,2.400280,3.635627,4.874886,6.060691];


in units of eV for the transition levels. This means that the first transition (+2/+1) occurs below the VBM and the last transition (-3/-4) above the CBM. Finally, this is vizualize calling the function plotTLD. One just needs to set the band gap (experimental or else) in param.gap as follows

param.gap = 5.47
plotTLD(q,trlvl,param)


which produces the following plot