import csv
import numpy as np
from numpy import sin, cos, tan, pi, exp
import sys
"""
概要: ブレント法を用いて半導体のフェルミ準位を計算するスクリプト。
詳細説明:
このスクリプトは、電子、正孔、イオン化ドナー、イオン化アクセプタの電荷バランス方程式
(Ne + NAm - Nh - NDp = 0) の根をブレント法 (Brent's method) を用いて探索します。
ブレント法は、二分法、割線法、逆二次補間法を組み合わせたロバストな求根アルゴリズムです。
アルゴリズムおよびその仮想的な実装については、以下のWikipedia記事を参照しています。
https://ja.wikipedia.org/wiki/%E3%83%96%E3%83%AC%E3%83%B3%E3%83%88%E6%B3%95
関連リンク:
:doc:`equation-brent_usage`
"""
# physical parameeters
e = 1.60218e-19 # C";
kB = 1.380658e-23 # JK<sup>-1</sup>";
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
#delta = (EFmax - EFmin) * 1.0e-1
delta = eps
nmaxiter = 100
iprintinterval = 1
[ドキュメント]
def fe(E, EF, T):
"""
概要: フェルミ・ディラック分布関数を計算する。
詳細説明:
与えられたエネルギーEにおけるフェルミ・ディラック分布関数の値を計算します。
:param E: エネルギー。
:type E: float
:param EF: フェルミ準位。
:type EF: float
:param T: 温度 (K)。
:type T: float
:returns: フェルミ・ディラック分布関数の値。
:rtype: float
"""
global Nc, Ec, kB
return 1.0 / (exp((E - EF) * ekBT) + 1.0)
[ドキュメント]
def Ne(EF, T):
"""
概要: 電子密度を計算する。
詳細説明:
伝導帯における電子密度を計算します。
この計算は、古典的なボルツマン近似に基づいています。
:param EF: フェルミ準位。
:type EF: float
:param T: 温度 (K)。
:type T: float
:returns: 電子密度。
:rtype: float
"""
global Nc, Ec, kB
return Nc * exp(-(Ec - EF) * ekBT)
[ドキュメント]
def Nh(EF, T):
"""
概要: 正孔密度を計算する。
詳細説明:
価電子帯における正孔密度を計算します。
この計算は、古典的なボルツマン近似に基づいています。
:param EF: フェルミ準位。
:type EF: float
:param T: 温度 (K)。
:type T: float
:returns: 正孔密度。
:rtype: float
"""
global Nv, Ev, kB
return Nv * exp(-(EF - Ev) * ekBT)
[ドキュメント]
def NDp(EF, T):
"""
概要: イオン化ドナー密度を計算する。
詳細説明:
ドナー準位から電子が放出されイオン化したドナーの密度を計算します。
:param EF: フェルミ準位。
:type EF: float
:param T: 温度 (K)。
:type T: float
:returns: イオン化ドナー密度。
:rtype: float
"""
global ND, ED, kB
return ND * (1.0 - fe(ED, EF, T))
[ドキュメント]
def NAm(EF, T):
"""
概要: イオン化アクセプタ密度を計算する。
詳細説明:
アクセプタ準位が電子を受け入れイオン化したアクセプタの密度を計算します。
:param EF: フェルミ準位。
:type EF: float
:param T: 温度 (K)。
:type T: float
:returns: イオン化アクセプタ密度。
:rtype: float
"""
global NA, EA, kB
return NA * fe(EA, EF, T)
[ドキュメント]
def main():
"""
概要: ブレント法を用いて半導体のフェルミ準位を計算し、結果を出力する。
詳細説明:
半導体中の電荷中性条件 (Ne + NAm - Nh - NDp = 0) を満たすフェルミ準位をブレント法で探索します。
初期値の選択、反復計算の実行、収束判定、および途中経過と最終結果の出力を行います。
ブレント法のアルゴリズムに従い、二分法、割線法、逆二次補間法を組み合わせて使用します。
:returns: 計算が成功し収束した場合は 1、初期条件エラーまたは非収束の場合は 0 を返す。
:rtype: int
"""
global EFmin, EFmax, delta, eps, nmaxiter, iprintinterval
print("Solution of EF by Brent's method")
print("")
[EFa, EFb] = [EFmin, EFmax]
dQa = Ne(EFa, T) + NAm(EFa, T) - Nh(EFa, T) - NDp(EFa, T)
dQb = Ne(EFb, T) + NAm(EFb, T) - Nh(EFb, T) - NDp(EFb, T)
if dQa * dQb > 0.0:
print("Error: Initial Emin and Emax should be chosen as dQmin * dQmax < 0")
return 0
if abs(dQa) < abs(dQb):
[EFa, EFb] = [EFb, EFa]
[dQa, dQb] = [dQb, dQa]
print(" EFmin = {:12.8f} dQmin = {:12.4g}".format(EFa, dQa))
print(" EFmax = {:12.8f} dQmax = {:12.4g}".format(EFb, dQb))
EFc = EFa
mflag = 1
for i in range(nmaxiter):
if abs(EFa - EFb) < eps:
EFh = (EFa + EFb) / 2.0
print(" Success: Convergence reached at EF = {}".format(EFh))
return 1
dQc = Ne(EFc, T) + NAm(EFc, T) - Nh(EFc, T) - NDp(EFc, T)
if abs(dQa - dQc) > 1.0e-10 and abs(dQb - dQc) > 1.0e-10:
#inverse quadratic interpolation
EFs = EFa * dQb * dQc / (dQa - dQb) / (dQa - dQc) \
+ EFb * dQa * dQc / (dQb - dQa) / (dQb - dQc) \
+ EFc * dQa * dQb / (dQc - dQa) / (dQc - dQb)
else:
#secant method
EFs = EFb - dQb * (EFb - EFa) / (dQb - dQa)
EF4 = (3.0*EFa + EFb) / 4.0
if not (EF4 < EFs < EFb or EFb < EFs < EF4):
mflag = 1
# print("condition 1")
elif mflag == 1 and abs(EFs - EFb) >= abs(EFb - EFc) / 2.0:
mflag = 1
# print("condition 2")
elif mflag == 0 and abs(EFs - EFb) >= abs(EFc - EFd) / 2.0:
mflag = 1
# print("condition 3")
elif mflag == 1 and abs(EFb - EFc) < delta:
mflag = 1
# print("condition 4")
elif mflag == 0 and abs(EFc - EFd) < delta:
mflag = 1
# print("condition 5")
else:
mflag = 0
# print("condition 6")
if mflag == 1:
EFs = (EFa + EFb) / 2.0 # bisection
dQs = Ne(EFs, T) + NAm(EFs, T) - Nh(EFs, T) - NDp(EFs, T)
EFd = EFc
dQd = dQc
EFc = EFb
# dQc = dQb
if dQa * dQs < 0.0:
# print("condition a: EFa,b,c,s=", EFa, EFb, EFc, EFs)
EFb = EFs
dQb = dQs
else:
# print("condition b: EFa,b,c,s=", EFa, EFb, EFc, EFs)
EFa = EFs
dQa = dQs
if abs(dQa) < abs(dQb):
[EFa, EFb] = [EFb, EFa]
[dQa, dQb] = [dQb, dQa]
Nes = Ne(EFs, T)
NAms = NAm(EFs, T)
Nhs = Nh(EFs, T)
NDps = NDp(EFs, T)
dQs = Nes + NAms - Nhs - NDps
if i % iprintinterval == 0:
print(" Iter {}: EFa,b = {:12.8f} - {:12.8f} dQa,b = {:12.4g} - {:12.8g}".format(
i, EFa, EFb, dQa, dQb))
print(" EFs = {:12.8f} dQs = {:12.4g}".format(i, EFs, dQs))
print(" Ne={:10.4e} Nh={:10.4e} NA-={:10.4e} ND+={:10.4e} dQ={:10.4e}".format(
Nes, Nhs, NAms, NDps, dQs))
else:
print(" Failed: Convergence did not reach")
return 0
if __name__ == "__main__":
main()