"""
非縮退半導体のEFの温度依存性を二分法で計算。
"""
import sys
import numpy as np
from numpy import sin, cos, tan, pi, exp
from matplotlib import pyplot as plt
# 定数
e = 1.60218e-19 # C
kB = 1.380658e-23 # JK-1
# semiconductor parameters (Default values)
Ev = 0.0 # 価電子帯上端エネルギー eV
Ec = 1.0 # 伝導帯下端エネルギー
Nv = 1.2e19 # 価電子帯有効状態密度
Nc = 2.1e18 # 伝導帯有効状態密度
EA = 0.05 # アクセプター準位
NA = 1.0e15 # アクセプター密度 cm-3
ED = 0.95 # ドナー準位 (Ec - 0.05)
ND = 1.0e16 # ドナー密度
# Temperature range
Tmin = 50 # K
Tmax = 1000
Tstep = 10
# 計算パラメータ
eps = 1.0e-5 # EFの誤差がepsより小さくなったら計算終了
nmaxiter = 200 # 二分法の最大繰り返し数
iprintinterval = 1 # 繰り返し中に途中経過を出力するサイクル数
# --- 物理関数群 ---
[ドキュメント]
def fe(E, EF, T):
"""Fermi-Dirac分布関数"""
global e, kB
if T == 0.0:
return 1.0 if E < EF else 0.0
return 1.0 / (exp((E - EF) * e / kB / T) + 1.0)
[ドキュメント]
def Ne_func(EF, T, Nc_val, Ec_val):
"""伝導帯電子密度 (非縮退近似)"""
global e, kB
if T == 0.0:
return 0.0
return Nc_val * exp(-(Ec_val - EF) * e / kB / T)
[ドキュメント]
def Nh_func(EF, T, Nv_val, Ev_val):
"""価電子帯正孔密度 (非縮退近似)"""
global e, kB
if T == 0.0:
return 0.0
return Nv_val * exp(-(EF - Ev_val) * e / kB / T)
[ドキュメント]
def NDp_func(EF, T, ND_val, ED_val):
"""イオン化ドナー密度"""
return ND_val * (1.0 - fe(ED_val, EF, T))
[ドキュメント]
def NAm_func(EF, T, NA_val, EA_val):
"""イオン化アクセプター密度"""
return NA_val * fe(EA_val, EF, T)
# --- メイン処理 ---
[ドキュメント]
def main():
global Ev, Ec, Nv, Nc, EA, NA, ED, ND
global Tmin, Tmax, Tstep, eps, nmaxiter, iprintinterval
# コマンドライン引数の処理
argv = sys.argv
if len(argv) > 1:
try:
if len(argv) >= 2: EA = float(argv[1])
if len(argv) >= 3: NA = float(argv[2])
if len(argv) >= 4: ED = float(argv[3])
if len(argv) >= 5: ND = float(argv[4])
if len(argv) >= 6: Ec = float(argv[5])
if len(argv) >= 7: Nv = float(argv[6])
if len(argv) >= 8: Nc = float(argv[7])
except ValueError:
print("Error: Arguments must be numeric.")
return
else:
print("usage: python EF-T-semiconductor.py EA NA ED ND Ec Nv Nc")
print(" ex: python EF-T-semiconductor.py 0.05 1.0e15 0.95 1.0e16 1.0 1.2e19 2.1e18")
# 引数なしでもデフォルト値で続行するため exit はしない
nT = int((Tmax - Tmin) / Tstep + 1)
print("Solution of EF by bisection method")
print("Ev=%f Ec=%f eV" % (Ev, Ec))
print("Nv=%e Nc=%e cm-3" % (Nv, Nc))
print("NA=%f cm-3 EA=%f eV" % (NA, EA))
print("ND=%f cm-3 ED=%f eV" % (ND, ED))
print("")
xT = []
xInvT = []
yEF = []
yNe = []
yNh = []
yNAm = []
yNDp = []
print(" T(K) \t EF(eV)\tNe(cm-3)\tNh(cm-3)\tNA+(cm-3)\tND-(cm-3)")
for iT in range(nT):
T = Tmin + iT * Tstep
# 初期範囲として、価電子帯下と伝導帯上を設定
EFmin_val = Ev - 1.0
EFmax_val = Ec + 1.0
# 電荷中性条件: dQ = n + Na- - p - Nd+ = 0 を解く
def get_dQ(target_EF):
return (Ne_func(target_EF, T, Nc, Ec) + NAm_func(target_EF, T, NA, EA) -
Nh_func(target_EF, T, Nv, Ev) - NDp_func(target_EF, T, ND, ED))
dQmin = get_dQ(EFmin_val)
dQmax = get_dQ(EFmax_val)
if dQmin * dQmax > 0.0:
print(f"Error at T={T}: Initial range does not bracket the root.")
continue
# 二分法開始
EFhalf = (EFmin_val + EFmax_val) / 2.0
for i in range(nmaxiter):
EFhalf = (EFmin_val + EFmax_val) / 2.0
dQhalf = get_dQ(EFhalf)
if abs(EFmin_val - EFhalf) < eps and abs(EFmax_val - EFhalf) < eps:
break
if dQmin * dQhalf < 0.0:
EFmax_val = EFhalf
else:
EFmin_val = EFhalf
dQmin = dQhalf
else:
print(f"Warning: T={T} did not converge.")
# 結果の保存（最終ステップの値を使用）
Neh = Ne_func(EFhalf, T, Nc, Ec)
Nhh = Nh_func(EFhalf, T, Nv, Ev)
NAmh = NAm_func(EFhalf, T, NA, EA)
NDph = NDp_func(EFhalf, T, ND, ED)
xT.append(T)
xInvT.append(1000.0 / T)
yEF.append(EFhalf)
yNe.append(Neh)
yNh.append(Nhh)
yNAm.append(NAmh)
yNDp.append(NDph)
print("%8.4f\t%10.6f\t%12.4g\t%12.4g\t%12.4g\t%12.4g" % (T, EFhalf, Neh, Nhh, NAmh, NDph))
# --- グラフ描画 ---
fig = plt.figure(figsize=(12, 5))
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)
# 左：EF vs T
ax1.plot(xT, yEF, label='EF', color='blue')
ax1.axhline(Ev, color='red', linestyle='--', label='Ev')
ax1.axhline(Ec, color='red', linestyle='--', label='Ec')
ax1.set_xlabel("T (K)")
ax1.set_ylabel("EF (eV)")
ax1.set_ylim([Ev - 0.2, Ec + 0.2])
ax1.legend()
# 右：Carrier Densities vs 1000/T (Arrhenius plot)
ax2.plot(xInvT, yNe, label='n (electrons)')
ax2.plot(xInvT, yNh, label='p (holes)')
ax2.plot(xInvT, yNAm, label='Na- (ionized acc.)', linestyle=':')
ax2.plot(xInvT, yNDp, label='Nd+ (ionized don.)', linestyle=':')
ax2.set_yscale("log")
ax2.set_xlabel("1000/T (K^-1)")
ax2.set_ylabel("Density (cm^-3)")
all_y = yNe + yNh + yNAm + yNDp
if all_y:
ax2.set_ylim([1e8, max(all_y) * 2])
ax2.legend()
plt.subplots_adjust(wspace=0.4)
plt.show(block=False)
input("Press ENTER to exit>>")
if __name__ == "__main__":
main()