3
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – Kirill Jul 17 '18 at 18:20
7
$\begingroup$

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)):
    result = integrate.quad(fermi_integral, eigenvalue[i], np.inf,
                            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)):
    result = integrate.quad(fermi_integral, eigenvalue[i], np.inf,
                            args=( fermi, T))
    print(result)
$\endgroup$
  • 4
    $\begingroup$ 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}}$$ $\endgroup$ – gammatester Jul 17 '18 at 17:38
  • 1
    $\begingroup$ @gammatester, that works. Although, one has to find what "large $x$" means. In this case it seems that 300 is large enough. $\endgroup$ – nicoguaro Jul 17 '18 at 18:38
  • 4
    $\begingroup$ 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$. $\endgroup$ – gammatester Jul 17 '18 at 18:56
5
$\begingroup$

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]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.