フェルミ・ディラック関数で電子密度計算をscipy.integrateでかんばろうとして、
z = lambda x: 1 / (1 + np.exp((x - E_F)/(k_b * T)))
total_n, total_n_err = spi.quad(z, 0, np.inf)
とすると、
0.0 0.0
となり実質的なエラー。docsに記述があり、どうやら計算アルゴリズム的に、積分範囲を適当に設定すると、フェルミ・ディラックみたいなほぼ0の値が裾を引くタイプの関数では想定通りの結果がでないらしい。
つまり、狙っている部分を囲むような有限積分で近似するのが方法としてよいようで。
フェルミ・ディラックの場合は$E_F$の周りに$k_b T$の何倍かを足しておいて、許容できる相対誤差まで動かしてみるのがよさそうだ。