cms.integration.EF_T_semiconductor のソースコード

import sys
import numpy as np
from numpy import sin, cos, tan, pi, exp
from matplotlib import pyplot as plt


"""
非縮退半導体のEFの温度依存性を二分法で計算。
"""


#定数
e  = 1.60218e-19      # C";
kB = 1.380658e-23     # JK<sup>-1</sup>";

# 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

# Temperature range
Tmin  =   50  # K
Tmax  = 1000
Tstep =   10
nT = int((Tmax - Tmin) / Tstep + 1)


# Treat argments
if __name__ == "__main__":
    argv = sys.argv
    if len(argv) <= 1:
        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()

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

#EFの誤差がepsより小さくなったら計算終了
eps   = 1.0e-5
#二分法の最大繰り返し数
nmaxiter = 200
#繰り返し中に途中経過を出力するサイクル数
iprintiterval = 1

# Fermi-Dirac function

[ドキュメント]
def fe(E, EF, T):
    global e, kB
    if T == 0.0:
        if E < EF:
            return 1.0
        else:
            return 0.0
    return 1.0 / (exp((E - EF) * e / kB / T) + 1.0)


# electron density

[ドキュメント]
def Ne(EF, T):
    global Nc, Ec, e, kB
    if T == 0.0:
        return 0.0
    return Nc * exp(-(Ec - EF) * e / kB / T)


# hole density

[ドキュメント]
def Nh(EF, T):
    global Nv, Ev, e, kB
    if T == 0.0:
        return 0.0
    return Nv * exp(-(EF - Ev) * e / kB / T)


# 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("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 = Ev - 1.0
        EFmax = Ec + 1.0
# まず、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))
                    break
            if dQmin * dQhalf < 0.0:
                EFmax = EFhalf
                dQmax = dQhalf
            else:
                EFmin = EFhalf
                dQmin = dQhalf
        else:
            print("  Failed: Convergence did not reach")
            return 0

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

    ax1 = fig.add_subplot(1, 2, 1)
    ax2 = fig.add_subplot(1, 2, 2)

    ax1.plot(xT, yEF, label = 'EF')
    ax1.plot([Tmin, Tmax], [Ev, Ev], color = 'red')
    ax1.plot([Tmin, Tmax], [Ec, Ec], color = 'red')
    ax1.set_xlabel("T (K)")
    ax1.set_ylabel("EF (eV)")
    ax1.set_ylim([Ev, Ec * 1.1])
    ax1.legend()

    ax2.plot(xInvT, yNe, label = 'Ne')
    ax2.plot(xInvT, yNh, label = 'Nh')
    ax2.plot(xInvT, yNAm, label = 'NA-')
    ax2.plot(xInvT, yNDp, label = 'ND+')
    ax2.set_yscale("log")    
    ax2.set_xlabel("1000/T (K^-1)")
    ax2.set_ylabel("N (cm^-3)")
    maxy = max([max(yNe), max(yNh), max(yNAm), max(yNDp)])
    ax2.set_ylim([1.0e8, maxy])
    ax2.legend()
    plt.subplots_adjust(wspace=0.4, hspace=0.6)

    plt.pause(0.1)
    

    print("Press ENTER to exit>>", end = '')
    input()




if __name__ == "__main__":
    main()