unable to calculate the chemical potential ranges for extrinsic dopant

[Atefeh-Yadegarifard]

Hello. I have been trying to dope CsPbI3 with Sn at the B-site. I have generated the competing phases with doped using e_above_hull = 0, which resulted in ['Sn', 'SnI4', 'Cs2SnI6']. After relaxing the structures, and trying to parse the data if I put only Cs2SnI6_EaH_0, Sn_EaH_0 and SnI4_EaH_0 in the competing_phases folder, the code asks for the relaxation data of the parent bulk phase as: "Could not find bulk phase for CsPbI3 in the supplied data. Found phases: []" However, if i put the vasp relaxtion result of CsPbI3_bulk into the competing_phases folder, this error occurs:

``--------------------------------------------------------------------------- KeyError Traceback (most recent call last) Cell In[2], line 2 1 cpa = CompetingPhasesAnalyzer("CsPbI3", extrinsic_species = "Sn") ----> 2 cpa.from_vaspruns(path='./competing_phases/', 3 folder = "vasp_std", 4 csv_fname = "./competing_phases/cspbi3_Sn_doped_energies.csv") 5 df= cpa.calculate_chempots (csv_fname = "./competing_phases/cspbi3_Sn_chempots.csv")

File ~/miniconda3/lib/python3.11/site-packages/doped/chemical_potentials.py:1199, in CompetingPhasesAnalyzer.from_vaspruns(self, path, folder, csv_fname) 1190 d = { 1191 "formula": v["pretty_formula"], 1192 "kpoints": kpoints, (...) 1195 "energy": final_energy, 1196 } 1197 temp_data.append(d) -> 1199 formation_energy_df = _calculate_formation_energies(temp_data, self.elemental_energies) 1200 if csv_fname is not None: 1201 formation_energy_df.to_csv(csv_fname, index=False)

File ~/miniconda3/lib/python3.11/site-packages/doped/chemical_potentials.py:181, in _calculate_formation_energies(data, elemental) 179 for d in data: 180 for e in elemental: --> 181 d[e] = Composition(d["formula"]).as_dict()[e] 183 formation_energy_df = pd.DataFrame(data) 184 formation_energy_df["formation_energy"] = formation_energy_df["energy_per_fu"]

KeyError: 'Sn'

I would like to ask how to solve this issue. Thank you.

[kavanase]

Hi @Atefeh-Yadegarifard, thanks very much for flagging this! 😄 This was occurring due to an update in pymatgen. I've fixed it now and will be included in the new version of doped to be released later today/tomorrow. I'll notify you when that's done!

[kavanase]

Hi @Atefeh-Yadegarifard, thanks again for notifying us of this 😄 It's now been fixed in the new version of doped (2.1.0, just released), which you can install with pip install --upgrade doped. Please reopen an issue if you run into any problems!

[Atefeh-Yadegarifard]

Dear Sean, Thanks for the update. My problem in calculating the chemical potential ranges was solved. However, now the package shows another bug. Previously I was able to successfully parse my defect calculations and generate a transition level diagram. However, after the update it is not working and I am seeing this error: "The lattice constants for defect and perfect models are different", even though the supercells used for all calculations have the same size. I will be thankful if you can help me fix this error as well. Thanks

[kavanase]

Hi Atefeh, Sorry about that issue! I haven't seen this pop up on our tests. Would you be able to share some defect & bulk vasprun.xml(.gz) so I can reproduce this and fix? Otherwise you could maybe send some screenshots or copy and paste the error traceback? Thanks!

[Atefeh-Yadegarifard]

Hi Sean, Thanks for the reply. This is 3 vasp runs from bulk and defect calculations: vaspruns.tar.gz

and this is the error traceback :

Parsing SnPb-2...

SupercellError Traceback (most recent call last) Cell In[3], line 7 5 print(f"Parsing {i}...") 6 defect_path = f"./cspbi3_doped_various_charges_LOCPOT/{i}/vasp_std" ----> 7 Sn_Pb_dict[i] = analysis.defect_entry_from_paths(defect_path, bulk_path, dielectric)

File ~/miniconda3/lib/python3.11/site-packages/doped/analysis.py:690, in defect_entry_from_paths(defect_path, bulk_path, dielectric, charge_state, initial_defect_structure, skip_corrections, bulk_bandgap_path, **kwargs) 687 skip_corrections = True 689 if not skip_corrections: --> 690 dp.apply_corrections() 692 # check that charge corrections are not negative 693 summed_corrections = sum( 694 val 695 for key, val in dp.defect_entry.corrections.items() 696 if any(i in key.lower() for i in ["freysoldt", "kumagai", "fnv", "charge"]) 697 )

File ~/miniconda3/lib/python3.11/site-packages/doped/analysis.py:1454, in DefectParser.apply_corrections(self) 1449 # try run Kumagai (eFNV) correction if required info available: 1450 if ( 1451 self.defect_entry.calculation_metadata.get("bulk_site_potentials", None) is not None 1452 and self.defect_entry.calculation_metadata.get("defect_site_potentials", None) is not None 1453 ): -> 1454 self.defect_entry.get_kumagai_correction(verbose=False, error_tolerance=self.error_tolerance) 1456 elif self.defect_entry.calculation_metadata.get( 1457 "bulk_locpot_dict" 1458 ) and self.defect_entry.calculation_metadata.get("defect_locpot_dict"): 1459 self.defect_entry.get_freysoldt_correction(verbose=False, error_tolerance=self.error_tolerance)

