Tutorial: Understanding and Computational Modelling of Defects in Semiconductors (with VASP) 

Thanks @Glaedr_4_life, really glad to hear it! 😃

Hi Jaegwan! Sorry for the slow reply, I was at APS March and on hols the last 2 weeks. To answer your questions: 1. Yes we need to use spin polarisation (ISPIN = 2) even if the material is non-magnetic. All odd-electron defects (e.g. V_Cd^-1 in CdTe) have an unpaired electron and so will be spin-polarised (and usually for this we would set NUPDOWN = 1 to help enforce this). For defects with even numbers of electrons in a non-magnetic host, in most cases spin polarisation is unnecessary, as all electrons will likely be paired, however for certain defects this may not actually be the case. For example, in our work on V_Cd in CdTe, we find a low-energy bi-polaron state for V_Cd^0, where you have two separate unpaired/spin-polarised hole polaron states, on two of the neighbouring Te atoms, which you don't find if you turn off spin polarisation with ISPIN = 1. Often the best practice is to allow spin polarisation for all defects, and then for even-electron defects after the initial relaxations, if their magnetisation (in the VASP OUTCAR) is zero, then can switch to ISPIN = 1 to be more efficient. 2. In general, setting MAGMOM is only recommended to be used if your host material is magnetic (in which case it should match the expected anti-/ferromagnetic etc ordering), or if you're trying to enforce localisation of an electron/hole polaron on a specific atomic site at/near the defect (in which case you would set it to match the expected electron & spin configuration for your defect supercell). Otherwise (i.e. most cases), we leave MAGMOM unset and use the VASP default

@ Thank you for your kind answer! Have a great day :)

Hi! Sorry I didn't see this before, it got caught in RU-vid's spam filter as it was duplicated ~10 times. 1) Yes this point can be confusing due to mixed terminology in the field. In both cases it is indeed the defect formation enthalpy, however often this is referred to as the defect formation 'energy' in the literature and orally. This is because most entropy contributions to the energy of defects (e.g. vibrational entropy) are usually considered negligible, such that ΔE = ΔH -TΔS ≃ ΔH. However, a key point here is that this neglects the configurational entropy contribution of defects, which is quite significant and in fact is the main reason defects form in all materials. I discuss this at 31:46 in the video. So in reality it's the formation enthalpy, which approximately matches the formation energy without the configurational entropy contribution, but people often still just refer to this as the 'defect formation energy' - confusing I know! 2) Nope this applies to all point defects. It is valid for extended and complex defects as well, however care needs to be made in these cases to adjust N (the number of potential sites for this defect) appropriately. N will be significantly reduced for extended or complex defects, due to the reduction in site degeneracy/multiplicity (i.e. number of potential sites to form). This is briefly discussed in pubs.rsc.org/en/Content/ArticleLanding/2022/FD/D2FD00043A, pubs.rsc.org/en/content/articlelanding/2017/ta/c6ta09155e, and chemrxiv.org/engage/chemrxiv/article-details/63d2a2c41fb2a8767ee1e06f. 3) E_F represents the Fermi level in the actual material, which depends on the formation energies of all defects in all charge states in that material (this is why we often call it the 'self-consistent Fermi level'), and we solve for it later when post-processing our defect calculations. It doesn't correspond to the Fermi level in any of the supercell calculations. This is discussed a little in: iopscience.iop.org/article/10.1088/2515-7655/aba081 4) There are different methods to calculate E_corr, the two most popular of which are the FNV (for isotropic systems) and eFNV (for anisotropic) methods. This is automatically calculated in our defect package doped: github.com/SMTG-UCL/doped (as it is in most defect computational packages). More discussion on this can be found in iopscience.iop.org/article/10.1088/2515-7655/aba081 and journals.aps.org/rmp/abstract/10.1103/RevModPhys.86.253 5) Yes! That term is independent of temperature. 6) Sure! We give that formula and discuss in our paper here: pubs.rsc.org/en/Content/ArticleLanding/2022/FD/D2FD00043A

Hi, I'd suggest you check out the tutorials on the doped code website: doped.readthedocs.io It shows the full process of generating and parsing defect and chemical potential calculations

Thank you so much

