diff --git a/src/nomad_simulations/properties/band_structure.py b/src/nomad_simulations/properties/band_structure.py
index 9f171388..1f269ec5 100644
--- a/src/nomad_simulations/properties/band_structure.py
+++ b/src/nomad_simulations/properties/band_structure.py
@@ -121,7 +121,6 @@ class ElectronicEigenvalues(BaseElectronicEigenvalues):
     def __init__(
         self, m_def: Section = None, m_context: Context = None, **kwargs
     ) -> None:
-        # ! `n_bands` need to be set up during initialization of the class
         super().__init__(m_def, m_context, **kwargs)
         self.name = self.m_def.name
 
@@ -207,20 +206,25 @@ def resolve_homo_lumo_eigenvalues(
         sorted_value = sorted_value.magnitude
 
         # Binary search ot find the transition point between `occupation = 2` and `occupation = 0`
-        # TODO add tolerance
+        tolerance = 1e-3  # TODO add tolerance from config fields
         homo = None
         lumo = None
         low_occ = 0
         high_unocc = sorted_value.shape[-1] - 1
         while low_occ <= high_unocc:
             mid = (low_occ + high_unocc) // 2
-            # Check if `occupation` changes at this index
-            if sorted_occupation[mid] > 0 and sorted_occupation[mid + 1] == 0:
+            # check if occupation[mid] and [mid+1] is 0
+            if sorted_occupation[mid] > 0 and (
+                sorted_occupation[mid + 1] >= -tolerance
+                and sorted_occupation[mid + 1] <= tolerance
+            ):
                 homo = sorted_value[mid] * sorted_value_unit
                 lumo = sorted_value[mid + 1] * sorted_value_unit
                 break
+            # check if the occupation[mid] is finite but [mid+1] is as well
             elif sorted_occupation[mid] > 0:
                 low_occ = mid + 1
+            # check if the occupation[mid] is 0
             else:
                 high_unocc = mid - 1
 
@@ -246,6 +250,15 @@ def extract_band_gap(self) -> Optional[ElectronicBandGap]:
         return band_gap
 
     def extract_fermi_surface(self, logger: BoundLogger) -> Optional[FermiSurface]:
+        """
+        Extract the Fermi surface for metal systems and using the `FermiLevel.value`.
+
+        Args:
+            logger (BoundLogger): The logger to log messages.
+
+        Returns:
+            (Optional[FermiSurface]): The extracted Fermi surface section to be stored in `Outputs`.
+        """
         # Check if the system has a finite band gap
         homo, lumo = self.highest_occupied, self.lowest_unoccupied
         if (homo and lumo) and (lumo - homo).magnitude > 0:
@@ -260,7 +273,7 @@ def extract_fermi_surface(self, logger: BoundLogger) -> Optional[FermiSurface]:
         fermi_level_value = fermi_level.value.magnitude
 
         # Extract values close to the `fermi_level.value`
-        tolerance = 1e-8
+        tolerance = 1e-8  # TODO change this for a config field
         fermi_indices = np.logical_and(
             self.value.magnitude >= (fermi_level_value - tolerance),
             self.value.magnitude <= (fermi_level_value + tolerance),
@@ -309,7 +322,6 @@ class ElectronicBandStructure(ElectronicEigenvalues):
     def __init__(
         self, m_def: Section = None, m_context: Context = None, **kwargs
     ) -> None:
-        # ! `n_bands` need to be set up during initialization of the class
         super().__init__(m_def, m_context, **kwargs)
         self.name = self.m_def.name