Questions about the code

[ajgoyal]

Questions:

1. How are the site concentrations calculated by the code?
2. How does the code determine the defect concentration if there are more than one type of point defect on a specific site? For example in NaCl, Vacancy at Chlorine site (v_Cl) as well as a substitutional defect, say Br on Cl (Br_Cl).
3. I believe the underlying assumption of the code is that the defect formation energy does not change with changing temperature and it is fixed to 0K value? If yes, I think it will be useful to state this explicitly either in usage notes and/or README.md description.

[alexsquires]

1. The site concentrations are calculated per-charge state, each concentration is then multiplied by the site degeneracy of defect species, (DefectSpecies.nsites). These are converted to per-volume values only really when reporting values to the user.
2. There are no site competition effects included, defects are expected to be dilute, no defect "sees" any other defects.
3. Yes, this true, and I will include it in the usage notes. Though I believe the code could be fairly straightforwardly patched to do this, maybe for a future version when this analysis becomes more prevalent in the community.

[bjmorgan]

@alexsquires do we say we are using a classic Zhang & Northrup model? (this relates to site competition and temperature independent defect formation energies)

[alexsquires]

No we don't, but we could add this explicity.

[ajgoyal]

Thank you for clarifying, this is helpful. I do have some follow up questions:

1. So as per the paper equation 5 (concentration), there is N{x} (site density) and g{x}^{charge} (charge state degeneracy). I believe in the code DefectSpecies.nsites represent N{x} and DefectChargeState.degeneracy represent g{x}^{charge}?

2. Is the code meant only for vacancies or it can be extended to interstitials and substitutional defects?

3. Yes, I understand that the defect formation energies are calculated under the dilute defect model. However, while solving for the defect concentrations, all the defects are pooled in to find the equilibrium Fermi level. So, if one consider say both vacancy and substitution (whose energies are calculated using dilute model), wouldn't the defects still compete if they form on a same lattice site?

[alexsquires]

1. Yes exactly

2. The identity of the defect is unimportant. No information about the "type" of defect is carried forward to the code, the only per-defect information required by the code is the degeneracies, and their per-charge state formation energies. See for example a recent paper which used py-sc-fermi to calculate the self-consistent Fermi energy which included vacancies, interstitials, antisites and extrinsic dopants: https://pubs.acs.org/doi/full/10.1021/acsenergylett.2c01961 (the SI has the py-sc-fermi methodology details)

3. In reality they will indeed compete, but If the concentration of defects is so low that defect-defect interactions are negligible, then it is also reasonable to assume that defect site competition is negligible.

[ajgoyal]

Thanks for sharing the reference

1. Yes, this is exactly what I was looking for. My request would be if it is not too much trouble for the @authors? To add Ga-doped ZnSb2O6 as an example in the tutorials along with the already given NaCl example. Ga-doped ZnSb2O6 will highlight the capability of the code more than what a simple NaCl example could do.
2. I agree. I would suggest that it would be helpful to mention in the README.md "... defect concentrations under dilute limit at thermal equilibrium ...".

It may help the user to be careful in case code predicts concentrations of the order of the site concentrations ~10^22 cm^-3 at high temperatures. For example, taking the NaCl example from the tutorial, with inputs as below, can give concentrations of vacancies in the order of ~10^23 cm^-3, higher than the site concentration (2.17*10^22 cm^-3).

v_Na_0 = DefectChargeState(charge= 0, degeneracy=1, energy=0.2)
v_Na_minus_1 = DefectChargeState(charge= -1, degeneracy= 1, energy = 0.2)
v_Cl_0 = DefectChargeState(charge= 0, degeneracy= 1, energy = 1.0)
v_Cl_1 = DefectChargeState(charge= 1, degeneracy= 1, energy = -1.0)

v_Na = DefectSpecies(name = "v_Na", nsites = 1, charge_states = {0 : v_Na_0, -1: v_Na_minus_1})
v_Cl = DefectSpecies(name = "v_Cl", nsites = 1, charge_states = {0 : v_Cl_0, 1: v_Cl_1})

