Update FermiDos.get_doping() to be more robust (#3879)

* Update `FermiDos.get_doping()` to be more robust.

* Update `tol` defaults for DOS methods, and add small notes to `Vasprun` DOS parsing docstrings

* pre-commit auto-fixes

* Use `scipy.special.expit` function for Fermi-Dirac distribution, slightly faster, cleaner and no masking required to avoid overflow warnings

* Update `FermiDos.get_doping()` to be more robust.

* Update `tol` defaults for DOS methods, and add small notes to `Vasprun` DOS parsing docstrings

* pre-commit auto-fixes

* Use `scipy.special.expit` function for Fermi-Dirac distribution, slightly faster, cleaner and no masking required to avoid overflow warnings

* Merge pymatgen master (pt 2)

* pre-commit auto-fixes

* pre-commit auto-fixes

* Fix default `tol` update

* Remove yaml files to force re-eval

* Ensure LF line endings for yamls

---------

Signed-off-by: Seán Kavanagh <51478689+kavanase@users.noreply.github.com>

src/pymatgen/electronic_structure/dos.py

@@ -12,6 +12,7 @@
 from scipy.constants import value as _constant
 from scipy.ndimage import gaussian_filter1d
 from scipy.signal import hilbert
+from scipy.special import expit
 from scipy.stats import wasserstein_distance
 
 from pymatgen.core import Structure, get_el_sp
@@ -75,7 +76,7 @@ def __str__(self) -> str:
 
     def get_interpolated_gap(
         self,
-        tol: float = 0.001,
+        tol: float = 1e-4,
         abs_tol: bool = False,
         spin: Spin | None = None,
     ) -> tuple[float, float, float]:
@@ -120,13 +121,11 @@ def get_interpolated_gap(
         end = get_linear_interpolated_value(terminal_dens, terminal_energies, tol)
         return end - start, end, start
 
-    def get_cbm_vbm(
-        self,
-        tol: float = 0.001,
-        abs_tol: bool = False,
-        spin: Spin | None = None,
-    ) -> tuple[float, float]:
-        """Find the CBM and VBM.
+    def get_cbm_vbm(self, tol: float = 1e-4, abs_tol: bool = False, spin: Spin | None = None) -> tuple[float, float]:
+        """
+        Expects a DOS object and finds the CBM and VBM eigenvalues.
+
+        `tol` may need to be increased for systems with noise/disorder.
 
         Args:
             tol (float): Tolerance in occupations for determining the gap.
@@ -168,13 +167,12 @@ def get_cbm_vbm(
 
         return self.x[i_gap_end], self.x[i_gap_start]
 
-    def get_gap(
-        self,
-        tol: float = 0.001,
-        abs_tol: bool = False,
-        spin: Spin | None = None,
-    ) -> float:
-        """Find the band gap.
+    def get_gap(self, tol: float = 1e-4, abs_tol: bool = False, spin: Spin | None = None) -> float:
+        """
+        Expects a DOS object and finds the band gap, using the determined
+        VBM and CBM eigenvalues.
+
+        `tol` may need to be increased for systems with noise/disorder.
 
         Args:
             tol (float): Tolerance in occupations for determining the gap.
@@ -306,7 +304,7 @@ def get_interpolated_value(self, energy: float) -> dict[Spin, float]:
 
     def get_interpolated_gap(
         self,
-        tol: float = 0.001,
+        tol: float = 1e-4,
         abs_tol: bool = False,
         spin: Spin | None = None,
     ) -> Tuple3Floats:
@@ -348,13 +346,11 @@ def get_interpolated_gap(
 
         return end - start, end, start
 
-    def get_cbm_vbm(
-        self,
-        tol: float = 0.001,
-        abs_tol: bool = False,
-        spin: Spin | None = None,
-    ) -> tuple[float, float]:
-        """Find the conduction band minimum (CBM) and valence band maximum (VBM).
+    def get_cbm_vbm(self, tol: float = 1e-4, abs_tol: bool = False, spin: Spin | None = None) -> tuple[float, float]:
+        """
+        Expects a DOS object and finds the CBM and VBM eigenvalues.
+
+        `tol` may need to be increased for systems with noise/disorder.
 
         Args:
             tol (float): Tolerance in occupations for determining the gap.
@@ -393,7 +389,7 @@ def get_cbm_vbm(
 
     def get_gap(
         self,
-        tol: float = 0.001,
+        tol: float = 1e-4,
         abs_tol: bool = False,
         spin: Spin | None = None,
     ) -> float:
@@ -487,6 +483,7 @@ def __init__(
         ecbm, evbm = self.get_cbm_vbm()
         self.idx_vbm = int(np.argmin(abs(self.energies - evbm)))
         self.idx_cbm = int(np.argmin(abs(self.energies - ecbm)))
+        self.idx_mid_gap = int(self.idx_vbm + (self.idx_cbm - self.idx_vbm) / 2)
         self.A_to_cm = 1e-8
 
         if bandgap:
@@ -502,7 +499,8 @@ def __init__(
             self.energies[idx_fermi:] += (bandgap - (ecbm - evbm)) / 2.0
 
     def get_doping(self, fermi_level: float, temperature: float) -> float:
-        """Calculate the doping (majority carrier concentration) at a given
+        """
+        Calculate the doping (majority carrier concentration) at a given
         Fermi level and temperature. A simple Left Riemann sum is used for
         integrating the density of states over energy & equilibrium Fermi-Dirac
         distribution.
@@ -518,15 +516,15 @@ def get_doping(self, fermi_level: float, temperature: float) -> float:
                 (P-type).
         """
         cb_integral = np.sum(
-            self.tdos[self.idx_cbm :]
-            * f0(self.energies[self.idx_cbm :], fermi_level, temperature)
-            * self.de[self.idx_cbm :],
+            self.tdos[self.idx_mid_gap :]
+            * f0(self.energies[self.idx_mid_gap :], fermi_level, temperature)
+            * self.de[self.idx_mid_gap :],
             axis=0,
         )
         vb_integral = np.sum(
-            self.tdos[: self.idx_vbm + 1]
-            * f0(-self.energies[: self.idx_vbm + 1], -fermi_level, temperature)
-            * self.de[: self.idx_vbm + 1],
+            self.tdos[: self.idx_mid_gap + 1]
+            * f0(-self.energies[: self.idx_mid_gap + 1], -fermi_level, temperature)
+            * self.de[: self.idx_mid_gap + 1],
             axis=0,
         )
         return (vb_integral - cb_integral) / (self.volume * self.A_to_cm**3)
@@ -1568,7 +1566,8 @@ def f0(E: float, fermi: float, T: float) -> float:
     Returns:
         float: The Fermi-Dirac occupation probability at energy E.
     """
-    return 1.0 / (1.0 + np.exp((E - fermi) / (_constant("Boltzmann constant in eV/K") * T)))
+    exponent = (E - fermi) / (_constant("Boltzmann constant in eV/K") * T)
+    return expit(-exponent)  # scipy logistic sigmoid function; expit(x) = 1/(1+exp(-x))
 
 
 def _get_orb_type_lobster(orb: str) -> OrbitalType | None: