From 773d1d645167b13f169c6861260d2b459d252905 Mon Sep 17 00:00:00 2001 From: Sean Kavanagh Date: Tue, 30 Jul 2024 11:54:06 -0400 Subject: [PATCH] Add checks for chempot compatibility with defect calculations in `DefectThermodynamics`, and add tests --- doped/thermodynamics.py | 113 +++++++++++++++++++++++++++++++---- tests/test_thermodynamics.py | 54 ++++++++++++++++- 2 files changed, 156 insertions(+), 11 deletions(-) diff --git a/doped/thermodynamics.py b/doped/thermodynamics.py index c8b817cc..0892dd63 100644 --- a/doped/thermodynamics.py +++ b/doped/thermodynamics.py @@ -191,6 +191,41 @@ def _parse_chempots(chempots: Optional[dict] = None, el_refs: Optional[dict] = N return chempots, chempots.get("elemental_refs") +def raw_energy_from_chempots(composition: Union[str, dict, Composition], chempots: dict) -> float: + """ + Given an input composition (as a ``str``, ``dict`` or ``pymatgen`` + ``Composition`` object) and chemical potentials dictionary, get the + corresponding raw energy of the composition (i.e. taking the energies given + in the ``'limits'`` subdicts of ``chempots``, in the ``doped`` chemical + potentials dictionary format). + + Args: + composition (Union[str, dict, Composition]): + Composition to get the raw energy of. + chempots (dict): + Chemical potentials dictionary. + + Returns: + Raw energy of the composition. + """ + if not isinstance(composition, Composition): + composition = Composition(composition) + + if "limits" not in chempots: + chempots, _el_refs = _parse_chempots(chempots) + + raw_energies_dict = dict(next(iter(chempots["limits"].values()))) + + if any(el.symbol not in raw_energies_dict for el in composition.elements): + raise ValueError( + f"The chemical potentials dictionary (with elements {list(raw_energies_dict.keys())} does not " + f"contain all the elements in the host composition " + f"({[el.symbol for el in composition.elements]})!" + ) + + return sum(raw_energies_dict.get(el.symbol, 0) * stoich for el, stoich in composition.items()) + + def group_defects_by_distance( entry_list: list[DefectEntry], dist_tol: float = 1.5 ) -> dict[str, dict[tuple, list[DefectEntry]]]: @@ -540,7 +575,7 @@ def _raise_VBM_band_gap_value_error(vals, type="VBM"): ) # order entries for deterministic behaviour (particularly for plotting) - self._sort_parse_and_check_entries(check_compatibility=check_compatibility) + self._sort_parse_and_check_entries() bulk_entry = self.defect_entries[0].bulk_entry if bulk_entry is not None: @@ -550,10 +585,11 @@ def _raise_VBM_band_gap_value_error(vals, type="VBM"): else: self.bulk_formula = None - def _sort_parse_and_check_entries(self, check_compatibility: bool = True): + def _sort_parse_and_check_entries(self): """ Sort the defect entries, parse the transition levels, and check the - compatibility of the bulk entries (if check_compatibility is True). + compatibility of the bulk entries (if ``self.check_compatibility`` is + ``True``). """ defect_entries_dict: dict[str, DefectEntry] = {} for entry in self.defect_entries: # rename defect entry names in dict if necessary ("_a", "_b"...) @@ -572,8 +608,9 @@ def _sort_parse_and_check_entries(self, check_compatibility: bool = True): with warnings.catch_warnings(): # ignore formation energies chempots warning when just parsing TLs warnings.filterwarnings("ignore", message="No chemical potentials") self._parse_transition_levels() - if check_compatibility: + if self.check_compatibility: self._check_bulk_compatibility() + self._check_bulk_chempots_compatibility(self._chempots) def as_dict(self): """ @@ -661,7 +698,11 @@ def _get_chempots(self, chempots: Optional[dict] = None, el_refs: Optional[dict] Parse chemical potentials, either using input values (after formatting them in the doped format) or using the class attributes if set. """ - return _parse_chempots(chempots or self.chempots, el_refs or self.el_refs) + chempots, el_refs = _parse_chempots(chempots or self.chempots, el_refs or self.el_refs) + if self.check_compatibility: + self._check_bulk_chempots_compatibility(chempots) + + return chempots, el_refs def _parse_transition_levels(self): r""" @@ -890,7 +931,7 @@ def _check_bulk_compatibility(self): """ Helper function to quickly check if all entries have compatible bulk calculation settings, by checking that the energy of - defect_entry.bulk_entry is the same for all defect entries. + ``defect_entry.bulk_entry`` is the same for all defect entries. By proxy checks that same bulk/defect calculation settings were used in all cases, from each bulk/defect combination already being checked when @@ -907,8 +948,8 @@ def _check_bulk_compatibility(self): f"eV. This can lead to inaccuracies in predicted formation energies! The bulk energies of " f"defect entries in `defect_entries` are:\n" f"{[(entry.name, entry.bulk_entry.energy) for entry in self.defect_entries]}\n" - f"You can suppress this warning by setting `check_compatibility=False` in " - "`DefectThermodynamics` initialisation." + f"You can suppress this warning by setting `DefectThermodynamics.check_compatibility = " + f"False`." ) def _check_bulk_defects_compatibility(self): @@ -918,7 +959,7 @@ def _check_bulk_defects_compatibility(self): Currently not used, as the bulk/defect compatibility is checked when parsing, and the compatibility across bulk calculations is checked with - _check_bulk_compatibility(). + ``_check_bulk_compatibility()``. """ # check each defect entry against its own bulk, and also check each bulk against each other reference_defect_entry = self.defect_entries[0] @@ -964,6 +1005,55 @@ def _check_bulk_defects_compatibility(self): f"{defect_entry.name}: \n{concatenated_warnings}" ) + def _check_bulk_chempots_compatibility(self, chempots: Optional[dict] = None): + r""" + Helper function to quickly check if the supplied chemical potentials + dictionary matches the bulk supercell used for the defect calculations, + by comparing the raw energies (from the bulk supercell calculation, and + that corresponding to the chemical potentials supplied). + + Args: + chempots (dict, optional): + Dictionary of chemical potentials to check compatibility with + the bulk supercell calculations (``DefectEntry.bulk_entry``\s), + in the ``doped`` format. + + If ``None`` (default), will use ``self.chempots`` (= 0 for all + chemical potentials by default). + This can have the form of ``{"limits": [{'limit': [chempot_dict]}]}`` + (the format generated by ``doped``\'s chemical potential parsing + functions (see tutorials)), or alternatively a dictionary of chemical + potentials for a single limit (``limit``), in the format: + ``{element symbol: chemical potential}``. + If manually specifying chemical potentials this way, you can set the + ``el_refs`` option with the DFT reference energies of the elemental phases, + in which case it is the formal chemical potentials (i.e. relative to the + elemental references) that should be given here, otherwise the absolute + (DFT) chemical potentials should be given. + """ + if chempots is None and self.chempots is None: + return + + bulk_entry = next(entry.bulk_entry for entry in self.defect_entries) + bulk_supercell_energy_per_atom = bulk_entry.energy / bulk_entry.composition.num_atoms + bulk_chempot_energy_per_atom = ( + raw_energy_from_chempots(bulk_entry.composition, chempots or self.chempots) + / bulk_entry.composition.num_atoms + ) + + if abs(bulk_supercell_energy_per_atom - bulk_chempot_energy_per_atom) > 0.025: + warnings.warn( # 0.05 eV intrinsic defect formation energy error tolerance, taking per-atom + # chempot error and multiplying by 2 to account for how this would affect antisite + # formation energies (extreme case) + f"Note that the raw (DFT) energy of the bulk supercell calculation (" + f"{bulk_supercell_energy_per_atom:.2f} eV/atom) differs from that expected from the " + f"supplied chemical potentials ({bulk_chempot_energy_per_atom:.2f} eV/atom) by >0.025 eV. " + f"This will likely give inaccuracies of similar magnitude in the predicted formation " + f"energies! \n" + f"You can suppress this warning by setting `DefectThermodynamics.check_compatibility = " + f"False`." + ) + def add_entries( self, defect_entries: Union[list[DefectEntry], dict[str, DefectEntry]], @@ -986,6 +1076,7 @@ def add_entries( entry (i.e. that all reference bulk energies are the same). (Default: True) """ + self.check_compatibility = check_compatibility if isinstance(defect_entries, dict): defect_entries = list(defect_entries.values()) @@ -996,7 +1087,7 @@ def add_entries( ) self._defect_entries += defect_entries - self._sort_parse_and_check_entries(check_compatibility=check_compatibility) + self._sort_parse_and_check_entries() @property def defect_entries(self): @@ -1061,6 +1152,8 @@ def chempots(self, input_chempots): (Default: None) """ self._chempots, self._el_refs = _parse_chempots(input_chempots, self._el_refs) + if self.check_compatibility: + self._check_bulk_chempots_compatibility(self._chempots) @property def el_refs(self): diff --git a/tests/test_thermodynamics.py b/tests/test_thermodynamics.py index f578a1fd..427cfdd4 100644 --- a/tests/test_thermodynamics.py +++ b/tests/test_thermodynamics.py @@ -18,6 +18,7 @@ import pandas as pd import pytest from monty.serialization import dumpfn, loadfn +from pymatgen.core.composition import Composition from doped.generation import _sort_defect_entries from doped.thermodynamics import DefectThermodynamics, get_fermi_dos, scissor_dos @@ -102,6 +103,14 @@ def setUp(self): self.Sb2O5_defect_thermo = deepcopy(self.orig_Sb2O5_defect_thermo) self.ZnS_defect_thermo = deepcopy(self.orig_ZnS_defect_thermo) + self.cdte_chempot_warning_message = ( + "Note that the raw (DFT) energy of the bulk supercell calculation (-3.37 eV/atom) differs " + "from that expected from the supplied chemical potentials (-3.50 eV/atom) by >0.025 eV. This " + "will likely give inaccuracies of similar magnitude in the predicted formation energies! " + "\nYou can suppress this warning by setting `DefectThermodynamics.check_compatibility = " + "False`." + ) + @classmethod def setUpClass(cls): cls.module_path = os.path.dirname(os.path.abspath(__file__)) @@ -376,6 +385,17 @@ def _check_defect_thermo( self._set_and_check_dist_tol(1.0, defect_thermo, 4) self._set_and_check_dist_tol(0.5, defect_thermo, 5) + # test mismatching chempot warnings: + print("Checking mismatching chempots") + mismatch_chempots = {el: -3 for el in Composition(defect_thermo.bulk_formula).as_dict()} + with warnings.catch_warnings(record=True) as w: + defect_thermo.chempots = mismatch_chempots + print([str(warning.message) for warning in w]) # for debugging + assert any( + "Note that the raw (DFT) energy of the bulk supercell calculation" in str(warning.message) + for warning in w + ) + def _check_CdTe_example_dist_tol(self, defect_thermo, num_grouped_defects): print(f"Testing CdTe updated dist_tol: {defect_thermo.dist_tol}") tl_df = defect_thermo.get_transition_levels() @@ -436,6 +456,13 @@ def test_initialisation_and_attributes(self): el_refs=self.V2O5_chempots["elemental_refs"], ) + print(f"Checking {name} with mismatching chempots") + with warnings.catch_warnings(record=True) as w: + _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 + def test_DefectsParser_thermo_objs(self): """ Test the `DefectThermodynamics` objects created from the @@ -1369,7 +1396,7 @@ def _check_formation_energy_methods(form_en_df_row, thermo_obj, fermi_level): el_refs={"Cd": -12, "Te": -1}, ) assert "Fermi level was not set" not in output - assert not w + assert len(w) == 1 # only mis-matching chempot warning assert manual_form_en_df.shape == (7, 10) # test sum of formation energy terms equals total and other formation energy df properties: self._check_form_en_df(manual_form_en_df, fermi_level=3, defect_thermo=self.CdTe_defect_thermo) @@ -1939,6 +1966,31 @@ def test_symmetry_degeneracy_unparsed(self): (e.g. transferring from old doped versions, from pymatgen-analysis- defects objects etc). """ + # TODO + + def test_incompatible_chempots_warning(self): + """ + Test that we get the expected warnings when we provide incompatible + chemical potentials for our DefectThermodynamics object. + """ + slightly_off_chempots = {"Cd": -1.0, "Te": -6} + for func in [ + self.CdTe_defect_thermo.get_equilibrium_concentrations, + self.CdTe_defect_thermo.get_dopability_limits, + self.CdTe_defect_thermo.get_doping_windows, + self.CdTe_defect_thermo.get_formation_energies, + self.CdTe_defect_thermo.plot, + ]: + _result, _tl_output, w = _run_func_and_capture_stdout_warnings( + func, chempots=slightly_off_chempots + ) + assert any(str(warning.message) == self.cdte_chempot_warning_message for warning in w) + + with warnings.catch_warnings(record=True) as w: + self.CdTe_defect_thermo.chempots = slightly_off_chempots + 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 def belas_linear_fit(T): #