计算科学 - 使用 Python 集成费米分布 - 吾爱随笔录

使用 Python 集成费米分布

计算科学 Python scipy 一体化

2021-12-01 23:26:36

我想用这个等式计算我的半导体的载流子浓度：

n (x) = \frac{m^{*}}{π ℏ^{2}} \int_{E_{k}}^{\infty} \frac{1}{1 + \exp (\frac{E - E_{f}}{k_{B} T})} d E

$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$

我正在使用这种简单的 scipy 方法：

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)

但是，我遇到了

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

我的结果总是(0.0, 0.0)

我想我必须使用另一种方法，但我被困住了，我希望你能给我一些有用的意见。

3个回答

您的问题之一是您使用的单位制。只需更改单位即可改善结果

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)

这些是结果

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

我仍然遇到溢出问题，但这可能是因为您正在评估非常大的数字的函数。我刚刚更换np.inf了 10 个并获得了相同的结果。

编辑

考虑到@gammatester 的建议，溢出警告不再出现。以下对我有用。

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)

我知道这已经得到了令人满意的回答，但让我只是添加这个温和的提醒，以检查积分是否可以解析地完成，然后再用数值方法处理它！在这种情况下，不定积分

f (E) = \frac{1}{1 + \exp [(E - E_{f}) / k_{B} T]}

$f(E) = \frac{1}{1+\exp[(E-E_f)/k_B T]}$ 这是：

F (E) = - k_{B} T \ln {1 + \exp [- (E - E_{f}) / k_{B} T]}

$F(E) = -k_B T \ln \left\{ 1+\exp[-(E-E_f)/k_B T] \right\}$ 和

F (\infty) = 0

$F(\infty)=0$ ，所以想要的答案是

- F (E_{k})

$-F(E_k)$ . 以下python代码

将 numpy 导入为 np

特征值 = np.array([0.9, 1.3])
费米 = 1.0
T = 300
kB = 8.6173303e-5

结果 = kB*T*np.log(1+np.exp(-(特征值-费米)/(kB*T)))
打印（结果）

给出准确的结果，以便与数值结果进行比较。

[1.00534649e-01 2.35891393e-07]

我主张将scipy.stats.logistic.sf用作费米-狄拉克分布。

Fermi-Dirac 等价于后勤生存函数（也称为complementary cumulative distribution function）。
使用上述方法的好处是，您可以立即访问通过 scipy 中的分布函数接口访问的许多不同方法 + 它通过检查输入参数来处理许多极端情况。就像这里scipy.special.expit用来计算这些极端情况的情况一样。

很好的是，人们可以从 FD 分布中抽取一些能量周围的样本。

其它你可能感兴趣的问题

上一篇hp-FEM 简介下一篇需要一个简单的网格格式（FEA）和一个生成网格的工具