Add checks for chempot compatibility with defect calculations in `Def…

…ectThermodynamics`, and add tests

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

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): #