import csv
import numpy as np
from numpy import sin, cos, tan, pi, exp
import sys
"""
半導体のフェルミエネルギーを2分法で求めるスクリプト。
詳細説明:
本スクリプトは、半導体における電荷中性条件 (伝導帯電子密度 + イオン化アクセプター密度 - 価電子帯正孔密度 - イオン化ドナー密度 = 0)
を満たすフェルミエネルギー (EF) を2分法を用いて計算します。
指定された物理定数、半導体パラメータ、ドーピング濃度、温度のもとでEFを探索し、その過程と結果を出力します。
関連リンク: :doc:`EF-bisection_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
Ev = 0.0
#伝導帯下端エネルギー
Ec = 1.0
#価電子帯有効状態密度
Nv = 1.2e19
#伝導帯有効状態密度
Nc = 2.1e18
#アクセプター準位
EA = 0.05
#アクセプター密度 cm-3
NA = 1.0e15
#ドナー準位
ED = Ec - 0.05
#ドナー密度
ND = 1.0e16
# parameters
#初期範囲として、価電子帯上端と伝導帯下端エネルギーを設定する
EFmin = Ev
#EFがepsより小さくなったら敬さん終了
EFmax = Ec
eps = 1.0e-5
#2分法の最大繰り返し数
nmaxiter = 200
#繰り返し中に途中経過を出力するサイクル数
iprintiterval = 1
# Fermi-Dirac function
[ドキュメント]
def fe(E, EF, T):
"""
フェルミ・ディラック分布関数を計算する。
詳細説明:
フェルミ準位EFと温度Tにおけるエネルギー準位Eの電子占有確率を計算する。
グローバル変数 ekBT が計算に使用される。
:param E: float: エネルギー準位 (eV)。
:param EF: float: フェルミ準位 (eV)。
:param T: float: 温度 (K)。
:returns: float: 電子占有確率。
"""
global Nc, Ec, kB
return 1.0 / (exp((E - EF) * ekBT) + 1.0)
# electron density
[ドキュメント]
def Ne(EF, T):
"""
伝導帯の電子密度を計算する。
詳細説明:
フェルミ準位EFと温度Tにおける伝導帯の電子密度をボルツマン近似を用いて計算する。
グローバル変数 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):
"""
価電子帯の正孔密度を計算する。
詳細説明:
フェルミ準位EFと温度Tにおける価電子帯の正孔密度をボルツマン近似を用いて計算する。
グローバル変数 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):
"""
イオン化ドナー密度を計算する。
詳細説明:
フェルミ準位EFと温度Tにおけるイオン化されたドナーの密度を計算する。
ドナーは電子を放出して正にイオン化する。
グローバル変数 ND, ED, ekBT が計算に使用される。
: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):
"""
イオン化アクセプター密度を計算する。
詳細説明:
フェルミ準位EFと温度Tにおけるイオン化されたアクセプターの密度を計算する。
アクセプターは電子を受け取って負にイオン化する。
グローバル変数 NA, EA, ekBT が計算に使用される。
:param EF: float: フェルミ準位 (eV)。
:param T: float: 温度 (K)。
:returns: float: イオン化アクセプター密度 (cm^-3)。
"""
global NA, EA, kB
return NA * fe(EA, EF, T)
[ドキュメント]
def main():
"""
2分法を用いて半導体のフェルミエネルギーを計算し、結果を出力する。
詳細説明:
電荷中性条件 (Ne + NAm - Nh - NDp = 0) を満たすフェルミエネルギー (EF) を2分法で探索する。
初期範囲 (EFmin, EFmax) が適切であるか確認し、指定された収束条件 (eps)
または最大繰り返し回数 (nmaxiter) に達するまで計算を繰り返す。
計算の途中経過と最終結果を標準出力に表示する。
グローバル変数 EFmin, EFmax, eps, nmaxiter, iprintinterval, T が計算に使用される。
:returns: int: 計算が成功した場合は1、失敗した場合は0。
"""
global EFmin, EFmax, eps, nmaxiter, iprintiterval
print("Solution of EF by bisection method")
print("")
# まず、EFmin,EFmaxにおけるΔQを計算し、それぞれが正・負あるいは負・生となっていることを確認する
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
# 2分法開始
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))
# EFの精度がepsより小さくなったら敬さん終了
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()