File ~/miniconda3/lib/python3.11/site-packages/doped/core.py:413, in DefectEntry.get_kumagai_correction(self, dielectric, defect_outcar, bulk_outcar, plot, filename, return_correction_error, error_tolerance, kwargs) 407 if dielectric is None: 408 raise ValueError( 409 "No dielectric constant provided, either as a function argument or in " 410 "defect_entry.calculation_metadata." 411 ) --> 413 efnv_correction_output = get_kumagai_correction( 414 defect_entry=self, 415 dielectric=dielectric, 416 defect_outcar=defect_outcar, 417 bulk_outcar=bulk_outcar, 418 plot=plot, 419 filename=filename, 420 kwargs, 421 ) 422 correction = efnv_correction_output if not plot and filename is None else efnv_correction_output[0] 423 self.corrections.update({"kumagai_charge_correction": correction.correction_energy})

File ~/miniconda3/lib/python3.11/site-packages/doped/corrections.py:420, in get_kumagai_correction(defect_entry, dielectric, defect_outcar, bulk_outcar, plot, filename, verbose, style_file, kwargs) 412 bulk_supercell.remove_oxidation_states() # pydefect needs structure without oxidation states 413 bulk_calc_results_for_eFNV = CalcResults( 414 structure=bulk_supercell, 415 energy=np.inf, 416 magnetization=np.inf, 417 potentials=bulk_site_potentials, 418 ) --> 420 efnv_correction = make_efnv_correction( 421 charge=defect_entry.charge_state, 422 calc_results=defect_calc_results_for_eFNV, 423 perfect_calc_results=bulk_calc_results_for_eFNV, 424 dielectric_tensor=dielectric, 425 defect_coords=defect_entry.sc_defect_frac_coords, # relaxed defect coords (except for vacancies) 426 kwargs, 427 ) 428 kumagai_correction_result = CorrectionResult( 429 correction_energy=efnv_correction.correction_energy, 430 metadata={"pydefect_ExtendedFnvCorrection": efnv_correction}, 431 ) 433 if verbose:

File ~/miniconda3/lib/python3.11/site-packages/pydefect/cli/vasp/make_efnv_correction.py:39, in make_efnv_correction(charge, calc_results, perfect_calc_results, dielectric_tensor, defect_coords, accuracy, unit_conversion) 23 def make_efnv_correction(charge: float, 24 calc_results: CalcResults, 25 perfect_calc_results: CalcResults, (...) 28 accuracy: float = defaults.ewald_accuracy, 29 unit_conversion: float = 180.95128169876497): 30 """ 31 Notes: 32 (1) The formula written in YK2014 need to be divided by 4pi in the SI unit. (...) 36 to make potential in V. 37 """ 38 sites, rel_coords, defect_coords = \ ---> 39 make_sites(calc_results, perfect_calc_results, defect_coords) 41 lattice = calc_results.structure.lattice 42 ewald = Ewald(lattice.matrix, dielectric_tensor, accuracy=accuracy)

File ~/miniconda3/lib/python3.11/site-packages/pydefect/cli/vasp/make_efnv_correction.py:65, in make_sites(calc_results, perfect_calc_results, defect_coords) 63 def make_sites(calc_results, perfect_calc_results, defect_coords): 64 if calc_results.structure.lattice != perfect_calc_results.structure.lattice: ---> 65 raise SupercellError("The lattice constants for defect and perfect " 66 "models are different") 67 structure_analyzer = DefectStructureComparator( 68 calc_results.structure, perfect_calc_results.structure) 69 if defect_coords is None:

SupercellError: The lattice constants for defect and perfect models are different

[kavanase]

Hi Atefeh, Thanks for sharing this. So this error message is correct, your lattice constants are significantly different for your bulk and defect supercells (nearly a 10% volume change):

This causes issues with the FNV/eFNV finite-size charge correction schemes, as for these the bulk and defect supercells should have the same lattice constants. This is why in doped the default INCAR settings have ISIF = 2 (i.e. keeping the supercell shape fixed), which is what we almost always do for defect calculations (see e.g. the YouTube tutorial).

If you want to use the finite-size charge corrections, I think it would be best to recalculate these defects with fixed supercells as is the norm (you can do this quickly using the functions shown in the doped generation tutorial)

[Atefeh-Yadegarifard]

Dear Sean, Thanks for checking my calculations. I will run them one more time and try again :)

[Atefeh-Yadegarifard]

Hello again! I have been trying to calculated the defect formation energy of the CsPbI3 doped with Sn. I was able to calculate the chemical potential ranges with doped. However, some 0 values were returned in the limit and it caused an error in calculating the defect formation energy graphs. I would like to ask how this error can be avoided ? the details of the error and the calculations are in the pdf file attached here. Thank you. chemical potential range issue.pdf

[kavanase]

Hi Atefeh, The chemical potential limits with values of 0 are fine, this just means that chemical potential facet has an element in equilibrium at that point in the phase diagram / chemical potential limit – in this case metallic lead (Pb). The issue with the plotting seems to be related to the format of your input chemical potentials (chemlims). I can't see the form of this input in the PDF you've attached. It should match the form shown in the parsing & plotting tutorial here (i.e. in the form output by the doped chemical potentials parsing functions), or you can also input these values manually as a dictionary alternatively.

If that doesn't work can you show what way you're inputting these values? Thanks!

[Atefeh-Yadegarifard]

Dear Sean, Hello. Thank you for your reply. As you said the problem was in the format of the chemlims. I was able to solve it by myself. Thanks again