# Integration of the Fermi distribution using Python

I want to calculate the carrier concentration of my semiconductor using this equation:

$$n(x) = \frac{m^*}{\pi\hbar^2}\int_{E_k}^{\infty}\frac{1}{1+\exp\left(\frac{E-E_f}{k_BT}\right)} \mathrm{d}E$$

I'm using this straightforward scipy approach:

import numpy as np
import scipy.constants as phys
import scipy.integrate as integrate

eigenvalue = [0.9 * phys.electron_volt, 1.3 * phys.electron_volt]
fermi = 1.0 * phys.electron_volt
T = 300

def fermi_integral(E, fermi, T):
return 1 / (1 + np.exp((E - fermi) / (phys.Boltzmann * T)))

for i in range(len(eigenvalue)):
result = integrate.quad(fermi_integral, eigenvalue[i], np.inf, args=(fermi, T))
print(result)


However, I'm running into

RuntimeWarning: overflow encountered in exp
return 1 / (1 + np.exp((E - fermi) / (phys.Boltzmann * T)))


and my result is always (0.0, 0.0)

I guess I have to use another approach, but I'm stuck and I hope you can give me some helpful input.

• Numerical integration methods usually assume that your function changes meaningfully on the scale of the interval you give it, which is something like $(0,1)$ here. But your function is almost identically zero from $100\times k_B T$ onwards, so trying to evaluate it at a default value of e.g. $\frac12$ causes an error. If there were no overflow, it would still give a bad result because the function values would be indistinguishable from zero on all the default values quad tries at first. That's why rescaling can fix this as in nicoguaro's answer. – Kirill Jul 17 '18 at 18:20

One of your problems is the system of units that you are using. Just changing the units improves the results

import numpy as np
import scipy.integrate as integrate

eigenvalue = [0.9, 1.3]
fermi = 1.0
T = 300
kB = 8.6173303e-5

def fermi_integral(E, fermi, T):
return 1 / (1 + np.exp((E - fermi) / (kB * T)))

for i in range(len(eigenvalue)):
args=( fermi, T))
print(result)


These are the results

(0.10053464900138948, 1.2260194855049155e-09)
(2.3589139312840435e-07, 4.1489077967047375e-11)


I still get the overflow problem, but that is probably because you are evaluating the functions for really large numbers. I just replaced np.inf for 10 and obtained the same results.

### Edit

Taking into account @gammatester suggestion the overflow warning does not appear anymore. The following works for me.

import numpy as np
import scipy.integrate as integrate

eigenvalue = [0.9, 1.3]
fermi = 1.0
T = 300
kB = 8.6173303e-5

def fermi_integral(E, fermi, T):
if E < fermi:
return 1 / (1 + np.exp((E - fermi) / (kB * T)))
else:
return np.exp(-(E - fermi) / (kB * T)) / (1 + np.exp(-(E - fermi) / (kB * T)))

for i in range(len(eigenvalue)):
args=( fermi, T))
print(result)

• Or you can transform the integand so that no overflow occurs for large $x:$ $$\frac{1}{1+e^x}=\frac{e^{-x}}{1+e^{-x}}$$ – gammatester Jul 17 '18 at 17:38
• @gammatester, that works. Although, one has to find what "large $x$" means. In this case it seems that 300 is large enough. – nicoguaro Jul 17 '18 at 18:38
• You can always use it, it is valid for all $x;$ but it would overflow for large negative $x$. I would use it if $x>0$ and the original for $x\le0,$ i.e. if $E>E_f$. – gammatester Jul 17 '18 at 18:56

I know this is already answered satisfactorily, but let me just add this gentle reminder to check whether the integral can be done analytically, before crunching it with a numerical method! In this case the indefinite integral of $$f(E) = \frac{1}{1+\exp[(E-E_f)/k_B T]}$$ is this: $$F(E) = -k_B T \ln \left\{ 1+\exp[-(E-E_f)/k_B T] \right\}$$ and $F(\infty)=0$, so the desired answer is $-F(E_k)$. The following python code

import numpy as np

eigenvalue = np.array([0.9, 1.3])
fermi = 1.0
T = 300
kB = 8.6173303e-5

result = kB*T*np.log(1+np.exp(-(eigenvalue-fermi)/(kB*T)))
print(result)



gives the exact result(s), for comparison with the numerical one(s).

[1.00534649e-01 2.35891393e-07]