2分法で半導体のフェルミエネルギーを求める

Download script from .\EF-bisection.py

import csv
import numpy as np
from numpy import sin, cos, tan, pi, exp
import sys


"""
2分法で半導体のフェルミエネルギーを求める
"""

# physical parameeters
e = 1.60218e-19 # C";
kB = 1.380658e-23 # JK-1";

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
EFmax = Ec
#EFがepsより小さくなったら敬さん終了
eps = 1.0e-5
#2分法の最大繰り返し数
nmaxiter = 200
#繰り返し中に途中経過を出力するサイクル数
iprintiterval = 1

# Fermi-Dirac function
def fe(E, EF, T):
    global Nc, Ec, kB
    return 1.0 / (exp((E - EF) * ekBT) + 1.0)

# electron density
def Ne(EF, T):
    global Nc, Ec, kB
    return Nc * exp(-(Ec - EF) * ekBT)

# hole density
def Nh(EF, T):
    global Nv, Ev, kB
    return Nv * exp(-(EF - Ev) * ekBT)

# ionized donor density
def NDp(EF, T):
    global ND, ED, kB
    return ND * (1.0 - fe(ED, EF, T))

# ionized acceptor density
def NAm(EF, T):
    global NA, EA, kB
    return NA * fe(EA, EF, T)


def main():
    global EFmin, EFmax, eps, nmaxiter, iprintinterval

    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()