Hi! So most parts of the doped (and ShakeNBreak) codes are agnostic to the underlying DFT/ML code being used, so you can still use it to generate defect structures, guess appropriate charge states, analyse symmetries, analyse relaxed structures, create DefectEntry objects and analyse defect thermodynamics etc. Doped uses pymatgen objects under the hood for the structures, so you can use its IO functions to output to CIF/XYZ etc which can be easily converted into/from Wien2k structure files. Hope that helps!

Hi! This question is discussed here in the doped Issues: github.com/SMTG-Bham/doped/issues/64 In a disordered material these sites are all symmetry-inequivalent. You can't just calculate e.g. one vacancy site for a disordered system.

Hi @maram-uL, thanks for your interest! Yes ShakeNBreak can be used with slab and cluster models, though there may be a couple additional considerations, such as manually setting the oxidation states of the elements (ShakeNBreak will try guessing these automatically and will print this info for you, but for such materials they are more likely to be off). Also, for dislocations (extended rather than point defects), the typical bond distortions employed in ShakeNBreak will not work directly out of the box. You can use the ShakeNBreak functions to rattle these structures, which already gives a better chance of finding the ground-state arrangements (works in ~30-40% of known reconstruction cases). Alternatively you could look at modifying the code somewhat to apply these bond distortions to multiple atoms around the dislocations? Hope that helps!

Thank you for your response. I'll definitely go through the code and return with more specific questions (through other channels) to adapt it for my case study.

Hi! To get the oxygen chemical potential, you need to calculate oxygen in its standard state (i.e. as O2 gas, not an isolated O atom). Typically this is done by calculating an O2 molecule in a supercell. This is automatically implemented in doped, where it generates this structure and the corresponding input files for you. Would recommend looking at our chemical potentials tutorial here: doped.readthedocs.io/en/latest/chemical_potentials_tutorial.html

Thank you so much :)) have a great day

Hi, yes lm-decomposed orbital projections for the band structure / DOS are exemplified here on the docs: smtg-bham.github.io/easyunfold/examples/example_nabis2.html#lm-decomposed-orbital-projections

Thank you very much for your response! I understand how to plot the d orbital projections from the tutorial. When I use the command 'easyunfold unfold plot-projections --atoms "Cr, Cr" --orbitals="dxy,dyz,dxz|dz2, dx2" ' i get the pDOS for each of the projections. But what I really want is THE pDOS of the sum of the dxy,dyz,dxz (t2g orbitals) and the sum of the dz2,dx2 (eg orbitals) in the same graph. Can you tell me please if this can be done? I don't see it clear in the tutorial. Thank you very much for your valuable contribution.

​@@luisfelipeleonpinzon4622 ah I see. This isn't automatically implemented in the code at present. In theory it should just be a small update to the orbital selection part of the code (if the user provides "t2g", then grouping the dxy/dxz/dyz orbital contributions and labelling as t2g), so if you have a little python experience maybe you could try this? If so you could maybe open a pull request on the GitHub for this? Or an Issue if you're not sure how to implement (the main developers are quite busy these days...)

Hahaha never copped the resemblance 😂 Thanks!

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 Harish, the `icet` code have a good background and tutorial page on doing this here, would recommend checking it out! icet.materialsmodeling.org/advanced_topics/sqs_generation.html Also the data upload linked in the paper (pubs.acs.org/doi/10.1021/acs.jpcc.3c05204) has the notebooks I used to generate the SQS structures. Hope that helps!

thanks. i will check it out

this link is showing error. I want to generate the SQS of Ca(Ti0⋅25Zr0.25Hf0.25Sn0.25)O3. can you guide me how to generate this sqs?

Hi! Currently easyunfold is only directly compatible with VASP and CASTEP. You could try updating the functions slightly to read in Quantum Espresso wave function data for it to work, as it's just an I/O issue

Hi Thomas, yes it was mostly similar to this. Unfortunately time constrained so just a quick overview of the code. Best to check out the docs for more details! I'll try make time to do a more in depth KZbin tutorial of the code soon 🤞

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 for parsing the relaxed structures and energies from CASTEP, this is doable but requires some manual python handling from the user, to create the doped DefectEntry objects for the defect supercells. For this, the CASTEP interfaces with ASE (and thus pymatgen) would be useful! We are hoping to make it more compatible with other codes in the future, but this may be a good while away.

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.

@ Thanks for your response

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

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 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!

