diff --git a/docs/Dev_ToDo.md b/docs/Dev_ToDo.md index cd299b06..8eebcc00 100644 --- a/docs/Dev_ToDo.md +++ b/docs/Dev_ToDo.md @@ -19,7 +19,7 @@ ## Housekeeping - Tutorials general structure clean-up -- Clean up repo, removing old unnecessary git blobs +- Remnant TODOs in code - Docs: - Update note at end of thermo tutorial to link to py-sc-fermi/doped interface. @@ -86,6 +86,7 @@ ## SK To-Do for next update: - `doped` repo/docs cleanup `TODO`s above +- Clean up repo, removing old unnecessary git blobs - Note in chempots tutorial that LaTeX table generator website can also be used with the `to_csv()` function to generate LaTeX tables for the competing phases. - Quick-start tutorial suggested by Alex G - Add chempot grid plotting tool, shown in `JOSS_plots` using Alex's chemical potential grid, and test (and remove TODO from JOSS plots notebook). diff --git a/doped/VASP_sets/RelaxSet.yaml b/doped/VASP_sets/RelaxSet.yaml index 72752df7..6aaa1ddd 100644 --- a/doped/VASP_sets/RelaxSet.yaml +++ b/doped/VASP_sets/RelaxSet.yaml @@ -13,7 +13,7 @@ INCAR: LORBIT: 11 # lm-decomposed orbital projections; useful for DOS analysis (e.g. sumo-dosplot) LREAL: AUTO NCORE: 16 - NEDOS: 2000 # dense energy mesh output in OUTCAR; useful for DOS analysis (e.g. sumo-dosplot) + NEDOS: 3000 # dense energy mesh output in OUTCAR; useful for DOS analysis (e.g. sumo-dosplot) NSW: 100 PREC: Accurate SIGMA: 0.05 diff --git a/doped/analysis.py b/doped/analysis.py index bf524076..d6d37212 100644 --- a/doped/analysis.py +++ b/doped/analysis.py @@ -21,6 +21,7 @@ from pymatgen.analysis.structure_matcher import ElementComparator, StructureMatcher from pymatgen.core.sites import PeriodicSite from pymatgen.core.structure import Structure +from pymatgen.electronic_structure.dos import FermiDos from pymatgen.ext.matproj import MPRester from pymatgen.io.vasp.inputs import Poscar from pymatgen.io.vasp.outputs import Procar, Vasprun @@ -1275,6 +1276,8 @@ def get_defect_thermodynamics( band_gap: Optional[float] = None, dist_tol: float = 1.5, check_compatibility: bool = True, + bulk_dos: Optional[FermiDos] = None, + skip_check: bool = False, ) -> DefectThermodynamics: r""" Generates a DefectThermodynamics object from the parsed ``DefectEntry`` @@ -1342,6 +1345,26 @@ def get_defect_thermodynamics( Whether to check the compatibility of the bulk entry for each defect entry (i.e. that all reference bulk energies are the same). (Default: True) + bulk_dos (FermiDos or Vasprun or PathLike): + ``pymatgen`` ``FermiDos`` for the bulk electronic density of states (DOS), + for calculating Fermi level positions and defect/carrier concentrations. + Alternatively, can be a ``pymatgen`` ``Vasprun`` object or path to the + ``vasprun.xml(.gz)`` output of a bulk DOS calculation in VASP. + Can also be provided later when using ``get_equilibrium_fermi_level()``, + ``get_quenched_fermi_level_and_concentrations`` etc, or set using + ``DefectThermodynamics.bulk_dos = ...`` (with the same input options). + + Usually this is a static calculation with the `primitive` cell of the bulk + material, with relatively dense `k`-point sampling (especially for materials + with disperse band edges) to ensure an accurately-converged DOS and thus Fermi + level. ``ISMEAR = -5`` (tetrahedron smearing) is usually recommended for best + convergence wrt `k`-point sampling. Consistent functional settings should be + used for the bulk DOS and defect supercell calculations. + (Default: None) + skip_check (bool): + Whether to skip the warning about the DOS VBM differing from the defect + entries VBM by >0.05 eV. Should only be used when the reason for this + difference is known/acceptable. (Default: False) Returns: doped DefectThermodynamics object (``DefectThermodynamics``) @@ -1361,6 +1384,8 @@ def get_defect_thermodynamics( band_gap=band_gap, dist_tol=dist_tol, check_compatibility=check_compatibility, + bulk_dos=bulk_dos, + skip_check=skip_check, ) def __repr__(self): diff --git a/doped/thermodynamics.py b/doped/thermodynamics.py index 57203531..610cafbc 100644 --- a/doped/thermodynamics.py +++ b/doped/thermodynamics.py @@ -417,6 +417,8 @@ def __init__( band_gap: Optional[float] = None, dist_tol: float = 1.5, check_compatibility: bool = True, + bulk_dos: Optional[FermiDos] = None, + skip_check: bool = False, ): r""" Create a DefectThermodynamics object, which can be used to analyse the @@ -494,6 +496,26 @@ def __init__( Whether to check the compatibility of the bulk entry for each defect entry (i.e. that all reference bulk energies are the same). (Default: True) + bulk_dos (FermiDos or Vasprun or PathLike): + ``pymatgen`` ``FermiDos`` for the bulk electronic density of states (DOS), + for calculating Fermi level positions and defect/carrier concentrations. + Alternatively, can be a ``pymatgen`` ``Vasprun`` object or path to the + ``vasprun.xml(.gz)`` output of a bulk DOS calculation in VASP. + Can also be provided later when using ``get_equilibrium_fermi_level()``, + ``get_quenched_fermi_level_and_concentrations`` etc, or set using + ``DefectThermodynamics.bulk_dos = ...`` (with the same input options). + + Usually this is a static calculation with the `primitive` cell of the bulk + material, with relatively dense `k`-point sampling (especially for materials + with disperse band edges) to ensure an accurately-converged DOS and thus Fermi + level. ``ISMEAR = -5`` (tetrahedron smearing) is usually recommended for best + convergence wrt `k`-point sampling. Consistent functional settings should be + used for the bulk DOS and defect supercell calculations. + (Default: None) + skip_check (bool): + Whether to skip the warning about the DOS VBM differing from the defect + entries VBM by >0.05 eV. Should only be used when the reason for this + difference is known/acceptable. (Default: False) Key Attributes: defect_entries (list): @@ -521,6 +543,14 @@ def __init__( entry (i.e. that all reference bulk energies are the same). bulk_formula (str): The reduced formula of the bulk structure (e.g. "CdTe"). + bulk_dos (FermiDos): + ``pymatgen`` ``FermiDos`` for the bulk electronic density of states + (DOS), used for calculating Fermi level positions and defect/carrier + concentrations. + skip_check (bool): + Whether to skip the warning about the DOS VBM differing from the defect + entries VBM by >0.05 eV. Should only be used when the reason for this + difference is known/acceptable. """ if isinstance(defect_entries, dict): if not defect_entries: @@ -534,6 +564,7 @@ def __init__( self._chempots, self._el_refs = _parse_chempots(chempots, el_refs) self._dist_tol = dist_tol self.check_compatibility = check_compatibility + self.skip_check = skip_check # get and check VBM/bandgap values: def _raise_VBM_band_gap_value_error(vals, type="VBM"): @@ -574,6 +605,8 @@ def _raise_VBM_band_gap_value_error(vals, type="VBM"): f"(calculation_metadata attributes). Please specify the {name} in the function input." ) + self.bulk_dos = bulk_dos # use setter method, needs to be after setting VBM + # order entries for deterministic behaviour (particularly for plotting) self._sort_parse_and_check_entries() @@ -626,13 +659,17 @@ def as_dict(self): "vbm": self.vbm, "band_gap": self.band_gap, "dist_tol": self.dist_tol, + "check_compatibility": self.check_compatibility, + "bulk_formula": self.bulk_formula, + "bulk_dos": self.bulk_dos, + "skip_check": self.skip_check, } @classmethod def from_dict(cls, d): """ Reconstitute a DefectThermodynamics object from a dict representation - created using as_dict(). + created using ``as_dict()``. Args: d (dict): dict representation of DefectThermodynamics. @@ -644,6 +681,7 @@ def from_dict(cls, d): "ignore", "Use of properties is" ) # `message` only needs to match start of message defect_entries = [DefectEntry.from_dict(entry_dict) for entry_dict in d.get("defect_entries")] + bulk_dos = FermiDos.from_dict(d.get("bulk_dos")) if d.get("bulk_dos") else None return cls( defect_entries, @@ -652,6 +690,9 @@ def from_dict(cls, d): vbm=d.get("vbm"), band_gap=d.get("band_gap"), dist_tol=d.get("dist_tol"), + check_compatibility=d.get("check_compatibility"), + bulk_dos=bulk_dos, + skip_check=d.get("skip_check"), ) def to_json(self, filename: Optional[PathLike] = None): @@ -1181,6 +1222,39 @@ def el_refs(self, input_el_refs): """ self._chempots, self._el_refs = _parse_chempots(self._chempots, input_el_refs) + @property + def bulk_dos(self): + """ + Get the ``pymatgen`` ``FermiDos`` for the bulk electronic density of + states (DOS), for calculating Fermi level positions and defect/carrier + concentrations, if set. + + Otherwise, returns None. + """ + return self._bulk_dos + + @bulk_dos.setter + def bulk_dos(self, input_bulk_dos: Union[FermiDos, Vasprun, PathLike]): + r""" + Set the ``pymatgen`` ``FermiDos`` for the bulk electronic density of + states (DOS), for calculating Fermi level positions and defect/carrier + concentrations. + + Should be a ``pymatgen`` ``FermiDos`` for the bulk electronic DOS, a + ``pymatgen`` ``Vasprun`` object or path to the ``vasprun.xml(.gz)`` + output of a bulk DOS calculation in VASP. + Can also be provided later when using ``get_equilibrium_fermi_level()``, + ``get_quenched_fermi_level_and_concentrations`` etc. + + Usually this is a static calculation with the `primitive` cell of the bulk + material, with relatively dense `k`-point sampling (especially for materials + with disperse band edges) to ensure an accurately-converged DOS and thus Fermi + level. ``ISMEAR = -5`` (tetrahedron smearing) is usually recommended for best + convergence wrt `k`-point sampling. Consistent functional settings should be + used for the bulk DOS and defect supercell calculations. + """ + self._bulk_dos = self._parse_fermi_dos(input_bulk_dos, skip_check=self.skip_check) + @property def defect_names(self): """ @@ -1433,8 +1507,11 @@ def get_equilibrium_concentrations( ) def _parse_fermi_dos( - self, bulk_dos: Union[PathLike, Vasprun, FermiDos], skip_check: bool = False + self, bulk_dos: Optional[Union[PathLike, Vasprun, FermiDos]] = None, skip_check: bool = False ) -> FermiDos: + if bulk_dos is None: + return None + if isinstance(bulk_dos, FermiDos): fdos = bulk_dos # most similar settings to Vasprun.eigenvalue_band_properties: @@ -1526,7 +1603,7 @@ def get_equilibrium_fermi_level( used for the bulk DOS and defect supercell calculations. ``bulk_dos`` can also be left as ``None`` (default), if it has previously - been provided and parsed, and thus is set as the ``self.fermi_dos`` attribute. + been provided and parsed, and thus is set as the ``self.bulk_dos`` attribute. chempots (dict): Dictionary of chemical potentials to use for calculating the defect formation energies (and thus concentrations and Fermi level). @@ -1590,11 +1667,12 @@ def get_equilibrium_fermi_level( corresponding electron and hole concentrations (in cm^-3) if ``return_concs=True``. """ if bulk_dos is not None: - self.fermi_dos = self._parse_fermi_dos(bulk_dos, skip_check=skip_check) - elif not hasattr(self, "fermi_dos"): + self.bulk_dos = self._parse_fermi_dos(bulk_dos, skip_check=skip_check) + + if self.bulk_dos is None: # none provided, and none previously set raise ValueError( "No bulk DOS calculation (`bulk_dos`) provided or previously parsed to " - "`DefectThermodynamics.fermi_dos`, which is required for calculating carrier " + "`DefectThermodynamics.bulk_dos`, which is required for calculating carrier " "concentrations and solving for Fermi level position." ) @@ -1618,7 +1696,7 @@ def _get_total_q(fermi_level): conc_df = _add_effective_dopant_concentration(conc_df, effective_dopant_concentration) qd_tot = (conc_df["Charge"] * conc_df["Concentration (cm^-3)"]).sum() qd_tot += get_doping( - fermi_dos=self.fermi_dos, fermi_level=fermi_level + self.vbm, temperature=temperature + fermi_dos=self.bulk_dos, fermi_level=fermi_level + self.vbm, temperature=temperature ) return qd_tot @@ -1629,7 +1707,7 @@ def _get_total_q(fermi_level): eq_fermi_level: float = brentq(_get_total_q, -1.0, self.band_gap + 1.0) # type: ignore if return_concs: e_conc, h_conc = get_e_h_concs( - self.fermi_dos, eq_fermi_level + self.vbm, temperature # type: ignore + self.bulk_dos, eq_fermi_level + self.vbm, temperature # type: ignore ) return eq_fermi_level, e_conc, h_conc @@ -1718,7 +1796,7 @@ def get_quenched_fermi_level_and_concentrations( used for the bulk DOS and defect supercell calculations. ``bulk_dos`` can also be left as ``None`` (default), if it has previously - been provided and parsed, and thus is set as the ``self.fermi_dos`` attribute. + been provided and parsed, and thus is set as the ``self.bulk_dos`` attribute. chempots (dict): Dictionary of chemical potentials to use for calculating the defect formation energies (and thus concentrations and Fermi level). @@ -1818,14 +1896,15 @@ def get_quenched_fermi_level_and_concentrations( raise ValueError(f"Invalid keyword arguments: {', '.join(kwargs.keys())}") if bulk_dos is not None: - self.fermi_dos = self._parse_fermi_dos(bulk_dos, skip_check=kwargs.get("skip_check", False)) - elif not hasattr(self, "fermi_dos"): + self.bulk_dos = self._parse_fermi_dos(bulk_dos, skip_check=kwargs.get("skip_check", False)) + + if self.bulk_dos is None: # none provided, and none previously set raise ValueError( "No bulk DOS calculation (`bulk_dos`) provided or previously parsed to " - "`DefectThermodynamics.fermi_dos`, which is required for calculating carrier " + "`DefectThermodynamics.bulk_dos`, which is required for calculating carrier " "concentrations and solving for Fermi level position." ) - orig_fermi_dos = deepcopy(self.fermi_dos) # can get modified during annealing loops + orig_fermi_dos = deepcopy(self.bulk_dos) # can get modified during annealing loops chempots, el_refs = self._get_chempots( chempots, el_refs @@ -1835,11 +1914,11 @@ def get_quenched_fermi_level_and_concentrations( _no_chempots_warning() annealing_dos = ( - self.fermi_dos + self.bulk_dos if delta_gap == 0 else scissor_dos( delta_gap, - self.fermi_dos, + self.bulk_dos, verbose=kwargs.get("verbose", False), tol=kwargs.get("tol", 1e-8), ) @@ -1861,7 +1940,7 @@ def get_quenched_fermi_level_and_concentrations( # gap not 0 ) assert not isinstance(annealing_fermi_level, tuple) # float w/ return_concs=False, for typing - self.fermi_dos = orig_fermi_dos # reset to original DOS for quenched calculations + self.bulk_dos = orig_fermi_dos # reset to original DOS for quenched calculations annealing_defect_concentrations = self.get_equilibrium_concentrations( chempots=chempots, diff --git a/examples/CdTe/CdTe_LZ_thermo_wout_meta.json.gz b/examples/CdTe/CdTe_LZ_thermo_wout_meta.json.gz index 78dc28fe..7fc88142 100644 Binary files a/examples/CdTe/CdTe_LZ_thermo_wout_meta.json.gz and b/examples/CdTe/CdTe_LZ_thermo_wout_meta.json.gz differ diff --git a/tests/test_analysis.py b/tests/test_analysis.py index 062fe921..c451221f 100644 --- a/tests/test_analysis.py +++ b/tests/test_analysis.py @@ -16,6 +16,7 @@ import pytest from monty.serialization import dumpfn, loadfn from pymatgen.core.structure import Structure +from pymatgen.electronic_structure.dos import FermiDos from test_thermodynamics import custom_mpl_image_compare from doped.analysis import ( @@ -341,7 +342,11 @@ def test_DefectsParser_CdTe(self): ) # integration test using parsed CdTe thermo and chempots for plotting: - default_thermo = default_dp.get_defect_thermodynamics() + default_thermo = default_dp.get_defect_thermodynamics( + bulk_dos=os.path.join(self.CdTe_EXAMPLE_DIR, "CdTe_prim_k181818_NKRED_2_vasprun.xml.gz"), + ) # test providing bulk DOS + assert isinstance(default_thermo.bulk_dos, FermiDos) + assert np.isclose(default_thermo.bulk_dos.get_cbm_vbm()[1], 1.65, atol=1e-2) return default_thermo.plot(chempots=self.CdTe_chempots, limit="CdTe-Te") diff --git a/tests/test_thermodynamics.py b/tests/test_thermodynamics.py index 427cfdd4..6793edcd 100644 --- a/tests/test_thermodynamics.py +++ b/tests/test_thermodynamics.py @@ -19,6 +19,7 @@ import pytest from monty.serialization import dumpfn, loadfn from pymatgen.core.composition import Composition +from pymatgen.electronic_structure.dos import FermiDos from doped.generation import _sort_defect_entries from doped.thermodynamics import DefectThermodynamics, get_fermi_dos, scissor_dos @@ -178,6 +179,10 @@ def setUpClass(cls): cls.orig_ZnS_defect_thermo = loadfn(os.path.join(data_dir, "ZnS/ZnS_thermo.json.gz")) + # cls.CdTe_fermi_dos = get_fermi_dos( + # os.path.join(cls.CdTe_EXAMPLE_DIR, "CdTe_prim_k181818_NKRED_2_vasprun.xml.gz") + # ) # not used twice yet + class DefectThermodynamicsTestCase(DefectThermodynamicsSetupMixin): def _compare_defect_thermo_and_dict(self, defect_thermo, defect_dict): @@ -456,13 +461,34 @@ def test_initialisation_and_attributes(self): el_refs=self.V2O5_chempots["elemental_refs"], ) - print(f"Checking {name} with mismatching chempots") + print("Checking CdTe with mismatching chempots") with warnings.catch_warnings(record=True) as w: - _defect_thermo = DefectThermodynamics(self.CdTe_defect_dict, chempots={"Cd": -1.0, "Te": -6}) + defect_thermo = DefectThermodynamics(self.CdTe_defect_dict, chempots={"Cd": -1.0, "Te": -6}) print([str(warning.message) for warning in w]) # for debugging assert len(w) == 1 # only chempot incompatibility warning assert str(w[0].message) == self.cdte_chempot_warning_message + assert defect_thermo.bulk_dos is None + defect_thermo.bulk_dos = os.path.join( + self.CdTe_EXAMPLE_DIR, "CdTe_prim_k181818_NKRED_2_vasprun.xml.gz" + ) + assert isinstance(defect_thermo.bulk_dos, FermiDos) + assert np.isclose(defect_thermo.bulk_dos.get_cbm_vbm()[1], 1.65, atol=1e-2) + + defect_thermo.bulk_dos = defect_thermo.bulk_dos # test setting with FermiDos input + assert np.isclose(defect_thermo.bulk_dos.get_cbm_vbm()[1], 1.65, atol=1e-2) + + print("Checking CdTe with fermi dos") + with warnings.catch_warnings(record=True) as w: + defect_thermo = DefectThermodynamics( + self.CdTe_defect_dict, + bulk_dos=os.path.join(self.CdTe_EXAMPLE_DIR, "CdTe_prim_k181818_NKRED_2_vasprun.xml.gz"), + ) + print([str(warning.message) for warning in w]) # for debugging + assert len(w) == 0 # no warning (e.g. VBM mismatch) + assert isinstance(defect_thermo.bulk_dos, FermiDos) + assert np.isclose(defect_thermo.bulk_dos.get_cbm_vbm()[1], 1.65, atol=1e-2) + def test_DefectsParser_thermo_objs(self): """ Test the `DefectThermodynamics` objects created from the