"""
半導体中のフェルミ準位を二分法で計算するモジュール。
このスクリプトは、電荷中性条件を満たすフェルミ準位（EF）を二分法を用いて数値的に決定します。
電子、正孔、イオン化されたドナーおよびアクセプターのキャリア密度を計算し、
それらの合計電荷がゼロになるEFを見つけます。
:doc:`equation-bisection_usage`
"""
import csv
import numpy as np
from numpy import sin, cos, tan, pi, exp
import sys
# physical parameeters
e = 1.60218e-19 # C";
kB = 1.380658e-23 # JK^-1";
T = 300.0
ekBT = e / kB / T
# semiconductor parameters
Ev = 0.0
Ec = 1.0
Nv = 1.2e19
Nc = 2.1e18
EA = 0.05
NA = 1.0e15
ED = Ec - 0.05
ND = 1.0e16
# parameters
EFmin = Ev
EFmax = Ec
eps = 1.0e-5
nmaxiter = 200
iprintiterval = 1
# Fermi-Dirac function
[ドキュメント]
def fe(E, EF, T):
"""
フェルミ・ディラック分布関数を計算する。
電子があるエネルギー準位Eを占有する確率をフェルミ・ディラック統計に基づいて計算します。
計算にはグローバル変数`ekBT` (e / kB / T) を使用します。
:param E: float: エネルギー準位 (eV)。
:param EF: float: フェルミ準位 (eV)。
:param T: float: 温度 (K)。
:returns: float: エネルギー準位Eが占有される確率。
"""
global Nc, Ec, kB # これらのグローバル変数は実際にはこの関数内で直接は使用されていないが、元のコードに存在するため残す。
return 1.0 / (exp((E - EF) * ekBT) + 1.0)
# electron density
[ドキュメント]
def Ne(EF, T):
"""
伝導帯における電子密度を計算する。
伝導帯の下端Ecとフェルミ準位EFから電子密度を計算します。
ここでは古典近似（Boltzmann近似）が使用されており、グローバル変数`Nc`, `Ec`, `ekBT`を使用します。
:param EF: float: フェルミ準位 (eV)。
:param T: float: 温度 (K)。
:returns: float: 電子密度 (cm^-3)。
"""
global Nc, Ec, kB # これらのグローバル変数は実際にはこの関数内で直接は使用されていないが、元のコードに存在するため残す。
return Nc * exp(-(Ec - EF) * ekBT)
# hole density
[ドキュメント]
def Nh(EF, T):
"""
価電子帯における正孔密度を計算する。
価電子帯の上端Evとフェルミ準位EFから正孔密度を計算します。
ここでは古典近似（Boltzmann近似）が使用されており、グローバル変数`Nv`, `Ev`, `ekBT`を使用します。
:param EF: float: フェルミ準位 (eV)。
:param T: float: 温度 (K)。
:returns: float: 正孔密度 (cm^-3)。
"""
global Nv, Ev, kB # これらのグローバル変数は実際にはこの関数内で直接は使用されていないが、元のコードに存在するため残す。
return Nv * exp(-(EF - Ev) * ekBT)
# ionized donor density
[ドキュメント]
def NDp(EF, T):
"""
イオン化されたドナー密度を計算する。
ドナー準位EDにある電子が伝導帯へ励起され、イオン化したドナーの密度を計算します。
フェルミ・ディラック分布関数`fe`とグローバル変数`ND`, `ED`を使用します。
:param EF: float: フェルミ準位 (eV)。
:param T: float: 温度 (K)。
:returns: float: イオン化ドナー密度 (cm^-3)。
"""
global ND, ED, kB # これらのグローバル変数は実際にはこの関数内で直接は使用されていないが、元のコードに存在するため残す。
return ND * (1.0 - fe(ED, EF, T))
# ionized acceptor density
[ドキュメント]
def NAm(EF, T):
"""
イオン化されたアクセプター密度を計算する。
アクセプター準位EAに電子が捕獲され、イオン化したアクセプターの密度を計算します。
フェルミ・ディラック分布関数`fe`とグローバル変数`NA`, `EA`を使用します。
:param EF: float: フェルミ準位 (eV)。
:param T: float: 温度 (K)。
:returns: float: イオン化アクセプター密度 (cm^-3)。
"""
global NA, EA, kB # これらのグローバル変数は実際にはこの関数内で直接は使用されていないが、元のコードに存在するため残す。
return NA * fe(EA, EF, T)
[ドキュメント]
def main():
"""
半導体中のフェルミ準位を二分法で計算し、結果を出力するメイン関数。
この関数は、`EFmin`と`EFmax`で定義された範囲内で、電荷中性条件
(電子密度 + イオン化アクセプター密度 - 正孔密度 - イオン化ドナー密度 = 0)
を満たすフェルミ準位 (EF) を二分法を用いて探索します。
探索の進行状況と最終結果を標準出力に表示します。
初期条件が不適切であったり、最大反復回数を超過した場合はエラーメッセージを出力します。
:returns: int: 計算の成功状態。1は成功、0は失敗または初期条件エラーを示す。
"""
global EFmin, EFmax, eps, nmaxiter, iprintiterval
print("Solution of EF by bisection method")
print("")
# 初期フェルミ準位範囲での電荷不均衡を計算
dQmin = Ne(EFmin, T) + NAm(EFmin, T) - Nh(EFmin, T) - NDp(EFmin, T)
dQmax = Ne(EFmax, T) + NAm(EFmax, T) - Nh(EFmax, T) - NDp(EFmax, T)
print(" EFmin = {:12.8f} dQmin = {:12.4g}".format(EFmin, dQmin))
print(" EFmax = {:12.8f} dQmax = {:12.4g}".format(EFmax, dQmax))
if dQmin * dQmax > 0.0:
print("Error: Initial Emin and Emax should be chosen as dQmin * dQmax < 0")
return 0
for i in range(nmaxiter):
EFhalf = (EFmin + EFmax) / 2.0
# 中間点での各キャリア密度と電荷不均衡を計算
Neh = Ne(EFhalf, T)
NAmh = NAm(EFhalf, T)
Nhh = Nh(EFhalf, T)
NDph = NDp(EFhalf, T)
dQhalf = Neh + NAmh - Nhh - NDph
print(" Iter {}: EFhalf = {:12.8f} dQhalf = {:12.4g}".format(i, EFhalf, dQhalf))
print(" Ne={:10.4e} Nh={:10.4e} NA-={:10.4e} ND+={:10.4e} dQ={:10.4e}".format(Neh, Nhh, NAmh, NDph, dQhalf))
# 収束判定
if abs(EFmin - EFhalf) < eps and abs(EFmax - EFhalf) < eps:
print(" Success: Convergence reached at EF = {}".format(EFhalf))
return 1
# 二分法の更新ステップ
if dQmin * dQhalf < 0.0:
EFmax = EFhalf
dQmax = dQhalf
else:
EFmin = EFhalf
dQmin = dQhalf
else:
# 最大反復回数に達した場合
print(" Failed: Convergence did not reach")
return 0
if __name__ == "__main__":
main()