code: def defect_concentration(self, chemical_potentials, temperature=300, fermi_level=0.0):
          
          n = self.multiplicity * 1e24 / self.defect.bulk_structure.volume
          conc = n * np.exp(-1.0 * self.formation_energy(chemical_potentials, fermi_level=fermi_level) /
                            (kb * temperature))
  
          return conc
docs: Get the defect concentration for a temperature and Fermi level.
          Args:
              temperature:
                  the temperature in K
              fermi_level:
                  the fermi level in eV (with respect to the VBM)
          Returns:
              defects concentration in cm^-3