code: def solve_for_fermi_energy(self, temperature, chemical_potentials, bulk_dos):
          
  
          fdos = FermiDos(bulk_dos, bandgap=self.band_gap)
  
          def _get_total_q(ef):
  
              qd_tot = sum([
                  d['charge'] * d['conc']
                  for d in self.defect_concentrations(
                      chemical_potentials=chemical_potentials, temperature=temperature, fermi_level=ef)
              ])
              qd_tot += fdos.get_doping(fermi=ef + self.vbm, T=temperature)
              return qd_tot
  
          return bisect(_get_total_q, -1., self.band_gap + 1.)
docs: Solve for the Fermi energy self-consistently as a function of T
          and p_O2
          Observations are Defect concentrations, electron and hole conc
          Args:
              bulk_dos: bulk system dos (pymatgen Dos object)
              gap: Can be used to specify experimental gap.
                  Will be useful if the self consistent Fermi level
                  is > DFT gap
          Returns:
              Fermi energy