import sys
import time
from math import exp, sqrt
import numpy as np
from scipy import integrate # 数値積分関数 integrateを読み込む
from matplotlib import pyplot as plt
"""
金属における有限温度での電子密度計算における数値積分の精度検証を行うスクリプト。
詳細説明:
本スクリプトは、金属の有限温度Tにおける電子密度Nの計算を目的とし、
特に数値積分の精度と性能を確認します。
フェルミエネルギーEF近傍の様々な積分範囲において、
`scipy.integrate.quad` および `scipy.integrate.romberg` 関数を用いて数値積分を実行し、
その結果、計算時間、および可能な場合の解析解との比較を出力します。
状態密度関数 `De(E)` とフェルミ・ディラック分布関数 `fe(E, T, EF)` の積 `Defe(E, T, EF)` を
被積分関数として使用します。
関連リンク: :doc:`N-integration-metal_usage`
"""
#定数
pi = 3.14159265358979323846
h = 6.6260755e-34 # Js";
hbar = 1.05459e-34 # "Js";
c = 2.99792458e8 # m/s";
e = 1.60218e-19 # C";
kB = 1.380658e-23 # JK<sup>-1</sup>";
# Temperature
T = 300.0 # K
# Fermi energy
EF = 5.0 # eV
# EF近傍の nrange * kBTの範囲で積分
nrange = 6.0
# quad/romberg積分の相対積分精度、romberg積分の最大分割数
eps = 1.0e-8
maxorder = 100
#計算時間を計測する繰り返し数
ncycle = 300
# Treat argments
if __name__ == "__main__":
argv = sys.argv
if len(argv) >= 2:
T = float(argv[1])
if len(argv) >= 3:
EF = float(argv[2])
ekBT = e / (kB * T) # eV単位のエネルギーを E/(kBT) に変換する係数
dE = nrange * kB * T / e
[ドキュメント]
def De(E):
"""
自由電子モデルにおける状態密度関数を計算する。
詳細説明:
エネルギー E の平方根に比例する状態密度を返します。
これは3次元自由電子ガスモデルで用いられる一般的な形です。
:param E: エネルギー (float)
:returns: 状態密度の値 (float)
"""
return sqrt(E)
[ドキュメント]
def fe(E, T, EF):
"""
フェルミ・ディラック分布関数を計算する。
詳細説明:
電子が特定のエネルギー準位 E を占める確率を示します。
温度 T が0Kの場合にはステップ関数として振る舞い、E < EF で1.0、E >= EF で0.0を返します。
T > 0Kの場合にはフェルミ・ディラック統計に従います。
:param E: エネルギー (float)
:param T: 温度 (float)
:param EF: フェルミエネルギー (float)
:returns: フェルミ・ディラック分布の値 (float)
"""
global kBT # この行は使用されていないが、元のコードにコメントアウトされていた可能性があるため残す
if T == 0.0:
if E < EF:
return 1.0
else:
return 0.0
return 1.0 / (exp((E - EF) * ekBT) + 1.0)
[ドキュメント]
def Defe(E, T, EF):
"""
状態密度関数とフェルミ・ディラック分布関数の積を計算する。
詳細説明:
数値積分の被積分関数として使用され、特定のエネルギー E における状態密度 `De(E)` に
フェルミ・ディラック分布 `fe(E, T, EF)` を乗じた値を返します。
これにより、エネルギー E における電子の占有状態を考慮した状態密度が表現されます。
:param E: エネルギー (float)
:param T: 温度 (float)
:param EF: フェルミエネルギー (float)
:returns: `De(E) * fe(E, T, EF)` の積の値 (float)
"""
return De(E) * fe(E, T, EF)
[ドキュメント]
def main():
"""
メインプログラムのエントリーポイント。
詳細説明:
コマンドライン引数から温度 T とフェルミエネルギー EF を設定し、
状態密度関数とフェルミ・ディラック分布関数の積を、異なる積分範囲で
`scipy.integrate.quad` および `scipy.integrate.romberg` を用いて数値積分します。
各積分の結果、計算時間、および解析解との比較（可能な場合）を出力し、
数値積分の精度と効率を評価します。
:returns: なし
"""
print("T = ", T, "k EF = ", EF, " eV")
print("")
Ianaly = 2.0 / 3.0 * pow(EF, 1.5)
print("Analytical integration for De(E) from ", 0, " to ", EF, " = ", Ianaly)
ret = integrate.quad(De, 0.0, EF, epsrel = eps)
print("quad : ret=", ret)
ret = integrate.romberg(De, 0.0, EF, rtol = eps, divmax = maxorder)
print("romberg: ret=", ret)
print("")
# 積分範囲
(Emin, Emax) = (0.0, EF + dE)
print("(1) Integration range: ", Emin, " - ", Emax, " eV")
t0 = time.time()
for i in range(ncycle):
ret = integrate.quad(lambda E: Defe(E, T, EF), Emin, Emax, epsrel = eps)
dt = time.time() - t0
print("quad : ret=", ret, " time=", dt, " s")
t0 = time.time()
for i in range(ncycle):
ret = integrate.romberg(lambda E: Defe(E, T, EF), Emin, Emax, rtol = eps, divmax = maxorder)
dt = time.time() - t0
print("romberg: ret=", ret, " time=", dt, " s")
print("")
# 積分範囲
(Emin, Emax) = (0.0, EF - dE)
print("(2) Integration range: ", Emin, " - ", Emax, " eV")
Ianaly = 2.0 / 3.0 * (pow(Emax, 1.5) - pow(Emin, 1.5))
print("Analytical integration for De(E) from ", Emin, " to ", Emax, " = ", Ianaly)
t0 = time.time()
for i in range(ncycle):
ret = integrate.quad(lambda E: Defe(E, T, EF), Emin, Emax, epsrel = eps)
dt = time.time() - t0
print("quad : ret=", ret, " time=", dt, " s")
t0 = time.time()
for i in range(ncycle):
ret = integrate.romberg(lambda E: Defe(E, T, EF), Emin, Emax, rtol = eps, divmax = maxorder)
dt = time.time() - t0
print("romberg: ret=", ret, " time=", dt, " s")
print("")
# 積分範囲
(Emin, Emax) = (EF - dE, EF + dE)
print("(3) Integration range: ", Emin, " - ", Emax, " eV")
t0 = time.time()
for i in range(ncycle):
ret = integrate.quad(lambda E: Defe(E, T, EF), Emin, Emax, epsrel = eps)
dt = time.time() - t0
print("quad : ret=", ret, " time=", dt, " s")
t0 = time.time()
for i in range(ncycle):
ret = integrate.romberg(lambda E: Defe(E, T, EF), Emin, Emax, rtol = eps, divmax = maxorder)
dt = time.time() - t0
print("romberg: ret=", ret, " time=", dt, " s")
print("")
if __name__ == '__main__':
main()