Hi Salem, all our data and input files for this project are available on the open-access repository Zenodo (zenodo.org/record/4683140#.ZFYlSuyZNhE) as well as Materials Cloud (archive.materialscloud.org/record/2021.133). Hope that's useful! 😃

Hi! Thanks for pointing this out! We noticed this issue recently (due to pymatgen's default oxidation state functions not returning a value for Au), and implemented a fix in the *develop* branch of ShakeNBreak on GitHub (github.com/SMTG-UCL/ShakeNBreak/compare/a69bb5a3d88f...f3064f08556e). This will be merged to the main brach and updated on PyPI (for standard installation: pip install ShakeNBreak) in the next few days (we're just waiting to add a link to a Nature Physics News & Views article on SnB that we've been told will be published this week). For now, you can install the current *develop* version of SnB with: git clone github.com/SMTG-UCL/ShakeNBreak cd ShakeNBreak git checkout develop pip install -e . Alternatively if you can wait a few days, I'll reply to this comment again when the PyPI version has been updated, and then you'll be able to update just with: pip install --upgrade ShakeNBreak

Hi! Just to let you know, we've now released this updated version of ShakeNBreak with this fix (and some other performance improvements), so it can now be installed using pip install --upgrade shakenbreak. Thanks again for noting this! 😃

@ Thanks very much for your response!

Hi Dark! Thanks for the kind words and sorry for the delay in replying, I was away at MRS Spring and on holidays after. For Quantum Espresso, we have examples of the defect structure generation here: shakenbreak.readthedocs.io/en/latest/ShakeNBreak_Example_Workflow.html#optional-input-file-generation-for-other-codes and at: shakenbreak.readthedocs.io/en/latest/Generation.html. The parsing and analysis parts of the code should be independent of the code used. Let us know if you have any issues with this! 😃

@ Thank you very much.

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!

​@ 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.

Thanks @MA Fadla! 😃

Hi! Sorry for the slow reply, I was at APS March + hols the last 2 weeks. Hmmm, so for the formula n = g*N*exp(-ΔH/KbT), usually we use this after post-processing our defect calculations to get the formation energy ΔH of a **single point defect**. Then N is the density of available sites for that defect (so e.g. for an oxygen vacancy this would be the density of oxygen atomic sites in the structure) - most defect calculation tools like our doped code (github.com/SMTG-UCL/doped) handle this automatically, and g is the configurational degeneracy of the relaxed defect structure. We have a discussion on this in the 'Tei in CdTe: site degeneracies' section of: pubs.rsc.org/en/Content/ArticleLanding/2022/FD/D2FD00043A.

For a supercell with multiple vacancies, it depends on what exactly you're trying to calculate. Is it a complex defect with V_Mo + 3*V_O? There will be a major entropic cost to forming such a defect complex involving 4 point defects (which should be included in ΔH), but the same formalism would apply with N being the density of available sites for this defect complex, and g the configurational degeneracy of the relaxed structure for this complex defect

Or is it a sub-stoichiometric phase (i.e. with vacancy concentrations on the order of a few percent)? In this case, we are no longer dealing with the dilute limit of point defects, and so the formalism for calculating this case differs from that of isolated point defects

Hi! Sorry for the slow reply, I was at APS March + hols the last 2 weeks. Hmmm, so for the formula n = g*N*exp(-ΔH/KbT), usually we use this after post-processing our defect calculations to get the formation energy ΔH of a **single point defect**. Then N is the density of available sites for that defect (so e.g. for an oxygen vacancy this would be the density of oxygen atomic sites in the structure) - most defect calculation tools like our doped code (github.com/SMTG-UCL/doped) handle this automatically, and g is the configurational degeneracy of the relaxed defect structure. We have a discussion on this in the 'Tei in CdTe: site degeneracies' section of: pubs.rsc.org/en/Content/ArticleLanding/2022/FD/D2FD00043A.

For a supercell with multiple vacancies, it depends on what exactly you're trying to calculate. Is it a complex defect with V_Mo + 3*V_O? There will be a major entropic cost to forming such a defect complex involving 4 point defects (which should be included in ΔH), but the same formalism would apply with N being the density of available sites for this defect complex, and g the configurational degeneracy of the relaxed structure for this complex defect =

Or is it a sub-stoichiometric phase (i.e. with vacancy concentrations on the order of a few percent)? In this case, we are no longer dealing with the dilute limit of point defects, and so the formalism for calculating this case differs from that of isolated point defects

Hi! Explained now in your comments in the "Metastable Defects..." video: kzbin.info/www/bejne/fqrMhqtrfbd8aqs

Hi! Sorry I didn't see this before, it got caught in KZbin'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 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 :)

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?

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! 😀

@ 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

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/

@ Thank you so much for kind reply!

Thanks @Hashtag Peace! 😃

Wow! Trippy idea! I can imagine there'd be a lot of interesting applications of this nonlocal info/energy/signal transfer... I guess it would depend on how the separation of the entangled electrons affects the quasiparticle interaction between them; e.g. might work easier for superconducting electron pairs rather than (Frenkel) excitons which rely on a real-space Coulomb interaction?

That's great to hear Mubbasilkhan!

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

@ 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

Thanks @Glaedr_4_life, really glad to hear it! 😃

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!

when you say DFT energies you mean through scf calculations?

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 KZbin 🤔). Re-uploaded here: kzbin.info/www/bejne/fIjdaKGjbtaihMk as KZbin 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.

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 😉)

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 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 (kzbin.info/www/bejne/fIjdaKGjbtaihMk) has closed captions, so this might be useful for helping clarify what's being said!

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.

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