Hi @Nigel Hew, thanks for your kind words! 1) For the chemical potentials, it's done the same way as PyCDT or other defect packages do it, it's just that the chemical potential for elements in their standard states is formally zero by convention. So yes within the underlying calculation to get defect formation energies, we use mu_Zn = bulk energy/atom etc, but in chemistry/materials science papers it is always the 'formal chemical potential' (= energy/atom relative to the DFT energy/atom of the elemental phase) that we report. This is because the absolute energy/atom for a bulk calculation will depend on system settings such as DFT functional, but the relative energies (which give formation energies and chemical potential limits for example) should not depend so much on these. 2) The transition levels are 'thermodynamic charge transition levels' (aka. defect levels), but do not directly correspond to the electronic states of the defect. This is a subtle distinction, and so I would recommend reading this great review paper to explain the difference: Freysoldt, C.; Grabowski, B.; Hickel, T.; Neugebauer, J.; Kresse, G.; Janotti, A.; Van de Walle, C. G. First-Principles Calculations for Point Defects in Solids. , (1), 253-305. doi.org/10.1103/RevModPhys.86.253. The difference is also readily seen from looking at Figures 1/2 (defect electronic states) and 4 (defect transition levels) in our recent paper on Cd vacancies in CdTe: pubs.acs.org/doi/10.1021/acsenergylett.1c00380 3) Unfortunately no, it was an in-house presentation from a few years ago.

@ 1) Thanks for the clarification. 2) Thanks, I will check out those papers. 3) No worries.

Hi Alief, the KPOINTS file for VASP sets the kpoint grid for the 1st Brillouin Zone of the unit cell structure (POSCAR) you provide. k-spacing refers to the density of this grid, or specifically the distance between these kpoints in reciprocal space. For more on this, I'd recommend looking at the KPOINTS (www.vasp.at/wiki/index.php/KPOINTS) and KSPACING (www.vasp.at/wiki/index.php/KSPACING) vaspwiki pages. Also, the kgrid tool developed by my research group (github.com/WMD-group/kgrid) can be very useful for determining a reasonable kpoint grid for your structure, but as always its best to checkt the convergence with respect to kpoints, which can be done rapidly using a tool I developed here: github.com/kavanase/vaspup2.0 Hope that helps!

Hi! Yes you can calculate the defect formation energy using any atomistic simulation code. CASTEP is quite similar to VASP, using a plane-wave basis set, and so all the same principles discussed here (such as the supercell approach, image charge corrections etc.) will apply. While good for GGA DFT calculations, I have heard that CASTEP is slower than VASP for hybrid DFT calculations, so this will be something to keep in mind.

