Allow bulk_dos to be set in DefectThermodynamics init, or as prop…

…erty (with setter method), and test – quicker init times for `FermiSolver`

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

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

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

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,

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