dos = DOS(dos=dos, edos=energy, bandgap=4, nelect=20, spin_polarised=False)
defect_system = DefectSystem(defect_species=[v_Na, v_Cl],dos=dos, volume=46.10,temperature=1700)

defect_system.report()
SC Fermi level :      0.6000346708439592  (eV)
Concentrations:
n (electrons)  : 4081163.333937821 cm^-3
p (holes)      : 1.5750771773798547e+20 cm^-3
v_Na           : 3.3837571876937696e+23 cm^-3
v_Cl           : 3.327033507046666e+23 cm^-3

Breakdown of concentrations for each defect charge state:
---------------------------------------------------------
v_Na       : Charge Concentration(cm^-3) Total
           :  0  5.538396e+21          1.64 
           : -1  3.328373e+23          98.36 
---------------------------------------------------------
v_Cl       : Charge Concentration(cm^-3) Total
           :  0  2.353551e+19          0.01 
           :  1  3.326798e+23          99.99 

Thank you!

[bjmorgan]

In reality they will indeed compete, but If the concentration of defects is so low that defect-defect interactions are negligible, then it is also reasonable to assume that defect site competition is negligible.

If you are in a regime where you need to worry about defect competition, you are almost certainly in a regime where you shouldn't be using dilute limit (concentration independent) defect formation energies, and defect–defect interactions will have a much stronger effect on your results (because they affect defect formation energies, which go into the exponent in the defect concentration term).

[alexsquires]

Thanks for sharing the reference

1. Yes, this is exactly what I was looking for. My request would be if it is not too much trouble for the @authors? To add Ga-doped ZnSb2O6 as an example in the tutorials along with the already given NaCl example. Ga-doped ZnSb2O6 will highlight the capability of the code more than what a simple NaCl example could do.

I would feel slightly uncomfortable doing so as I would not want to imply that this work was performed by the developers of py-sc-fermi; on the "homepage" of the documentation there is a section that lists papers (and therefore the systems) that have used py-sc-fermi - would it be acceptable to you if we stated that the NaCl example was a toy-system, and explicitly directed the reader towards the literature for more complex examples?

2. I agree. I would suggest that it would be helpful to mention in the README.md "... defect concentrations under dilute limit at thermal equilibrium ...".

I am happy to do this, yes.

It may help the user to be careful in case code predicts concentrations of the order of the site concentrations ~10^22 cm^-3 at high temperatures. For example, taking the NaCl example from the tutorial, with inputs as below, can give concentrations of vacancies in the order of ~10^23 cm^-3, higher than the site concentration (2.17*10^22 cm^-3).

Yes, I had actually previously raised this as something I would like to implement in the code, I would be happy to expedite its implementation for the review.

[alexsquires]

@ajgoyal after a bit of discussion internally, we've update the DefectSystem object with some more functionality. We've added a method DefectSystem.site_percentages which will give the site percentage of each defect species, in addition the DefectSpecies.report() method will now report the site concentrations as well, i.e.

If you want to check this out, its on the new branch site_percentages.

Hopefully this should help guide the user on the limits on the model.

[ajgoyal]

I would feel slightly uncomfortable doing so as I would not want to imply that this work was performed by the developers of py-sc-fermi; on the "homepage" of the documentation there is a section that lists papers (and therefore the systems) that have used py-sc-fermi - would it be acceptable to you if we stated that the NaCl example was a toy-system, and explicitly directed the reader towards the literature for more complex examples?

In my opinion having a simple but concrete example then just a toy example is more helpful for the learners. Sometimes, papers do the site the tool but they do not give sufficient details to reproduce the results. But, anyway, I will leave up to the @authors. Thank you!

[ajgoyal]

@ajgoyal after a bit of discussion internally, we've update the DefectSystem object with some more functionality. We've added a method DefectSystem.site_percentages which will give the site percentage of each defect species, in addition the DefectSpecies.report() method will now report the site concentrations as well, i.e.

If you want to check this out, its on the new branch site_percentages.

Hopefully this should help guide the user on the limits on the model.

Thanks for making the changes. Yes, this will be helpful.

[rkurchin]

Mentioning the ongoing JOSS review in this issue to easily keep tabs on it (though it looks like this will probably be closed soon?).