5
$\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.

Improve this question
$\endgroup$
1
  • 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
    Commented Jul 17, 2018 at 18:20

3 Answers 3

Reset to default
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)
Improve this answer
$\endgroup$
3
  • 5
    $\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
    Commented Jul 17, 2018 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 ♦
    Commented Jul 17, 2018 at 18:38
  • 5
    $\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
    Commented Jul 17, 2018 at 18:56
6
$\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]
Improve this answer
$\endgroup$
2
$\begingroup$

I would advocate on using the scipy.stats.logistic.sf as the Fermi-Dirac distribution.

The Fermi-Dirac is equivalent to the logistics survival function (otherwise known as complementary cumulative distribution function).
The benefits of using the above is that you'll have immediate access to many different methods accessible through distribution function interfaces in scipy + it handles many corner cases by examining input arguments. As is the case here where it uses scipy.special.expit to calculate these corner cases.

What is nice, is that one can draw samples from the FD distribution around some energy.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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