Tutorial: Modelling Point Defects in Semiconductors with VASP (Audio Fix)

Citable DOI: 10.5281/zenodo.10981906
Reuploaded due to KZbin error in audio/video sync in final 30 mins; original video here: • Tutorial: Understandin...
*Now with bonus 25 minutes intro to defect structure searching and ShakeNBreak*
Tutorial talk on understanding and computationally modelling defects in semiconductors (using plane wave DFT with VASP). Given jointly with Joe Willis at the Scanlon Materials Theory Group meeting(s).
doped: doped.readthed...
ShakeNBreak: shakenbreak.re...
See our open-access papers on modelling semiconductor defects here:
iopscience.iop...
doi.org/10.102...
ShakeNBreak: www.nature.com...
Slides available here:
speakerdeck.co...
speakerdeck.co...
For other talks on KZbin, have a look at my channel!
For other research articles see:
bit.ly/3pBMxOG
Slides for my talks can be found here:
speakerdeck.co...
Other references:
Matter Preview of Defect Structure Searching: www.sciencedir...
Metastable defects : doi.org/10.103...
Recombination at V_Cd in CdTe (case study): pubs.acs.org/d...

[@SolStryker]

1:32:10 Baseline params for defect calculations 1:39:32 Introducing charged states 1:43:19 Testing calculation params with supercell 1:49:53 Finite-size corrections 1:57:34 Polarons

[@AbdurRahim-hi6un]

Awesome. Plz make more videos about VASP. Specially, charge transfer mechanism, Bader charge analysis

[@mubbasilkhanpathan3665]

Hello I am currently pursuing my masters in physics. My main focus is on thin film solar cell simulations in SCAPS and all the basic information for defects and chemical potential was a great help thanks for this video

That's great to hear Mubbasilkhan!

[@johnlee3122]

13:22 Would you be kind enough to explain how you got the -8.0 for the E^O ? I'm currently trying to calculate energy of creating surface oxygen vacancy on ceria, where I'm kinda stuck on how I should calculate the energy of single oxygen element. :(

[@MRF77]

2:01:08 Here we see the hole charge density because of how you setup the magmom for V_In case. But what if your metal is ferromagnetic? i.e. let's say you have Fe2O3 instead of In2O3, in this case do you still set magmom for Fe atoms to be 0 (in order to show charge state of atoms surrounding defect)? 1:29:30 I'm very bad at chemical intuition! So it's hard interpret where there are 7 combinations for each substitution, which one to potentially rule out! Is there perhaps a "guaranteed"/methodical approach to rule out, when unsure based on intuition? 1:39:30 I'm not sure which tag represents the number of valance electrons of an atom in POTCAR (to calculate NELECT)? Could you please specify? Lastly, let's say we have a neutral In vacancy (V_In) an In2O3 (as opposed to V_O). Q1) Then the vacancy charge should have a charge of 0-(+3) = -3, right? Q2) Now If we want to dope it with 3 Li atoms at each V_In site to preserve the charge neutrality (i.e. substitution), would you suggest testing the formation energy for 1 Li and 2 Li atoms in the V_In site too? If I substitute let's say only 1 Li+ atom, do I need to create 2 holes (by Ploaron method) using the surrounding Oxygen atom? Thank you so much both for all your effort!

Hi Abdullah. For Q.1, if your material is ferromagnetic, you would set MAGMOM (and NUPDOWN) according to the expected spin configuration of all the magnetic atoms in your system, so this would not be 0 for Fe in this case. You would have MAGMOM as e.g. 1 for each Fe atom (e.g. low-spin Fe3+), matching the bulk groundstate spin configuration, and then from this reference point we adjust MAGMOM to match our expected defect spin configuration (e.g. hole on Oxygen -> change MAGMOM of one of the neighbouring O atoms). With this result we will still be able to visualise the hole charge density as shown in the video, by generating the VASP PARCHG for this hole band. (www.vasp.at/wiki/index.php/Band_decomposed_charge_densities). This is actually the case for a recent study of ours on defects in the LMNO battery cathode material, which is ferrimagnetic: chemrxiv.org/engage/chemrxiv/article-details/63d2a2c41fb2a8767ee1e06f Q.2: There's not really any 'guaranteed' method for choosing the charge states to calculate, beyond just calculating all potential charge states, and then seeing which ones end up being stable in the bandgap when we finally parse the results and plot our defect formation energy diagram. There's no harm in doing additional charge states that end up being unstable, beyond the fact that it uses up some computational resources. For picking charge states, it's best to be safe and include more rather than less charge states (i.e include charge states that we think are possible yet unlikely). E.g. we find some "unusual" charge states to be stable in Sb2Se3 due to a strong ability to rearrange and form new bonds: arxiv.org/abs/2302.04901 The standard approach is to rationalise likely defect charge states in terms of the oxidation states of the elements involved, and the Madelung potential of the defect site (basically if it's a cation site it will favour more positive charge states, and vice versa for an anion site). For e.g. Cd_on_Te, Cd likes to be +2 and Te in CdTe is a -2 site, so the 'fully-ionised' charge state is +4 (+2 on a -2 site). So we would want to calculate from 0 to +4. For Te_on_Cd, Te can be -2, +2, +4 and +6, and is now on a cation site so will favour more positive charge states. So fully-ionised is now -4 (-2 on a +2 site), but Te could also be +4 (which is a defect charge state of +2), but very unlikely to be +6 (very extreme charge state), so we'd want to calculate from -4 to +2. I've done the calculations for these defects, and it turns out that Cd_on_Te is stable from 0 to +2, and Te_on_Cd from -2 to +2. So in each case the actual number of stable charge states is slightly less than our test range (two less in each case), but was good to be safe and include all those potential charge states to be sure. It's summed up well in this paper (iopscience.iop.org/article/10.1088/2515-7655/aba081): "The range of possible charge states is usually inferred from the oxidation states of atoms involved in the defect. For example, the charge state of SnZn can be 0 or +2 which correspond to Sn(II) and Sn(IV), respectively. However, the chemistry occurring at defect sites can be unexpected, e.g. cation-cation bonding, so an unbiased search over a wide range of charge states is often necessary to identify the accessible configurations." This is something that comes with experience! From doing the calculations yourself and from reading other papers. The charge state range output by doped is usually a pretty good start.