Hi Samar, more details on these types of calculations can be found in pubs.acs.org/doi/10.1021/acsenergylett.1c00380, pubs.rsc.org/en/Content/ArticleLanding/2022/FD/D2FD00043A and arxiv.org/abs/2207.09862. Basically the idea is to calculate the potential energy surface (PES) as the defect structure changes from the ground-state of one defect charge to the ground-state of the other, giving the PESs shown in the video. We do this by interpolating the structures between these defect charge states, as we would do for a nudged elastic band calculation (you can use the pytmatgen tools for this: pymatgen.org/pymatgen.core.structure.html#pymatgen.core.structure.IStructure.interpolate), then calculate the DFT energies of these structures given the datapoints shown in the figure, which we can then fit a curve to in order to generate the PESs shown. From these we can determine the optical transition energies as discussed in the video and in pubs.acs.org/doi/10.1021/acsenergylett.1c00380, as well as calculating the non-radiative carrier recombination of these defects - which is actually the main limiting factor to the efficiency of emerging solar materials - with more details on this given in the papers I linked above, and the theory is described by Alkauskas et al in link.aps.org/doi/10.1103/PhysRevB.90.075202.

@SeanRKavanagh thanks a lot! I actually attended your MRS talk about shake and break.. Great one!

@@samarfawzy7240 ah I'm delighted to hear that! Thank you very much for the kind words! 😃 Hopefully I'll catch you at a different conference and see your name on some defect papers in the future!

I just want to know how to define energy in CPLAP? Total energy of BaSnO3 is arround 33.70 eV by VASP. How to scale it 11.46 eV?

Hi Alireza, thanks for the compliment! In metallic systems, the situation is somewhat different. Defects will only adopt neutral charge states, as the presence of free electrons will act to neutralise the defect and there is no bandgap for localised defect electronic states to occur. This means that the defect formation energy will be independent of the Fermi level, as is the case here for neutral defects in semiconductors. I'd highly recommend this review paper for more information: Freysoldt, C.; Grabowski, B.; Hickel, T.; Neugebauer, J.; Kresse, G.; Janotti, A.; Van de Walle, C. G. First-Principles Calculations for Point Defects in Solids. , (1), 253-305. doi.org/10.1103/RevModPhys.86.253.

@@romanaafroz1684 you can use either absolute energies from VASP (or your DFT program of choice) or relative formation energies / chemical potentials (i.e. with respect to elemental phases) with CPLAP. As long as you are consistent across your dataset, the program will work as expected and the results will be the same.

@@romanaafroz1684 Alternatively, you can use the pymatgen phase diagram module for this

Hi! Explained now in your comments in the "Metastable Defects..." video: ru-vid.com/video/%D0%B2%D0%B8%D0%B4%D0%B5%D0%BE-HyiUx5HRI4I.html

It depends on the host material and the element forming the vacancy. In semiconductors vacancies can adopt multiple charge states as discussed in this video (and in pubs.acs.org/doi/10.1021/acsenergylett.1c00380). If you use our defects code doped (doped.readthedocs.io), it has a charge state guessing algorithm built in (which we find performs quite well). These charge states are automatically generated when you generate your defects with doped as shown in the tutorial on the website (and you can look at the underlying code for more insight).

Hi! Glad you found it useful! 😃 We try to explain some of this in the video. The defect formation energy depends on the value of the Fermi level (i.e. the x-coordinate in the transition level diagram plots), as shown in the formation energy equation at 8:44 (see the reference given there for more info) and discussed at 18:11. So once we've calculated the other terms in that equation, the Fermi level is the only variable, so we can plot the formation energies against it; discussed by Joe at 24:10. Then the transition levels are simply the intersection points between these lines, as shown at 37:24. In practice, this is usually automated in whatever software you're using to calculate your defects. We use doped (github.com/SMTG-UCL/doped) as our defect code, but there are also others out there

@ Thank you very much for your response. I have a kind request please make a video on the doped package, how can we use it. it will be very convenient for us. if possible

Hi! Yes this is something I was planning to do very soon. We've recently update doped and the documentation, so we have tutorials on using it here: doped.readthedocs.io/en/latest/. But yes I plan on making a tutorial video specifically for doped soon!

Thank you looking forward to hear from you

Hi Joy Roy, glad you found it useful! Sorry about the audio issue, we didn't realise when recording that it would come out fuzzy. Our re-uploaded version linked in the description (ru-vid.com/video/%D0%B2%D0%B8%D0%B4%D0%B5%D0%BE-FWz7nm9qoNg.html) has closed captions, so this might be useful for helping clarify what's being said!

Hi Eric! Yep both tutorials were amalgamated in this video, you can see part 2 starts at 55:56 (where I'm suddenly wearing different clothes 😉)

Thanks for pointing this out Amir! I don't know why that happened, I checked and the audio is fine on the original (and if I download the video from RU-vid 🤔). Re-uploaded here: ru-vid.com/video/%D0%B2%D0%B8%D0%B4%D0%B5%D0%BE-FWz7nm9qoNg.html as RU-vid won't let me edit the video itself.

@ Thank you for taking care of it.

Hi Nigel, for complexes firstly, dealing with chemical potentials is relatively straightforward, and we follow a Law of Mass Action approach, such that the chemical potentials for defect complexes are just the sum of the chemical potentials of the included point defects. For multiple dopants, it's the same approach discussed here, but you would also need to consider competing phases involving both species (e.g. if Cd and Se were co-dopants in your system, you'd need to consider e.g. CdSe as a competing phase). Here I'm assuming you're talking about co-doping. If not and you're just considering multiple potential dopants for a system, you don't need to do anything extra from the approach discussed here as they're independent. For your specific example, yes you'd use the alloy for the bulk properties. For the chemical potentials, you should just treat it as bulk Hg18Cd14Te32 being the host material, with the competing phases to this setting your chemical potential limits. Being a ternary system (Hg-Cd-Te), there are likely to be more competing phases than for binary compounds, but the maths is the same, e.g. µ(Hg18Cd14Te32) = 18µ(Hg) + 14µ(Cd) + 32µ(Te)

@ That was a very prompt reply, thank you. Just a quick follow-up to see if I got this right. I would just follow the same approach as you guys did here for determining the chemical potential dopant limits, for HgTe doped with Cd (replacing the Hg atom). Fortunately, HgCdTe only has HgTe and CdTe as competing phases. Since I would be treating bulk Hg18Cd14Te32 as the host material, the term in 15:19, sum_i n_i mu_i would be equal to +1mu_Hg to remove 1 Hg atom. Is that correct?

@@nigelhew3324 Yes that's correct! And the value of µ_Hg here would be determined by the phase stability region of bulk Hg18Cd14Te32, bordered by HgTe and CdTe, via the simultaneous equations discussed at that point and slightly earlier in the video (14:07). Our defects code doped (github.com/SMTG-UCL/doped) has tools for doing this automatically from VASP calculations, and the PyCDT paper has a nice overview of the formalism for chemical potentials here: www.sciencedirect.com/science/article/pii/S0010465518300079#sec2.3

@ Thank you for clarifying. Yes, I've been using doped recently and find it very helpful.