Merge pull request #80 from SMTG-Bham/develop

Develop

CHANGELOG.rst

@@ -1,6 +1,14 @@
 Change Log
 ==========
 
+v.2.4.6
+----------
+- Update ``Defect``, ``DefectEntry`` and ``DefectThermodynamics`` properties/methods to be even more
+  efficient with calculations of formation energies and concentrations. Gives ~10x speedup in Fermi
+  solving and concentration calculations (e.g. from 2 hours to 12 minutes for 2D chempot vs temp CdTe grid
+  in thermodynamics tutorial).
+- Avoid unnecessary ``DeprecationWarning``s from latest ``spglib`` release.
+
 v.2.4.5
 ----------
 - Enforce ``shakenbreak>=2.3.4`` requirement.

README.md

@@ -62,8 +62,12 @@ As shown in the `doped` tutorials, it is highly recommended to use the [`ShakeNB
 
 - B. E. Murdock et al. **_Li-Site Defects Induce Formation of Li-Rich Impurity Phases: Implications for Charge Distribution and Performance of LiNi<sub>0.5-x</sub>M<sub>x</sub>Mn<sub>1.5</sub>O<sub>4</sub> Cathodes (M = Fe and Mg; x = 0.05–0.2)_** [_Advanced Materials_](https://doi.org/10.1002/adma.202400343) 2024
 - A. G. Squires et al. **_Oxygen dimerization as a defect-driven process in bulk LiNiO2<sub>2</sub>_** [_ChemRxiv_](https://doi.org/10.26434/chemrxiv-2024-lcmkj) 2024
+- Y. Fu & H. Lohan et al. **_Factors Enabling Delocalized Charge-Carriers in Pnictogen-Based
+Solar Absorbers: In-depth Investigation into CuSbSe<sub>2</sub>_** [_arXiv_](https://doi.org/10.48550/arXiv.2401.02257) 2024
+- S. Hachmioune et al. **_Exploring the Thermoelectric Potential of MgB4: Electronic Band Structure, Transport Properties, and Defect Chemistry_** [_Chemistry of Materials_](https://doi.org/10.1021/acs.chemmater.4c00584) 2024
+- J. Hu et al. **_Enabling ionic transport in Li3AlP2 the roles of defects and disorder_** [_ChemRxiv_](https://doi.org/10.26434/chemrxiv-2024-3s0kh) 2024
 - X. Wang et al. **_Upper efficiency limit of Sb<sub>2</sub>Se<sub>3</sub> solar cells_** [_Joule_](https://doi.org/10.1016/j.joule.2024.05.004) 2024
-- I. Mosquera-Lois et al. **_Machine-learning structural reconstructions for accelerated point defect calculations_** [_arXiv_](https://doi.org/10.48550/arXiv.2401.12127) 2024
+- I. Mosquera-Lois et al. **_Machine-learning structural reconstructions for accelerated point defect calculations_** [_npj Computational Materials_](https://doi.org/10.1038/s41524-024-01303-9) 2024
 - W. Dou et al. **_Giant Band Degeneracy via Orbital Engineering Enhances Thermoelectric Performance from Sb<sub>2</sub>Si<sub>2</sub>Te<sub>6</sub> to Sc<sub>2</sub>Si<sub>2</sub>Te<sub>6</sub>_** [_ChemRxiv_](https://doi.org/10.26434/chemrxiv-2024-hm6vh) 2024
 - K. Li et al. **_Computational Prediction of an Antimony-based n-type Transparent Conducting Oxide: F-doped Sb<sub>2</sub>O<sub>5</sub>_** [_Chemistry of Materials_](https://doi.org/10.1021/acs.chemmater.3c03257) 2023
 - X. Wang et al. **_Four-electron negative-U vacancy defects in antimony selenide_** [_Physical Review B_](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.108.134102) 2023

docs/conf.py

@@ -25,7 +25,7 @@
 author = 'Seán R. Kavanagh'
 
 # The full version, including alpha/beta/rc tags
-release = '2.4.5'
+release = '2.4.6'
 
 
 # -- General configuration ---------------------------------------------------

docs/index.rst

@@ -91,8 +91,11 @@ Studies using ``doped``, so far
 
 - B\. E. Murdock et al. **Li-Site Defects Induce Formation of Li-Rich Impurity Phases: Implications for Charge Distribution and Performance of LiNi** :sub:`0.5-x` **M** :sub:`x` **Mn** :sub:`1.5` **O** :sub:`4`  **Cathodes (M = Fe and Mg; x = 0.05–0.2)** `Advanced Materials <https://doi.org/10.1002/adma.202400343>`_ 2024
 - A\. G. Squires et al. **Oxygen dimerization as a defect-driven process in bulk LiNiO₂** `ChemRxiv <https://doi.org/10.26434/chemrxiv-2024-lcmkj>`_ 2024
+- Y\. Fu & H. Lohan et al. **Factors Enabling Delocalized Charge-Carriers in Pnictogen-Based Solar Absorbers: In-depth Investigation into CuSbSe<sub>2</sub>** `arXiv <https://doi.org/10.48550/arXiv.2401.02257>`_ 2024
+- S\. Hachmioune et al. **Exploring the Thermoelectric Potential of MgB4: Electronic Band Structure, Transport Properties, and Defect Chemistry** `Chemistry of Materials <https://doi.org/10.1021/acs.chemmater.4c00584>`_ 2024
+- J\. Hu et al. **Enabling ionic transport in Li3AlP2 the roles of defects and disorder** `ChemRxiv <https://doi.org/10.26434/chemrxiv-2024-3s0kh>`_ 2024
 - X\. Wang et al. **Upper efficiency limit of Sb₂Se₃ solar cells** `Joule <https://doi.org/10.1016/j.joule.2024.05.004>`_ 2024
-- I\. Mosquera-Lois et al. **Machine-learning structural reconstructions for accelerated point defect calculations** `arXiv <https://doi.org/10.48550/arXiv.2401.12127>`_ 2024
+- I\. Mosquera-Lois et al. **Machine-learning structural reconstructions for accelerated point defect calculations** `npj Computational Materials <https://doi.org/10.1038/s41524-024-01303-9>`_ 2024
 - W\. Dou et al. **Giant Band Degeneracy via Orbital Engineering Enhances Thermoelectric Performance from Sb₂Si₂Te₆ to Sc₂Si₂Te₆** `ChemRxiv <https://doi.org/10.26434/chemrxiv-2024-hm6vh>`_ 2024
 - K\. Li et al. **Computational Prediction of an Antimony-based n-type Transparent Conducting Oxide: F-doped Sb₂O₅** `Chemistry of Materials <https://doi.org/10.1021/acs.chemmater.3c03257>`_ 2024
 - X\. Wang et al. **Four-electron negative-U vacancy defects in antimony selenide** `Physical Review B <https://journals.aps.org/prb/abstract/10.1103/PhysRevB.108.134102>`_ 2023

doped/__init__.py

@@ -64,17 +64,16 @@ def _ignore_pmg_warnings():
     warnings.filterwarnings("ignore", message="Use get_magnetic_symmetry()")
     warnings.filterwarnings("ignore", message="Use of properties is now deprecated")
 
-    warnings.filterwarnings("ignore", message="get_vasp_input")  # deprecation warning introduced
-    # in pymatgen>2024.1.6, fixed in our PR: https://github.com/materialsproject/pymatgen/pull/3601
-    # (now merged) -- delete later if pymatgen requirement is updated beyond this
-
     # avoid warning about selective_dynamics properties (can happen if user explicitly set "T T T" (or
     # otherwise) for the bulk):
     warnings.filterwarnings("ignore", message="Not all sites have property")
 
     # ignore warning about structure charge that appears when getting Vasprun.as_dict():
     warnings.filterwarnings("ignore", message="Structure charge")
 
+    # SpglibDataset warning introduced in v2.4.1
+    warnings.filterwarnings("ignore", message="dict interface")
+
 
 _ignore_pmg_warnings()
 

doped/analysis.py

@@ -217,6 +217,7 @@ def defect_from_structures(
             Dictionary of bulk supercell Voronoi node information, for
             further expedited site-matching.
     """
+    warnings.filterwarnings("ignore", "dict interface")  # ignore spglib warning from v2.4.1
     try:
         def_type, comp_diff = get_defect_type_and_composition_diff(bulk_supercell, defect_supercell)
     except RuntimeError as exc:

doped/core.py

@@ -29,6 +29,10 @@
 
     from easyunfold.procar import Procar as EasyunfoldProcar
 
+# SpglibDataset warning introduced in v2.4.1, can be ignored for now
+warnings.filterwarnings("ignore", message="dict interface")
+
+
 _orientational_degeneracy_warning = (
     "The defect supercell has been detected to possibly have a non-scalar matrix expansion, "
     "which could be breaking the cell periodicity and possibly preventing the correct _relaxed_ "
@@ -148,6 +152,10 @@ class DefectEntry(thermo.DefectEntry):
     # codes etc)
     equivalent_supercell_sites: list[PeriodicSite] = field(default_factory=list)
     bulk_supercell: Optional[Structure] = None
+    _bulk_entry_energy: Optional[float] = None
+    _bulk_entry_hash: Optional[int] = None
+    _sc_entry_energy: Optional[float] = None
+    _sc_entry_hash: Optional[int] = None
 
     def __post_init__(self):
         """
@@ -167,32 +175,6 @@ def __post_init__(self):
                 f"{name_wout_charge}_{'+' if self.charge_state > 0 else ''}{self.charge_state}"
             )
 
-    def _check_if_multiple_finite_size_corrections(self):
-        """
-        Checks that there is no double counting of finite-size charge
-        corrections, in the defect_entry.corrections dict.
-        """
-        matching_finite_size_correction_keys = {
-            key
-            for key in self.corrections
-            if any(x in key for x in ["FNV", "freysoldt", "Freysoldt", "Kumagai", "kumagai"])
-        }
-        if len(matching_finite_size_correction_keys) > 1:
-            warnings.warn(
-                f"It appears there are multiple finite-size charge corrections in the corrections dict "
-                f"attribute for defect {self.name}. These are:"
-                f"\n{matching_finite_size_correction_keys}."
-                f"\nPlease ensure there is no double counting / duplication of energy corrections!"
-            )
-
-    @property
-    def corrected_energy(self) -> float:
-        """
-        The energy of the defect entry with `all` corrections applied.
-        """
-        self._check_if_multiple_finite_size_corrections()
-        return self.sc_entry.energy + sum(self.corrections.values())
-
     def to_json(self, filename: Optional[str] = None):
         """
         Save the ``DefectEntry`` object to a json file, which can be reloaded
@@ -829,13 +811,61 @@ def get_eigenvalue_analysis(
 
     def _get_chempot_term(self, chemical_potentials=None):
         chemical_potentials = chemical_potentials or {}
+        element_changes = {elt.symbol: change for elt, change in self.defect.element_changes.items()}
 
         return sum(
-            chem_pot * -self.defect.element_changes[Element(el)]
+            chem_pot * -element_changes[el]
             for el, chem_pot in chemical_potentials.items()
-            if Element(el) in self.defect.element_changes
+            if el in element_changes
         )
 
+    def get_ediff(self) -> float:
+        """
+        Get the energy difference between the defect and bulk supercell
+        calculations, including finite-size corrections.
+
+        Refactored from ``pymatgen-analysis-defects`` to be more efficient,
+        reducing compute times when looping over formation energy /
+        concentration functions.
+        """
+        if self.bulk_entry is None:
+            raise RuntimeError(
+                "Attempting to compute the energy difference without a defined bulk entry for the "
+                "DefectEntry object!"
+            )
+        return self.corrected_energy - self.bulk_entry_energy
+
+    @property
+    def corrected_energy(self) -> float:
+        """
+        Get the energy of the defect supercell calculation, with `all`
+        corrections (in ``DefectEntry.corrections``) applied.
+
+        Refactored from ``pymatgen-analysis-defects`` to be more efficient,
+        reducing compute times when looping over formation energy /
+        concentration functions.
+        """
+        self._check_if_multiple_finite_size_corrections()
+        return self.sc_entry_energy + sum(self.corrections.values())
+
+    def _check_if_multiple_finite_size_corrections(self):
+        """
+        Checks that there is no double counting of finite-size charge
+        corrections, in the defect_entry.corrections dict.
+        """
+        matching_finite_size_correction_keys = {
+            key
+            for key in self.corrections
+            if any(x in key for x in ["FNV", "freysoldt", "Freysoldt", "Kumagai", "kumagai"])
+        }
+        if len(matching_finite_size_correction_keys) > 1:
+            warnings.warn(
+                f"It appears there are multiple finite-size charge corrections in the corrections dict "
+                f"attribute for defect {self.name}. These are:"
+                f"\n{matching_finite_size_correction_keys}."
+                f"\nPlease ensure there is no double counting / duplication of energy corrections!"
+            )
+
     def formation_energy(
         self,
         chempots: Optional[dict] = None,
@@ -1011,6 +1041,7 @@ def equilibrium_concentration(
         vbm: Optional[float] = None,
         per_site: bool = False,
         symprec: Optional[float] = None,
+        formation_energy: Optional[float] = None,
     ) -> float:
         r"""
         Compute the `equilibrium` concentration (in cm^-3) for the
@@ -1098,6 +1129,13 @@ def equilibrium_concentration(
                 If ``symprec`` is set, then the point symmetries and corresponding
                 orientational degeneracy will be re-parsed/computed even if already
                 present in the ``DefectEntry`` object ``calculation_metadata``.
+            formation_energy (float):
+                Pre-calculated formation energy to use for the defect concentration
+                calculation, in order to reduce compute times (e.g. when looping over
+                many chemical potential / temperature / etc ranges). Only really intended
+                for internal ``doped`` usage. If ``None`` (default), will calculate the
+                formation energy using the input ``chempots``, ``limit``, ``el_refs``,
+                ``vbm`` and ``fermi_level`` arguments. (Default: None)
 
         Returns:
             Concentration in cm^-3 (or as fractional per site, if per_site = True) (float)
@@ -1135,22 +1173,24 @@ def equilibrium_concentration(
         if self.calculation_metadata.get("periodicity_breaking_supercell", False):
             warnings.warn(_orientational_degeneracy_warning)
 
-        formation_energy = self.formation_energy(  # if chempots is None, this will throw warning
-            chempots=chempots, limit=limit, el_refs=el_refs, vbm=vbm, fermi_level=fermi_level
-        )
+        if formation_energy is None:
+            formation_energy = self.formation_energy(  # if chempots is None, this will throw warning
+                chempots=chempots, limit=limit, el_refs=el_refs, vbm=vbm, fermi_level=fermi_level
+            )
 
         with np.errstate(over="ignore"):
             exp_factor = np.exp(
                 -formation_energy / (constants_value("Boltzmann constant in eV/K") * temperature)
             )
 
-        degeneracy_factor = (
-            reduce(lambda x, y: x * y, self.degeneracy_factors.values()) if self.degeneracy_factors else 1
-        )
-        if per_site:
-            return exp_factor * degeneracy_factor
+            degeneracy_factor = (
+                reduce(lambda x, y: x * y, self.degeneracy_factors.values())
+                if self.degeneracy_factors
+                else 1
+            )
+            if per_site:
+                return exp_factor * degeneracy_factor
 
-        with np.errstate(over="ignore"):
             return self.bulk_site_concentration * degeneracy_factor * exp_factor
 
     @property
@@ -1161,7 +1201,7 @@ def bulk_site_concentration(self):
         V_O in SrTiO3, returns the site concentration of (symmetry-equivalent)
         oxygen atoms in SrTiO3).
         """
-        volume_in_cm3 = self.defect.structure.volume * 1e-24  # convert volume in Å^3 to cm^3
+        volume_in_cm3 = self.defect.volume * 1e-24  # convert volume in Å^3 to cm^3
         return self.defect.multiplicity / volume_in_cm3
 
     def __repr__(self):
@@ -1187,9 +1227,54 @@ def __repr__(self):
     def __eq__(self, other):
         """
         Determine whether two ``DefectEntry`` objects are equal, by comparing
-        self.name and self.sc_entry.
+        ``self.name``, ``self.sc_entry``, ``self.bulk_entry`` and
+        ``self.corrections`` (i.e. name and energy match).
+        """
+        return (
+            self.name == other.name
+            and self.sc_entry == other.sc_entry
+            and self.bulk_entry == other.bulk_entry
+            and self.corrections == other.corrections
+        )
+
+    @property
+    def bulk_entry_energy(self):
+        r"""
+        Get the raw energy of the bulk supercell calculation (i.e.
+        ``bulk_entry.energy``).
+
+        Refactored from ``pymatgen-analysis-defects`` to be more efficient,
+        reducing compute times when looping over formation energy /
+        concentration functions.
+        """
+        if self.bulk_entry is None:
+            return None
+
+        if hasattr(self, "_bulk_entry_energy") and self._bulk_entry_hash == hash(self.bulk_entry):
+            return self._bulk_entry_energy
+
+        self._bulk_entry_hash = hash(self.bulk_entry)
+        self._bulk_entry_energy = self.bulk_entry.energy
+
+        return self._bulk_entry_energy
+
+    @property
+    def sc_entry_energy(self):
+        r"""
+        Get the raw energy of the defect supercell calculation (i.e.
+        ``sc_entry.energy``).
+
+        Refactored from ``pymatgen-analysis-defects`` to be more efficient,
+        reducing compute times when looping over formation energy /
+        concentration functions.
         """
-        return self.name == other.name and self.sc_entry == other.sc_entry
+        if hasattr(self, "_sc_entry_energy") and self._sc_entry_hash == hash(self.sc_entry):
+            return self._sc_entry_energy
+
+        self._sc_entry_hash = hash(self.sc_entry)
+        self._sc_entry_energy = self.sc_entry.energy
+
+        return self._sc_entry_energy
 
     def plot_site_displacements(
         self,
@@ -1283,7 +1368,10 @@ def _guess_and_set_struct_oxi_states(structure):
         if oxidation states could not be guessed.
     """
     bv_analyzer = BVAnalyzer()
-    with contextlib.suppress(ValueError):  # ValueError raised if oxi states can't be assigned
+    with contextlib.suppress(ValueError), warnings.catch_warnings():
+        # ValueError raised if oxi states can't be assigned, and SpglibDataset warning introduced in
+        # v2.4.1, can be ignored for now:
+        warnings.filterwarnings("ignore", message="dict interface")
         oxi_dec_structure = bv_analyzer.get_oxi_state_decorated_structure(structure)
         if all(
             np.isclose(int(specie.oxi_state), specie.oxi_state) for specie in oxi_dec_structure.species
@@ -1605,6 +1693,14 @@ def _set_oxi_state(self):
                 self.structure, timeout_1=5, timeout_2=5, break_early_if_expensive=True
             ):
                 self.structure = struct_w_oxi
+                if self.defect_type != core.DefectType.Interstitial:
+                    self._defect_site = min(
+                        self.structure.get_sites_in_sphere(
+                            self.site.coords,
+                            0.5,
+                        ),
+                        key=lambda x: x[1],
+                    )
             else:
                 self.oxi_state = "Undetermined"
                 return
@@ -1900,6 +1996,46 @@ def get_charge_states(self, padding: int = 1) -> list[int]:
 
         return charges
 
+    @property
+    def defect_site(self) -> PeriodicSite:
+        """
+        The defect site in the structure.
+
+        Re-written from ``pymatgen-analysis-defects`` version to
+        be far more efficient, when used in loops (e.g. for calculating
+        defect concentrations as functions of chemical potentials,
+        temperature etc.).
+        """
+        if self.defect_type == core.DefectType.Interstitial:
+            return self.site  # same as self.defect_site
+
+        # else defect_site is the closest site in ``structure`` to the provided ``site``:
+        if not hasattr(self, "_defect_site"):
+            self._defect_site = min(
+                self.structure.get_sites_in_sphere(
+                    self.site.coords,
+                    0.5,
+                ),
+                key=lambda x: x[1],
+            )
+
+        return self._defect_site
+
+    @property
+    def volume(self) -> float:
+        """
+        The volume (in ų) of the structure in which the defect is created
+        (i.e. ``Defect.structure``).
+
+        Ensures volume is only computed once when calculating defect
+        concentrations in loops (e.g. for calculating defect concentrations as
+        functions of chemical potentials, temperature etc.).
+        """
+        if not hasattr(self, "_volume"):
+            self._volume = self.structure.volume
+
+        return self._volume
+
 
 def doped_defect_from_pmg_defect(defect: core.Defect, bulk_oxi_states=False, **doped_kwargs):
     """