For Q.3; the valence electrons of each atom in the POTCAR pseudopotential are given on the second line of the POTCAR for each atom (i.e. the line under e.g. " PAW_PBE Cs_sv 08Apr2002"). This is automatically parsed by doped to then automatically determine NELECT and NUPDOWN for the INCAR. For the final Q; well the 'neutral' vacancy means it has a charge of 0. In the fully-ionised charge state however, it would be -3 as you say. (So for this defect we'd want to calculate all charge states from 0 to -3). If you were looking at the possibility of placing Li+ on the V_In site, then yes I'd definitely recommend calculating this sequentially from one to two to three Li. In that case, yes it would be good to explicitly set the location of the two holes on the surrounding oxygen with MAGMOM, to ensure VASP correctly localises the charge. Just to mention, I think it's quite unlikely that 3 Li ions will be stable on the V_In site, as that's a lot of atoms to try fit in to a very small space!

[@MRF77]

​@ Hi Sean, thanks very much for the comprehensive reply and sharing your paper! I'm still trying to figure out actually calculating chemical potential terms, where there's competing phases (like in In2O3, coz how do we even know beforehand those competing phases exist 13:57 ! dumb question, I know!). BTW, the github link in the slide at 14:50 doesn't seem to work. Please check whenever possible. From my last question though, I was referring to 6:22 where the charge of V_Cd2- is perhaps referred to as q = 0 - (+2) = - 2 for vacancy site of the Cd2+ ion. In this case, in order to make the overall charge neutral, two holes need to be created in the surrounding atoms of V_Cd2- site, but I wonder which atoms should I pick (out of four) spin-up (u) in magmom? Is there any preference or any of the two atoms would do? e.g. *duuu* resulting two net spin-ups. Now let's say our vacancy charge is positive, V^q = +2, and we'd want to make two electrons (as opposed to holes) instead. Should we follow similar approach with magmom, but with spin-down (d), e.g. *uddd* resulting two net spin-downs?

​@@MRF77 Hi Abdullah, for figuring out the chemical potential terms, we have to calculate the energies of potential competing phases with the same elements (e.g. compounds containing Fe and/or O for Fe2O3). We usually obtain these potential competing phases from a database, for example in doped it does this automatically for you by querying the Materials Project data and returning the possible competing phases (with a certain error threshold on the Materials-Project-calculated energies), see the example here: github.com/SMTG-UCL/doped/blob/master/examples/dope_chemical_potentials.ipynb The GitHub link you mentioned is for an internal wiki page (sorry this was initially recorded as an internal lecture only. I will reply to your email with the PDF of this page though). About V_Cd^-2, yes it's minus 2 charged in the _fully-ionised charge state_, however defects usually have multiple possible charge states, see e.g. Fig. 1 in pubs.acs.org/doi/10.1021/acsenergylett.1c00380. If we have V_Cd^-2 in our material, you're right that there needs to be some positive charge elsewhere in the material to counter-balance this and make it charge neutral. However, this does not have to be in the immediate surrounding environment of the defect, and could come from delocalised holes far away from the defect, or in fact other _positively-charged_ defects in the same material (e.g. Cd_i^2+ in CdTe). When you perform the defect supercell calculation in VASP or other electronic structure packages, it adds a neutralising opposite-charge-density background if you have a charged defect like this (briefly discussed at 1:37:22, more discussion in journals.aps.org/rmp/abstract/10.1103/RevModPhys.86.253. In the neutral charge state of the defect (V_Cd^0), this can be thought of as 'adding' two holes to our ionised V_Cd^-2 defect, and there are a couple of possibilities of for where these could be placed. We find that you can have a bipolaron with two isolated holes on two of the nearby Te atoms, or a Te-Te dimer can form, as discussed in: pubs.acs.org/doi/10.1021/acsenergylett.1c00380. If we know we have a hole-polaron solution for our defect, then often we have to enumerate the different possibilities of where these holes could be located, and compare their energies (i.e. testing each of **duuu**, **dudu** etc). This is relevant to our recent work on defect-structure searching: www.nature.com/articles/s41524-023-00973-1 where we developed the ShakeNBreak package (shakenbreak.readthedocs.io/en/latest/) to help automate this ground-state-searching problem.

[@asiyeshokri9544]

thank you so much for great explanations. My question is, by this correction, this way of ploting formation energy vs fermi energy is again connected to how much persent of defect we have in our crystal or not?

Thanks Asiye! Yes with this formation energy plot, we can firstly determine the self-consistent Fermi level in the material (discussed around 37:32), which then gives us the formation energies of defects in the material, which we can then use to calculate their concentration (i.e. percentage defects) using the equations shown at 34:13. In practice, we usually automate these calculations, and most defect computational packages will perform this analysis for you. The defect packages we discuss here are our doped (github.com/SMTG-UCL/doped) and ShakeNBreak (shakenbreak.readthedocs.io/en/latest/) tools. With doped you can perform this analysis, and do further processing and analysis of the Fermi level and defect concentrations with the py-sc-fermi package: py-sc-fermi.readthedocs.io/en/latest/

[@asiyeshokri9544]

@ Thank you so much for kind reply!

[@thebhamuji]

Wonderful

[@samarfawzy7240]

Thank you!

[@SelmaMayda]

At 1:39:21, you showed the electron number of In2O3 with oxygen vacancy. I could not understand the numbers. How many oxygen atoms did you remove from the system? It should be 2 oxygen atoms wrt numbers. Am I right?

Hi @SelmaMayda, no, here we're removing one oxygen atom (i.e. one oxygen vacancy) from the system. The absolute NELECT (number of electrons) in the calculation here is somewhat arbitrary, as it depends on both the POTCAR (pseudopotential) and supercell choices. So it's really the change in NELECT here that we're trying to exemplify; that changing the number of electrons in the supercell with NELECT is how you can change the defect charge state. So adding one electron (by adding +1 to NELECT) is adding an electron to the supercell, going from the neutral vacancy to the -1 charge vacancy, etc

[@harishabibi6324]

Thank you so much sir for such a informative video, can you please make a video on "how to calculate thermoelectric properties with boltztrap on VASP PLEASE"

Hi Haris, thanks for your nice comment! Unfortunately I'm not an expert on using Boltztrap though!

[@harishabibi6324]

@ it's oky sir, thank you.. I have one question. if you have time, can you please make a video on "how to calculate young modulus, elastic constant with VASP please for 2D materials

[@asiyeshokri9544]

Hello, thank you again for your great video. can you please explain how to have an optimization choose for nodes and cores on a computing machin? emagine having 160 atoms in a bulk and having HSE calculations to get the formation energy! I am using 24 nodes each having 24 cores with the help of input file seting you explaind. . best.

Hi Asiye! Determining the optimal number of nodes and cores for a particular calculation can be a bit tricky and depends on several factors, including the size of your system, the type of calculation you're running, and the specific architecture of your computing cluster. As a rule of thumb, larger systems and more computationally intensive calculations (like those involving hybrid functionals such as HSE) will benefit from using more nodes and cores. However, there is a limit to how much additional performance you can get from using more computational resources due to communication overhead between nodes. Given your system size of 160 atoms and the use of HSE calculations, using 24 nodes each with 24 cores sounds like a reasonable starting point. The most important INCAR parameters for calculation optimisation are NCORE and KPAR. NCORE depends on your supercomputer architecture, and usually for 24-core nodes, NCORE = 12 is best. Then KPAR depends on the number of kpoints for your calculation, and should divide into your number of kpoints. Usually for a HSE defect supercell calculation, we would use a value between 1 and 4 for KPAR, depending on the number of kpoints. Our defects calculation package doped (github.com/SMTG-UCL/doped) will automatically generate the defect structures and VASP input files for your defect calculations (as well as parsing & plotting the results), so I would recommend using this! Best of luck with your calculations!

[@NitinKMSP]

At 24:48, Is it formation energy or enthalpy? What is the difference between the two?

[@NitinKMSP]

In this video, it is given that, ΔH(D,q) = E(D,q) − E(H)+ Σniμi + q*E(F) + Ecorr In an article I find a similar expression,ΔEd(D,q) = E(D,q) − E(bulk) − ΣniEi − ΣniΔμi + q*μ(e) + Ecorr. Comparing the two, I have 3 doubts 1) How E(H) and E(bulk) are related? 2) How ΔH(D,q) and ΔEd(D,q) are related? 3) At 24:48, why didn't you include the enthalpy/energy of impurities added?

