"""
半導体のキャリア統計量を自由電子近似 (FEA) でシミュレーションするスクリプト。
このスクリプトは、VASP計算結果に基づき、半導体の温度依存性、フェルミ準位依存性、
および有効質量に関する特性を自由電子近似を用いて計算・プロットします。
コマンドライン引数により、'me' (有効質量評価), 'T' (温度依存性), 'EF' (フェルミ準位依存性) の
3つのモードで動作します。
関連リンク: :doc:`EF-T-semi_FEA_usage`
"""
import os
import sys
from math import exp, sqrt, log
import numpy as np
#from numpy import arange
from scipy import optimize # newton関数はscipy.optimizeモジュールに入っている
from scipy.interpolate import interp1d
from matplotlib import pyplot as plt
from tklib.tkfile import tkFile
from tklib.tkutils import terminate, getarg, getintarg, getfloatarg, save_csv
#from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
from tklib.tkapplication import tkApplication
from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me
from tklib.tktransport.tkTransport import integrate_Simpson_list, fe, fh, NC2meff_FEA
from tklib.tktransport.tkDOS_FEA import integrate_Simpson, integrate_Simpson_list, tkDOS
[ドキュメント]
def pint(str):
"""
文字列を整数に変換します。
:param str: 変換する文字列。
:type str: str
:returns: 変換された整数。
:rtype: int
"""
return int(str)
[ドキュメント]
def pfloat(str):
"""
文字列を浮動小数点数に変換します。
:param str: 変換する文字列。
:type str: str
:returns: 変換された浮動小数点数。
:rtype: float
"""
return float(str)
#================================
# Global variables
#================================
app = tkApplication()
argv, narg = app.get_args()
params = app.get_param_dic()
#params.debug = 0
#mode: 'me', 'EF', 'T'
mode = 'me'
outcsvfile = None
# read from DOSCAR
dos = tkDOS()
dos.meeff = 0.3
dos.mheff = 1.0
dos.EV = 0.0
dos.EC = 1.1
dos.EA = 0.05 # Acceptor level, eV
dos.NA = 0.0e17
dos.ED = 1.05 # Donor level
dos.ND = 0.8e17
dos.EF0 = 0.0 # initial EF to find EF
T0 = 300.0
# For mode = 'me'. E range for fitting to me*(DOS)
dE_me_lsq = 0.5 # measured from EC/EV, eV
# EF range for mode = 'EF'
T0 = 300 # K
dEFmin = -0.5 # measured from EV, eV
dEFmax = 0.5 # measured from EC, eV
dos.nEF = 51
dos.Estep = 0.01 # Integration step, eV
# Temperature range for mode = 'T'
Tmin = 300 # K
Tmax = 600
nT = 11
Tstep = None
# Integration range w.r.t. kBT
dos.Einteg0 = -3.0
dos.Einteg1 = 3.0
dos.nrange = 6.0
# Bisection method parameters
# Initial search range
dos.dE_bisec = 0.4
# EPS
dos.eps_bisec = 1.0e-7
# max iteration number
dos.nmaxiter_bisec = 200
# Output cycle
dos.iprintiterval_bisec = 1
figsize = (6, 6)
#=============================
# Treat argments
#=============================
[ドキュメント]
def usage():
"""
スクリプトの利用方法とコマンドライン引数の例を表示します。
この関数は、プログラムが不正な引数で実行された場合や、ヘルプ表示が必要な場合に呼び出されます。
'me' (有効質量評価), 'T' (温度依存性), 'EF' (フェルミ準位依存性) の各モードにおける
引数の形式と例を詳細に説明します。
"""
global mode, dos
global T0
global dEFmin, dEFmax
global dE_me_lsq
print("")
print("Usage: Variables in () are optional")
print(" python {} mode (file args)".format(sys.argv[0], ))
print(" (i) python {} me dE".format(sys.argv[0]))
print(" Plot DOS and fit for me*(DOS)")
print(" dE: Energy range from EC/EV for fitting to m*(DOS)")
print(" ex: python {} me {}".format(sys.argv[0], dE_me_lsq))
print(" (ii) python {} T EC EA NA ED ND Tmin Tmax nT)".format(sys.argv[0]))
print(" Calculate T dependences of ne, nh, ND+, NA-")
print(" Tmin, Tmax, nT: T range and mesh")
print(" ex: python {} T {} {} {} {} {} {} {} {}".format(sys.argv[0], dos.EC, dos.EA, dos.NA, dos.ED, dos.ND, Tmin, Tmax, nT))
print(" (iii) python {} EF T EC EA NA ED ND nEF)".format(sys.argv[0]))
print(" Calculate EF dependences of ne, nh, ND+, NA-")
print(" nEF: # of EF mesh")
print(" ex: python {} EF {} {} {} {} {} {} {}".format(sys.argv[0], T0, dos.EC, dos.EA, dos.NA, dos.ED, dos.ND, dos.nEF))
[ドキュメント]
def updatevars():
"""
コマンドライン引数を解析し、グローバル変数を更新します。
スクリプト実行時に渡された引数を読み込み、選択された動作モード（'me', 'EF', 'T'）に基づいて、
対応するグローバル変数（例: バンド端エネルギー、ドーピング濃度、温度範囲など）を設定します。
引数が不正または不足している場合は、使用方法を表示してプログラムを終了します。
"""
global mode, dos, outcsvfile
global T0
global dEFmin, dEFmax
global Tmin, Tmax, nT, Tstep
global dE_me_lsq
argv = sys.argv
if len(argv) == 1:
usage()
exit()
mode = getarg( 1, mode)
if mode == 'me':
dE_me_lsq = getfloatarg(2, dE_me_lsq)
elif mode == 'EF':
T0 = getfloatarg(2, T0)
dos.EC = getfloatarg(3, dos.EC)
dos.EA = getfloatarg(4, dos.EA)
dos.NA = getfloatarg(5, dos.NA)
dos.ED = getfloatarg(6, dos.ED)
dos.ND = getfloatarg(7, dos.ND)
dos.nEF = getintarg(8, dos.nEF)
dos.Eg = dos.EC - dos.EV
elif mode == 'T':
dos.EC = getfloatarg(2, dos.EC)
dos.EA = getfloatarg(3, dos.EA)
dos.NA = getfloatarg(4, dos.NA)
dos.ED = getfloatarg(5, dos.ED)
dos.ND = getfloatarg(6, dos.ND)
dos.Eg = dos.EC - dos.EV
Tmin = getfloatarg(7, Tmin)
Tmax = getfloatarg(8, Tmax)
nT = getintarg(9, nT)
Tstep = (Tmax - Tmin) / (nT - 1)
outcsvfile = f'EF-T-semi_FEA-{mode}.csv'
[ドキュメント]
def exec_T():
"""
温度依存性のキャリア統計量を計算し、結果をプロットしてCSVに保存します。
指定された温度範囲 (`Tmin` から `Tmax`) で、フェルミ準位 (EF)、
電子濃度 (Ne)、正孔濃度 (Nh)、ドナーイオン濃度 (ND+)、アクセプタイオン濃度 (NA-)、
ホール係数 (RH)、正味キャリア濃度 (Ns)、および活性化エネルギー (Ea) を計算します。
計算されたデータは、標準出力に表示され、CSVファイルに保存されます。
また、Matplotlibを使用して、これらの値の温度依存性を示す複数のグラフが生成されます。
"""
global mode, dos, T0
print("")
print("T dependence")
print(" Temperature range: {} - {}, {} K step".format(Tmin, Tmax, Tstep))
xT = [Tmin + iT * Tstep for iT in range(nT)]
xTinv, yEF, yNe, yNh, yNAm, yNDp, yRH, yNs, yNsabs, ydQ = dos.cal_T_lists(xT, print_level = 1)
yEaNe = []
yEaNeTh = []
yEaNh = []
yEaNhTh = []
print("")
print ("%8s %10s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s"
% ("T(K)", "EF(eV)", "Ea(Ne,meV)", "Ea(Ne,T1/2)", "Ea(Nh,meV)", "Ea(Nh,T1/2)", "Ne", "Nh", "ND+", "NA-", "RH", "Ns", "dQ"))
for i in range(nT):
if i == 0:
i1 = i
i2 = i+1
elif i == nT - 1:
i1 = i-1
i2 = i
else:
i1 = i-1
i2 = i+1
Th1 = sqrt(xT[i1])
Th2 = sqrt(xT[i2])
dT = 1.0 / (xTinv[i2] - xTinv[i1])
eps = 1.0e-300
EaNe = (log(yNe[i2]+eps) - log(yNe[i1]+eps)) * dT # meV
EaNeTh = (log(yNe[i2]/Th2+eps) - log(yNe[i1]/Th1+eps)) * dT # meV
EaNh = (log(yNh[i2]+eps) - log(yNh[i1]+eps)) * dT # meV
EaNhTh = (log(yNh[i2]/Th2+eps) - log(yNh[i1]/Th1+eps)) * dT # meV
yEaNe.append (-1000.0 * 1000.0 * kB / e * EaNe)
yEaNeTh.append(-1000.0 * 1000.0 * kB / e * EaNeTh)
yEaNh.append (-1000.0 * 1000.0 * kB / e * EaNh)
yEaNhTh.append(-1000.0 * 1000.0 * kB / e * EaNhTh)
print("%8.3f %10.4f %12.4g %12.4g %12.4g %12.4g " % (xT[i], yEF[i], yEaNe[i], yEaNeTh[i], yEaNh[i], yEaNhTh[i])
+ "%12g %12g %12g %12g %12g %12g %12g" % (yNe[i], yNh[i], yNDp[i], yNAm[i], yRH[i], yNs[i], ydQ[i]))
print("")
print("Save simulation data to [{}]".format(outcsvfile))
save_csv(outcsvfile,
["T(K)", "1000/T(K^-1)", "EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH", "Ns", "Ea(Ne,meV)", "Ea(Ne,T1/2)", "Ea(Nh,meV)", "Ea(Nh,T1/2)", "dQ"],
[xT, xTinv, yEF, yNe, yNh, yNDp, yNAm, yRH, yNs, yEaNe, yEaNeTh, yEaNh, yEaNhTh, ydQ])
#=============================
# Plot graphs
#=============================
fig = plt.figure(figsize = figsize)
ax1 = fig.add_subplot(2, 3, 1)
ax2 = fig.add_subplot(2, 3, 2)
ax3 = fig.add_subplot(2, 3, 4)
ax4 = fig.add_subplot(2, 3, 3)
ax5 = fig.add_subplot(2, 3, 6)
ax1.plot(xT, yEF, label = '$E_F$', color = 'black')
ax1.set_xlabel("T (K)")
ax1.set_ylabel("$E_F$ (eV)")
ax1.legend()
ax2.plot(dos.Elist, dos.DOSlist, label = 'DOS', color = 'black')
ax2.plot([dos.EF0, dos.EF0], [min(dos.DOSlist), max(dos.DOSlist)], color = 'red', linestyle = 'dashed')
ax2.set_xlabel("E (eV)")
ax2.set_ylabel("DOS (states/cm$^3$)")
ax2.legend()
ax3.plot(xT, yRH, label = 'RH', color = 'black', marker = 'o')
ax3.set_xlabel("T (K)")
ax3.set_ylabel("R$_H$ (C$^{-1}$cm$^{-3}$)")
ax3.legend()
ax4.plot(xTinv, yNe, label = '$N_e$', color = 'black', marker = 'o', linewidth = 0.5, markersize = 5.0)
ax4.plot(xTinv, yNh, label = '$N_h$', color = 'red', marker = 'o', linewidth = 0.5, markersize = 5.0)
ax4.plot(xTinv, yNsabs, label = '|$N_s$|', color = 'purple', linewidth = 0.5, linestyle = 'dashed', marker = '^', markersize = 2.0)
ax4.plot(xTinv, yNDp, label = '$N_D^+$', color = 'blue', marker = 'x', linewidth = 0.5, markersize = 5.0)
ax4.plot(xTinv, yNAm, label = '$N_A^-$', color = 'green', marker = 'x', linewidth = 0.5, markersize = 5.0)
ax4.set_yscale('log')
minNs = [min(yNe), min(yNh), min(yNsabs)]
maxNs = [max(yNe), max(yNh), max(yNsabs)]
ylim = [ max([min(minNs)*0.5, 1.0e8]), min([max(maxNs)*2.0, 1.0e23]) ]
ax4.set_ylim(ylim)
# ax4.plot([0, 0], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax4.set_xlabel("1000/T (K$^{-1}$)")
ax4.set_ylabel("N (cm$^{-3}$)")
ax4.legend()
ax5.plot(xTinv, yEaNe, label = '$E_a(N_e)$', linewidth = 0.5, color = 'black', marker = 'o', markersize = 5.0)
ax5.plot(xTinv, yEaNeTh, label = '$E_a(N_e*T^{-0.5}$)', linewidth = 0.5, color = 'black', marker = 'x', markersize = 5.0)
ax5.plot(xTinv, yEaNh, label = '$E_a(N_h)$', linewidth = 0.5, color = 'red', marker = 'o', markersize = 5.0)
ax5.plot(xTinv, yEaNhTh, label = '$E_a(N_h*T^{-0.5}$)', linewidth = 0.5, color = 'red', marker = 'x', markersize = 5.0)
ax5.set_xlabel("1000/T (K$^{-1}$)")
ax5.set_ylabel("Ea (meV)")
ax5.legend()
# Rearange the graph axes so that they are not overlapped
plt.tight_layout()
plt.pause(0.1)
print("")
usage()
print("")
print("Press ENTER to exit>>", end = '')
input()
[ドキュメント]
def exec_EF():
"""
フェルミ準位依存性のキャリア統計量を計算し、結果をプロットしてCSVに保存します。
指定された温度 (`T0`) とフェルミ準位範囲 (`dos.EFmin` から `dos.EFmax`) で、
電子濃度 (Ne)、正孔濃度 (Nh)、ドナーイオン濃度 (ND+)、アクセプタイオン濃度 (NA-)、
ホール係数 (RH)、正味キャリア濃度 (Ns)、および特性温度 (T0_Ne, T0_Nh) を計算します。
計算されたデータは、標準出力に表示され、CSVファイルに保存されます。
また、Matplotlibを使用して、これらの値のフェルミ準位依存性を示す複数のグラフが生成されます。
バンドギャップ中央での有効状態密度 (NC, NV) と有効質量 (me*, mh*) も近似計算し表示します。
"""
global mode, dos, T0
print("")
print("EF dependence at {} K".format(T0))
print(" EF range: {} - {}, {:12.6g} eV step, nE={}".format(dos.EFmin, dos.EFmax, dos.EFstep, dos.nEF))
xEF = [dos.EFmin + i * dos.EFstep for i in range(dos.nEF)]
yNh, yNe, yNAm, yNDp, yRH, yNs, yNsabs = dos.cal_EF_lists(T0, xEF, print_level = 1)
yT0Ne = []
yT0Nh = []
eps = 1.0e-300
print("")
print ("%5s\t%8s\t%12s\t%12s\t%12s\t%12s\t%12s\t%12s\t%12s" % ("EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH", "T0(Ne)", "T0(Nh)", "Ns"))
for i in range(dos.nEF):
EF = xEF[i]
if i == 0:
i1 = i
i2 = i+1
elif i == dos.nEF - 1:
i1 = i-1
i2 = i
else:
i1 = i-1
i2 = i+1
dEF = xEF[i2] - xEF[i1]
dlnNedEF = (log(yNe[i2]+eps) - log(yNe[i1]+eps)) / dEF
dlnNhdEF = (log(yNh[i2]+eps) - log(yNh[i1]+eps)) / dEF
if dlnNedEF == 0.0:
dlnNedEF = 1.0-4
T0Ne = 1.0 / kB / dlnNedEF * e
if dlnNhdEF == 0.0:
dlnNhdEF = 1.0-4
T0Nh = -1.0 / kB / dlnNhdEF * e
yT0Ne.append(T0Ne)
yT0Nh.append(T0Nh)
print ("%10.4f\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g"
% (EF, yNe[i], yNh[i], yNDp[i], yNAm[i], yRH[i], T0Ne, T0Nh, yNs[i]))
print("")
Emid = (dos.EV + dos.EC) / 2.0
imid = int((Emid - dos.EFmin) / dos.EFstep + 1.0e-5)
Eimid = dos.EFmin + imid * dos.EFstep
print("Effective density of states at the mid gap {} eV (i={}, E={:8.4g}):".format(Emid, imid, Eimid))
Ne = yNe[imid]
Nh = yNh[imid]
T0Ne = yT0Ne[imid]
T0Nh = yT0Ne[imid]
NC = Ne * exp((dos.EC - Eimid) * e / kB / T0Ne)
NV = Nh * exp((Eimid - dos.EV) * e / kB / T0Nh)
meeff = NC2meff_FEA(NC, T0)
mheff = NC2meff_FEA(NV, T0)
print(" NC = {:12.6g} cm-3 (T0 = {:8.4g} K)".format(NC, T0Ne))
print(" NV = {:12.6g} cm-3 (T0 = {:8.4g} K)".format(NV, T0Nh))
print(" me* = {:8.4g} me".format(meeff))
print(" mh* = {:8.4g} me".format(mheff))
yNeapprox = []
yNhapprox = []
for i in range(dos.nEF):
EF = xEF[i]
kexpc = (dos.EC - EF) * e / kB / T0
kexpv = (EF - dos.EV) * e / kB / T0
yNeapprox.append(dos.NC * exp(-kexpc))
yNhapprox.append(dos.NV * exp(-kexpv))
print("")
print("Save simulation data to [{}]".format(outcsvfile))
save_csv(outcsvfile,
["EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH", "T0(Ne)", "T0(Nh)", "Ns"],
[xEF, yNe, yNh, yNDp, yNAm, yRH, yT0Ne, yT0Nh, yNs])
#=============================
# Plot graphs
#=============================
fig = plt.figure(figsize = figsize)
ax2 = fig.add_subplot(2, 2, 1)
ax3 = fig.add_subplot(2, 2, 2)
ax4 = fig.add_subplot(2, 2, 3)
ax1 = fig.add_subplot(2, 2, 4)
ax2.plot(dos.Elist, dos.DOSlist, label = 'DOS', color = 'black')
ax2.plot([dos.EF0, dos.EF0], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax2.plot([dos.EV, dos.EV], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax2.plot([dos.EC, dos.EC], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax2.set_xlabel("E (eV)")
ax2.set_ylabel("DOS (states/cm$^3$)")
ax2.legend()
ax3.plot(xEF, yNe, label = '$N_e$', color = 'black', linewidth = 0.5) #, marker = 'o')
ax3.plot(xEF, yNh, label = '$N_h$', color = 'red', linewidth = 0.5) #, marker = 'x')
ax3.plot(xEF, yNsabs, label = '|$N_s$|', color = 'purple', linewidth = 0.5, linestyle = 'dashed') #, marker = '^')
ax3.plot(xEF, yNDp, label = '$N_D^+$', color = 'blue') #, marker = 'x')
ax3.plot(xEF, yNAm, label = '$N_A^-$', color = 'green') #, marker = 'x')
ax3.plot(xEF, yNeapprox, label = 'Ne(Boltz)', color = 'gray', linewidth = 0.5) #, marker = 'x')
ax3.plot(xEF, yNhapprox, label = 'Nh(Boltz)', color = 'purple', linewidth = 0.5) #, marker = 'x')
ax3.plot([dos.EV, dos.EV], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax3.plot([dos.EC, dos.EC], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax3.set_yscale('log')
ax3.set_xlabel("EF (eV)")
ax3.set_ylabel("N (cm$^{-3}$)")
ax3.legend()
ax4.plot(xEF, yRH, label = 'R$_H$', color = 'black', linewidth = 0.5, marker = 'o')
ax4.set_xlabel("EF (eV)")
ax4.set_ylabel("R$_H$ (C$^{-1}$cm$^{-3}$)")
ax4.legend()
ax1.plot(xEF, yT0Ne, label = 'T$_0$(Ne)', color = 'black') #, marker = 'o')
ax1.plot(xEF, yT0Nh, label = 'T$_0$(Nh)', color = 'red') #, marker = 'x')
ax1.plot([dos.EV, dos.EV], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax1.plot([dos.EC, dos.EC], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax1.set_ylim([0.0, 2.0*T0])
ax1.set_xlabel("EF (eV)")
ax1.set_ylabel("T$_0$ (K)")
ax1.legend()
# Rearange the graph axes so that they are not overlapped
plt.tight_layout()
plt.pause(0.1)
print("")
usage()
print("")
print("Press ENTER to exit>>", end = '')
input()
[ドキュメント]
def exec_me():
"""
状態密度 (DOS) をプロットし、有効質量 (m*) を計算します。
DOSリストから計算された状態密度と、その二乗値をプロットします。
また、DOSのエネルギー微分を用いて有効質量を計算し、そのエネルギー依存性もプロットします。
有効質量の計算は、DOSが自由電子近似の式にフィットすると仮定して行われます。
計算されたデータはグラフとして表示され、ユーザーの入力を待って終了します。
"""
global dos
print("")
print("Plot DOS and calculate m*")
print(" dE_me_lsq: {:8.4g} eV".format(dE_me_lsq))
print(" EV: {:8.4g} eV dEFmin: {:8.4g} eV".format(dos.EV, dEFmin))
print(" EC: {:8.4g} eV dEFmax: {:8.4g} eV".format(dos.EC, dEFmax))
plotE0 = min([dos.EV + dEFmin, dos.EV - dE_me_lsq])
plotE1 = max([dos.EC + dEFmax, dos.EC + dE_me_lsq])
print("Plot E range: {:8.4g} - {:8.4g} eV".format(plotE0, plotE1))
ne_from_dos = integrate_Simpson_list(dos.Elist, dos.DOSlist)
E_me = []
dos2 = []
Emid = (dos.EV + dos.EC) / 2.0
c = 0
for i in range(dos.nE):
E = dos.Elist[i]
# print("ne=", dos.Emin, E, dos.DOS(E), ne)
if E < plotE0 or plotE1 < E:
continue
E_me.append(E)
dos2.append(dos.DOSlist[i]**2)
me_me = []
nme = len(E_me)
for i in range(nme):
if i == 0:
i0 = 0
i1 = 1
elif i == nme - 1:
i0 = i - 1
i1 = i
else:
i0 = i - 1
i1 = i + 1
DC0 = (dos2[i1] - dos2[i0]) / (E_me[i1] - E_me[i0])
if DC0 < 0.0:
sign = -1.0
else:
sign = 1.0
DC0 = sqrt(abs(DC0))
# DC0 = sqrt(2.0) / pi / pi * pow(me * meff, 1.5) / hbar / hbar / hbar * 1.0e-6 * pow(e, 1.5)
meff = pow(DC0 / sqrt(2.0) * pi * pi * hbar * hbar * hbar / 1.0e-6 / pow(e, 1.5), 2.0/3.0) / me
meff *= sign
# print("m=", DC0, meff)
me_me.append(meff)
#=============================
# Plot graphs
#=============================
fig = plt.figure(figsize = figsize)
ax1 = fig.add_subplot(2, 2, 1)
ax2 = fig.add_subplot(2, 2, 2)
ax3 = fig.add_subplot(2, 2, 3)
ax4 = fig.add_subplot(2, 2, 4)
ax1.plot(dos.Elist, dos.DOSlist, label = 'DOS', color = 'black', linewidth = 0.5)
ax1.plot([dos.EF0, dos.EF0], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax1.plot([dos.EV, dos.EV], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax1.plot([dos.EC, dos.EC], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax1.set_xlabel("E (eV)")
ax1.set_ylabel("DOS (states/cm$^3$/eV)")
ax1.legend()
# print("len=", len(dos.Elist), len(ne_from_dos))
# print("ne=", ne_from_dos)
ax2.plot(dos.Elist, ne_from_dos, label = 'N', color = 'black', linewidth = 0.5)
ylim = [min(ne_from_dos), max(ne_from_dos)]
ax2.plot([dos.EF0, dos.EF0], ylim, label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax2.plot([dos.EV, dos.EV], ylim, label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax2.plot([dos.EC, dos.EC], ylim, label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax2.set_xlabel("E (eV)")
ax2.set_ylabel("N (states/cm$^3$)")
ax2.legend()
ax3.plot(E_me, dos2, label = 'DOS$^2$', color = 'black', linewidth = 0.5)
ylim = [min(dos2), max(dos2)]
ax3.plot([dos.EF0, dos.EF0], ylim, label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax3.plot([dos.EV, dos.EV], ylim, label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax3.plot([dos.EC, dos.EC], ylim, label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
# ax3.plot([dos.EV - dE_me_lsq, dos.EV - dE_me_lsq], ylim, label = '$E_{fit}$', color = 'green', linestyle = 'dashed', linewidth = 0.5)
# ax3.plot([dos.EC + dE_me_lsq, dos.EC + dE_me_lsq], ylim, label = '$E_{fit}$', color = 'green', linestyle = 'dashed', linewidth = 0.5)
ax3.set_xlabel("E (eV)")
ax3.set_ylabel("DOS$^2$ (states$^2$/cm$^6$)")
ax3.legend()
ax4.plot(E_me, me_me, label = '$m^*$', color = 'black', linewidth = 0.5)
ylim = [min(me_me), max(me_me)]
ax4.plot([dos.EF0, dos.EF0], ylim, label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax4.plot([dos.EV, dos.EV], ylim, label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax4.plot([dos.EC, dos.EC], ylim, label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
# ax4.plot([dos.EV - dE_me_lsq, dos.EV - dE_me_lsq], ylim, label = '$E_{fit}$', color = 'green', linestyle = 'dashed', linewidth = 0.5)
# ax4.plot([dos.EC + dE_me_lsq, dos.EC + dE_me_lsq], ylim, label = '$E_{fit}$', color = 'green', linestyle = 'dashed', linewidth = 0.5)
ax4.set_xlabel("E (eV)")
ax4.set_ylabel("$m^*$")
ax4.legend()
# Rearange the graph axes so that they are not overlapped
plt.tight_layout()
plt.pause(0.1)
print("")
usage()
print("")
print("Press ENTER to exit>>", end = '')
input()
[ドキュメント]
def main():
"""
スクリプトのメイン実行関数です。
コマンドライン引数を解析し、必要なグローバル変数を初期化します。
半導体のバンド構造とドーピング濃度に基づいて、状態密度 (DOS) 計算のための
エネルギー範囲とメッシュを設定します。
その後、選択された動作モード ('me', 'T', 'EF') に応じて、
対応する実行関数 (`exec_me`, `exec_T`, `exec_EF`) を呼び出します。
不正なモードが指定された場合は、エラーメッセージを表示し、利用方法を出力して終了します。
"""
global mode, dos, T0
global Vcell, file, E_raw, dos_raw
global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
updatevars()
print("")
print("===============================================================================")
print(" Calculate T dependences of semiconductor statistics using VASP results")
print("===============================================================================")
print("mode: ", mode)
dos.Eg = dos.EC - dos.EV
dos.NC = dos.meff2NC_FEA(dos.meeff, T0)
dos.NV = dos.meff2NC_FEA(dos.mheff, T0)
dos.DC0 = dos.meff2DC0_FEA(dos.meeff, T0)
dos.DV0 = dos.meff2DC0_FEA(dos.mheff, T0)
print("")
print("Eg: {} eV".format(dos.Eg))
print("Valence band:")
print(" EV : {} eV".format(dos.EV))
print(" mh*: {} me".format(dos.mheff))
print(" NV : {} me at {} K".format(dos.NV, T0))
print(" DV0: {} me at {} K".format(dos.DV0, T0))
print("Conduction band:")
print(" EC : {} eV".format(dos.EC))
print(" me*: {} me".format(dos.meeff))
print(" NC : {} me at {} K".format(dos.NC, T0))
print(" DC0: {} me at {} K".format(dos.DC0, T0))
print("")
print("Acceptor states:")
print(" {} cm^-3 at {} eV".format(dos.NA, dos.EA))
print("Donoar states:")
print(" {} cm^-3 at {} eV".format(dos.ND, dos.ED))
dos.EFmin = dos.EV + dEFmin
dos.EFmax = dos.EC + dEFmax
dos.EFstep = (dos.EFmax - dos.EFmin) / (dos.nEF - 1)
dos.Emin = dos.EV + dEFmin
dos.Emax = dos.EC + dEFmax
dos.nE = int((dos.Emax - dos.Emin) / dos.Estep + 1.00001)
dos.Elist = [dos.Emin + i * dos.Estep for i in range(dos.nE)]
dos.DOSlist = [dos.DOS(dos.Elist[i]) for i in range(dos.nE)]
print("")
print("DOS range: {:8.4g} - {:8.4g}, {:8.4g} steps, nE={}".format(dos.Emin, dos.Emax, dos.Estep, dos.nE))
print("")
print("Integration configuration")
print(" Integration E range: {} - {} eV, or EF+-{}*kBT eV".format(dos.Einteg0, dos.Einteg1, dos.nrange))
if mode == 'me':
exec_me()
elif mode == 'T':
exec_T()
elif mode == 'EF':
exec_EF()
else:
terminate("Error: Invalid mode [{}]".format(mode), usage)
if __name__ == '__main__':
main()