金属の有限温度 T における電子密度 N の計算。

Download script from EF-N-metal.py
Related files:

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 の計算。
"""


#定数
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-1";
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
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):
    return D0 * sqrt(E)

# FE分布関数
def fe(E, T, EF):
    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):
    return De(E) * fe(E, T, EF)



def main():
    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()