[@NitinKMSP]

Please help me. I am struggling in understanding this

Hi Nitin, yes this point can be confusing due to mixed terminology in the field. The two formulas you give both refer to the defect formation enthalpy, however often this is referred to as the defect formation 'energy' in the literature. 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! For your specific Qs: Those two equations you give are the exact same, just with different symbols/representations of the same terms. 1. They're the same (i.e. E(H) = E(Host) = E(Bulk)) 2. They're the same 3. I'm not totally sure what you're asking here, we do include the enthalpy of impurities added by computing the E(Defect supercell) - E(Host supercell) term, and accounting for the chemical potentials of any additional/removed atoms with Σniμi Hope this helps! 😀

[@nitink9879]

@ Hi Sean, thank you very much for clearing my doubt❤ At 3) I was talking about the extra term − ΣniEi in the equation in the article I was talking about (Reference: Spinney: Post-processing of first-principles calculations of point defects in semiconductors with Python). Also, they have used chemical potential of electron instead of fermi energy E_F which is not always valid, because chemical potential is equal to fermi energy only at 0 K.

Hi Nitin, you're welcome! 😃 So in those equations; these two terms are equal: Σniμi = ΣniEi + ΣniΔμi It's just a different representation of the same quantities. In the Spinney article they use ΣniEi as the reference energies of the elements and then adding the chemical potential differences with respect to these references (ΣniΔμi) to give the same energies as Σniμi

[@BabiLectures]

Is this code compatible with Quantum espresso? If No what effort is put in place to make it compatible?

Hi! The defect structure generation, and plotting/analysis/thermodynamics functions are all agnostic to the underlying DFT/forcefield code. However full direct calculation I/O is currently only supported for VASP. So you can generate your defect supercell and structure inputs for Quantum Espresso, but then the user is required to write the other input files controlling basis set etc. For parsing the relaxed structures and energies, the user should be able to do this without much issue, if they have some moderate python skills. We are hoping to make it more compatible with other codes in the future, but this may be a good while away.

[@BabiLectures]

@ Thanks for your response