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 の計算。
このスクリプトは、自由電子モデルに基づき、金属中の電子の有限温度における電子密度を計算します。
状態密度関数、フェルミ・ディラック分布関数を用いて、電子密度を数値積分で求めます。
また、計算結果と各関数をグラフにプロットして可視化します。
:doc:`EF-N-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>";
me = 9.1093897e-31 # kg";
# Temperature
T = 300.0 # K
# Fermi energy
EF = 5.0 # eV
# EF近傍の nrange * kBTの範囲で積分
nrange = 6.0
# quad/romberg積分の相対積分精度
eps = 1.0e-8
fontsize = 14
# 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
S = 1.0 / 2.0
D0 = (2.0 * S + 1.0) * 2.0 * pi * pow(2.0 * me, 1.5) / h**3 * pow(e, 1.5) * 1.0e-6 # eV/cm^3
[ドキュメント]
def De(E):
"""
自由電子モデルに基づく状態密度関数を計算します。
この関数は、エネルギー E における電子の状態密度 D_e(E) を返します。
D_e(E) = D0 * sqrt(E) の形式で表現されます。
:param E: float: エネルギー (eV)
:returns: float: エネルギー E における状態密度 (states/cm^3/eV)
"""
return D0 * sqrt(E)
[ドキュメント]
def fe(E, T, EF):
"""
フェルミ・ディラック分布関数を計算します。
与えられたエネルギー E、温度 T、フェルミエネルギー EF における電子の占有確率を計算します。
温度が0Kの場合、フェルミエネルギー未満では1、以上では0となるステップ関数として機能します。
:param E: float: エネルギー (eV)
:param T: float: 温度 (K)
:param EF: float: フェルミエネルギー (eV)
:returns: float: フェルミ・ディラック分布関数の値
"""
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):
"""
状態密度とフェルミ・ディラック分布関数の積を計算します。
この関数は、電子の総数（電子密度）を計算するための被積分関数 `D_e(E) * f_{FD}(E, T, E_F)` を提供します。
:param E: float: エネルギー (eV)
:param T: float: 温度 (K)
:param EF: float: フェルミエネルギー (eV)
:returns: float: `D_e(E) * f_{FD}(E, T, E_F)` の値 (states/cm^3/eV)
"""
return De(E) * fe(E, T, EF)
[ドキュメント]
def main():
"""
金属の電子密度を計算し、その結果をプロットします。
この関数は以下の処理を実行します。
1. 0Kにおける電子密度を解析的および数値的に計算し、出力します。
2. 有限温度における電子密度を数値的に計算し、出力します。
3. 状態密度 D_e(E)、0Kおよび有限温度での D_e(E) * f_{FD}(E,T)、
フェルミ・ディラック分布関数 f_{FD}(E,T) をエネルギーの関数としてプロットします。
4. プロット結果を画面に表示します。
コマンドライン引数により、温度 T およびフェルミエネルギー EF を変更できます。
:returns: None
"""
print("")
print(f"at T = 0 K, EF = {EF} eV")
Na = 2.0 / 3.0 * pow(EF, 1.5)
print(f" Analytical integration: N={Na} cm^-3")
N = integrate.quad(De, 0.0, EF, epsrel = eps)
print(f" Numerical integration: N={N[0]} cm^-3")
print("")
print(f"at T = {T} K, EF = {EF} eV")
# 積分範囲
Emax = EF + dE
N = integrate.quad(lambda E: Defe(E, T, EF), 0.0, Emax, epsrel = eps)
print(f" Numerical integration: N={N[0]} cm^-3")
# plot
xE = np.arange(0.0, Emax, 0.01)
yDe = [De(E) for E in xE]
yfe = [fe(E, T, EF) for E in xE]
yDe0 = [Defe(E, 0.0, EF) for E in xE]
yDefe = [Defe(E, T, EF) for E in xE]
fig, axis = plt.subplots(1, 1, figsize = (8, 6))
axis.tick_params(labelsize = fontsize)
# fFE(E)を見やすくするため、yDefeのスケールに合わせる
yfe = max(yDefe) * np.array(yfe)
axis.plot(xE, yDe, label = r'$D_e(E)$', linestyle = '-', color = 'black')
axis.plot(xE, yDe0, label = r'$D_e(E) \cdot f_{FD}(E,0)$', linestyle = '-', color = 'blue')
axis.plot(xE, yfe, label = r'$f_{FD}(E,T,E_F)$', linestyle = '-', color = 'red')
axis.plot(xE, yDefe, label = r'$D_e(E) \cdot f_{FD}(E,T)$', linestyle = '-', color = 'darkgreen')
axis.axvline(EF, label = r'$E_F$', linestyle = 'dashed', color = 'red', linewidth = 0.5)
axis.axhline(0.0, linestyle = 'dashed', color = 'red', linewidth = 0.5)
axis.fill_between(xE, yDefe, 0.0, facecolor='lime', alpha=0.5)
axis.set_xlabel(r'$E - E_C$ (eV)', fontsize = fontsize)
axis.set_ylabel(r'DOS (states/cm$^3$)', fontsize = fontsize)
axis.legend(fontsize = fontsize)
plt.tight_layout()
plt.show()
if __name__ == '__main__':
main()