Minor update to make FermiDos more robust

src/pymatgen/electronic_structure/dos.py

@@ -156,14 +156,13 @@ def get_cbm_vbm(self, tol: float = 1e-4, abs_tol: bool = False, spin: Spin | Non
 
         # Work backwards until tolerance is reached
         i_gap_start = i_fermi
-        while i_gap_start >= 1 and tdos[i_gap_start - 1] <= tol:
+        while i_gap_start >= 1 and tdos[i_gap_start] <= tol:
             i_gap_start -= 1
 
         # Work forwards until tolerance is reached
-        i_gap_end = i_gap_start
+        i_gap_end = i_gap_start + 1
         while i_gap_end < tdos.shape[0] and tdos[i_gap_end] <= tol:
             i_gap_end += 1
-        i_gap_end -= 1
 
         return self.x[i_gap_end], self.x[i_gap_start]
 
@@ -378,14 +377,13 @@ def get_cbm_vbm(self, tol: float = 1e-4, abs_tol: bool = False, spin: Spin | Non
 
         # Work backwards until tolerance is reached
         i_gap_start = i_fermi
-        while i_gap_start >= 1 and tdos[i_gap_start - 1] <= tol:
+        while i_gap_start >= 1 and tdos[i_gap_start] <= tol:
             i_gap_start -= 1
 
         # Work forwards until tolerance is reached
-        i_gap_end = i_gap_start
+        i_gap_end = i_gap_start + 1
         while i_gap_end < tdos.shape[0] and tdos[i_gap_end] <= tol:
             i_gap_end += 1
-        i_gap_end -= 1
 
         return self.energies[i_gap_end], self.energies[i_gap_start]
 
@@ -518,15 +516,15 @@ def get_doping(self, fermi_level: float, temperature: float) -> float:
                 (P-type).
         """
         cb_integral = np.sum(
-            self.tdos[self.idx_mid_gap :]
-            * f0(self.energies[self.idx_mid_gap :], fermi_level, temperature)
-            * self.de[self.idx_mid_gap :],
+            self.tdos[max(self.idx_mid_gap, self.idx_vbm + 1) :]
+            * f0(self.energies[max(self.idx_mid_gap, self.idx_vbm + 1) :], fermi_level, temperature)
+            * self.de[max(self.idx_mid_gap, self.idx_vbm + 1) :],
             axis=0,
         )
         vb_integral = np.sum(
-            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],
+            self.tdos[: min(self.idx_mid_gap, self.idx_cbm - 1) + 1]
+            * f0(-self.energies[: min(self.idx_mid_gap, self.idx_cbm - 1) + 1], -fermi_level, temperature)
+            * self.de[: min(self.idx_mid_gap, self.idx_cbm - 1) + 1],
             axis=0,
         )
         return (vb_integral - cb_integral) / (self.volume * self.A_to_cm**3)

tests/electronic_structure/test_dos.py

@@ -13,7 +13,7 @@
 
 from pymatgen.core import Element, Structure
 from pymatgen.electronic_structure.core import Orbital, OrbitalType, Spin
-from pymatgen.electronic_structure.dos import DOS, CompleteDos, FermiDos, LobsterCompleteDos
+from pymatgen.electronic_structure.dos import DOS, CompleteDos, Dos, FermiDos, LobsterCompleteDos
 from pymatgen.util.testing import TEST_FILES_DIR, PymatgenTest
 
 TEST_DIR = f"{TEST_FILES_DIR}/electronic_structure/dos"
@@ -103,6 +103,26 @@ def test_as_dict(self):
         assert isinstance(dos_dict["densities"]["1"][0], float)
         assert not isinstance(dos_dict["densities"]["1"][0], np.float64)
 
+    def test_get_vbm_cbm_doping(self):
+        dos = Dos(
+            energies=np.array([0.0, 0.5, 1.0, 1.5, 2.0]),
+            densities={
+                Spin.up: np.array([1.0, 2.0, 0.0, 3.0, 4.0]),
+                Spin.down: np.array([0.5, 1.0, 0.0, 1.5, 2.0]),
+            },
+            efermi=0.8,
+        )
+        fermi_dos = FermiDos(
+            dos,
+            structure=self.dos.structure,
+        )
+        assert fermi_dos.get_cbm_vbm() == (1.5, 0.5)
+        assert np.isclose(
+            fermi_dos.get_doping(fermi_level=1.0, temperature=300),
+            -1385561583858093.5,  # <0 because e doping; greater DOS in CBM than VBM here, and efermi set to mid-gap
+            rtol=1e-3,
+        )
+
 
 class TestCompleteDos(TestCase):
     def setUp(self):