"""
VASPのDOSCARファイルに基づき、キャリア密度とフェルミ準位を計算するスクリプト。
このスクリプトは、VASP計算で得られた状態密度 (DOS) データを用いて、半導体のキャリア輸送特性を評価します。
具体的には、電子密度、正孔密度、イオン化されたドナー・アクセプター密度、ホール係数などを、
温度の関数として、またはフェルミ準位の関数として計算します。
また、有効状態密度や特性温度などのパラメータも導出します。
関連リンク: :doc:`EF-T-DOS_usage`
"""
import os
import sys
from math import exp, sqrt, log
import numpy as np
from numpy import arange
from scipy import integrate # 数値積分関数 integrateを読み込む
from scipy import optimize # newton関数はscipy.optimizeモジュールに入っている
from scipy.interpolate import interp1d
from matplotlib import pyplot as plt
import csv # savecsv, read_csv関数で使用されるが、元のコードにはimportがなかったため追加。ルール1との兼ね合いで注意。ただし、これらの関数は現在未使用。
# constants
pi = 3.14159265358979323846
h = 6.6260755e-34 # Js";
hbar = 1.05459e-34 # "Js";
c = 2.99792458e8 # m/s";
e = 1.60218e-19 # C";
kB = 1.380658e-23 # JK<sup>-1</sup>";
me = 9.1093897e-31 # kg";
Debug = 0
#mode: 'EF', 'T'
mode = 'T'
# EF range for mode = 'EF'
T0 = 300 # K
dEFmin = -0.2 # measured from EV, eV
dEFmax = 0.2 # measured from EC, eV
nEF = 50
EFmin = None
EFmax = None
EFstep = None
# Temperature range for mode = 'T'
Tmin = 300 # K
Tmax = 600
nT = 11
Tstep = None
# Sample information (DOS.data is converted to states/eV/Vcell)
Vcell = 212.309 # A^3
# read from DOSCAR
file = "DOS.dat"
E_raw = None
dos_raw = None
nDOS = None
Emin = None
Emax = None
Estep = None
EV = None
EC = None
dEA = 0.32 # Acceptor level measured from EV, eV
NA = 0.8e17
dED = -0.28 # Donor level measured from EC, eV
ND = 0.5e17
EA = None # Acceptor states are not considered if EA = None
ED = None # Donor states are not considered if ED = None
# other parameters
EF0 = 0.4 # EF at charge neutral
Egth = 0.01 # Critera DOS value to search band edges
# integration parameters
# Integration range w.r.t. kBT
Einteg0 = -3.0
Einteg1 = 3.0
nrange = 6.0
# Relative accuracy of the quad functin
eps = 1.0e-5
nmaxdiv = 150
# Bisection method parameters
# Initial search range
dEbisec = 0.4
# EPS
epsbisec = 1.0e-7
# max iteration number
nmaxiter = 200
# Output cycle
iprintiterval = 1
# graph configuration
fontsize = 16
legend_fontsize = 12
#=============================
# Treat argments
#=============================
[ドキュメント]
def pfloat(str):
"""
文字列を浮動小数点数に変換する。
変換できない場合はNoneを返す。
:param str (str): 変換する文字列。
:returns: float または None: 変換された浮動小数点数、またはNone。
"""
try:
return float(str)
except:
return None
[ドキュメント]
def pint(str):
"""
文字列を整数に変換する。
変換できない場合はNoneを返す。
:param str (str): 変換する文字列。
:returns: int または None: 変換された整数、またはNone。
"""
try:
return int(str)
except:
return None
[ドキュメント]
def getarg(position, defval = None):
"""
コマンドライン引数を取得する。
指定された位置の引数が存在しない場合、デフォルト値を返す。
:param position (int): 取得する引数の位置（0はスクリプト名）。
:param defval (Any, optional): 引数が存在しない場合に返すデフォルト値。デフォルトはNone。
:returns: str または Any: 取得した引数、またはデフォルト値。
"""
try:
return sys.argv[position]
except:
return defval
[ドキュメント]
def getfloatarg(position, defval = None):
"""
コマンドライン引数を浮動小数点数として取得する。
`getarg`で取得した文字列を`pfloat`で浮動小数点数に変換する。
:param position (int): 取得する引数の位置。
:param defval (float, optional): 引数が存在しない場合に返すデフォルト値。デフォルトはNone。
:returns: float または None: 変換された浮動小数点数、またはNone。
"""
return pfloat(getarg(position, defval))
[ドキュメント]
def getintarg(position, defval = None):
"""
コマンドライン引数を整数として取得する。
`getarg`で取得した文字列を`pint`で整数に変換する。
:param position (int): 取得する引数の位置。
:param defval (int, optional): 引数が存在しない場合に返すデフォルト値。デフォルトはNone。
:returns: int または None: 変換された整数、またはNone。
"""
return pint(getarg(position, defval))
[ドキュメント]
def usage():
"""
スクリプトの正しい使用方法を標準出力に表示する。
コマンドライン引数の形式と例を示す。
"""
global mode
global Vcell, T0, EF0
global file, E_raw, dos_raw, nDOS, Emin, Emax, Estep
global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
print("")
print("Usage: Variables in () are optional")
print(" python {} mode (file args)".format(sys.argv[0], ))
print(" (i) python {} T file Tmin Tmax nT)".format(sys.argv[0]))
print(" ex: {} T {} {} {} {}".format(sys.argv[0], file, Tmin, Tmax, nT))
print(" (ii) python {} EF file nEF)".format(sys.argv[0]))
print(" ex: {} EF {} {}".format(sys.argv[0], file, nEF))
[ドキュメント]
def updatevars():
"""
コマンドライン引数に基づいてグローバル変数を更新する。
`sys.argv`を解析し、`mode`, `file`, `Tmin`, `Tmax`, `nT`, `nEF`などの
グローバル変数を設定する。引数が不十分な場合は`usage`を呼び出して終了する。
"""
global mode
global Vcell, T0, EF0
global file, E_raw, dos_raw, nDOS, Emin, Emax, Estep
global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
global Tmin, Tmax, nT, Tstep
argv = sys.argv
if len(argv) == 1:
usage()
exit()
mode = getarg( 1, mode)
file = getarg( 2, file)
if mode == 'EF':
nEF = getintarg(3, nEF)
if mode == 'T':
Tmin = getintarg(3, Tmin)
Tmax = getintarg(4, Tmax)
nT = getintarg(5, nT)
Tstep = (Tmax - Tmin) / (nT - 1)
#=============================
# other functions
#=============================
[ドキュメント]
def savecsv(outfile, header, datalist):
"""
データをCSVファイルに保存する。
ヘッダーとデータリストを行ごとに書き出す。
:param outfile (str): 出力ファイル名。
:param header (list[str]): CSVファイルのヘッダー行。
:param datalist (list[list[Any]]): 列ごとのデータを含むリストのリスト。
:returns: None
"""
try:
print("Write to [{}]".format(outfile))
f = open(outfile, 'w')
except:
# except IOError:
print("Error: Can not write to [{}]".format(outfile))
else:
fout = csv.writer(f, lineterminator='\n')
fout.writerow(header)
# fout.writerows(data)
for i in range(0, len(datalist[0])):
a = []
for j in range(len(datalist)):
a.append(datalist[j][i])
fout.writerow(a)
f.close()
[ドキュメント]
def read_csv(infile, xmin = None, xmax = None, delimiter = ','):
"""
CSVファイルからデータを読み込む。
指定されたファイルからヘッダーと2列の数値データを読み込み、
xminとxmaxの範囲でフィルタリングする。
:param infile (str): 入力ファイル名。
:param xmin (float, optional): x値の下限（Noneの場合はフィルタリングなし）。デフォルトはNone。
:param xmax (float, optional): x値の上限（Noneの場合はフィルタリングなし）。デフォルトはNone。
:param delimiter (str, optional): 列の区切り文字。デフォルトは','。
:returns: tuple[list[str], list[float], list[float]]: ヘッダー、xデータ、yデータのタプル。
"""
print("xrange=", xmin, xmax)
data = []
try:
infp = open(infile, "r")
f = csv.reader(infp, delimiter = delimiter)
header = next(f)
print("header=", header)
for j in range(len(header)):
data.append([])
for row in f:
x = pfloat(row[0])
if xmin is not None and xmin <= x <= xmax:
y = pfloat(row[1])
data[0].append(x)
data[1].append(y)
except:
print("Error: Can not read [{}]".format(infile))
exit()
return header, data[0], data[1]
[ドキュメント]
def FindBandEdges(E, DOS, EF0, Egth):
"""
DOSデータから価電子帯上端 (EV) と伝導帯下端 (EC) を見つける。
フェルミ準位の初期値EF0を基準に、DOS値がEgthを超える範囲を探索してEVとECを特定する。
:param E (numpy.ndarray): エネルギー値の配列。
:param DOS (numpy.ndarray): DOS値の配列。
:param EF0 (float): フェルミ準位の初期値 (eV)。
:param Egth (float): バンド端を決定するためのDOSのスレッショルド値。
:returns: tuple[float, float]: 価電子帯上端EVと伝導帯下端ECのタプル。
"""
EV = 1.0e30
EC = 1.0e30
i0 = -1
for i in range(len(E)):
if E[i] >= EF0:
i0 = i
break
for i in range(i0, 0, -1):
if DOS[i] > Egth:
EV = E[i+1]
break
for i in range(i0, len(E), +1):
if DOS[i] > Egth:
EC = E[i-1]
break
return EV, EC
[ドキュメント]
def rieman(x0, dx, y, xmin, xmax):
"""
リーマン和を用いて数値積分を行う。
xminからxmaxの範囲で、等間隔のデータyに対してリーマン和を計算する。
:param x0 (float): データの開始x値。
:param dx (float): x軸のステップ幅。
:param y (list[float]): 積分対象のy値のリスト。
:param xmin (float): 積分範囲の下限。
:param xmax (float): 積分範囲の上限。
:returns: float: 計算された積分値。
"""
S = 0.0
for i in range(len(y)):
x = x0 + i * dx
if x < xmin or xmax < x:
continue
S += y[i]
return S * dx
# define the DOS function
[ドキュメント]
def DOS(E):
"""
指定されたエネルギーにおける状態密度 (DOS) を補間して返す。
グローバル変数`E_raw`と`dos_raw`を用いて、`scipy.interpolate.interp1d`による3次スプライン補間を行う。
DOSのE範囲外が指定された場合はエラーで終了する。
:param E (float): エネルギー (eV)。
:returns: float: 補間されたDOS値 (states/eV/Vcell単位をcm^-3に変換済み)。
"""
global E_raw, dos_raw
global Emin, Emax, Estep, nDOS
if E < Emin or Emax < E:
print("Error in f_dos: Given E={} exceeds the DOS E range [{}, {}]".format(E, Emin, Emax))
exit()
fdos = interp1d(E_raw, dos_raw, kind = 'cubic')
return fdos(E)
[ドキュメント]
def fe(E, T, EF):
"""
フェルミ・ディラック分布関数を計算する。
エネルギーE、温度T、フェルミ準位EFにおける電子の占有確率を返す。
T=0Kの場合はステップ関数として処理される。
:param E (float): エネルギー (eV)。
:param T (float): 温度 (K)。
:param EF (float): フェルミ準位 (eV)。
:returns: float: フェルミ・ディラック分布関数の値。
"""
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)
[ドキュメント]
def fh(E, T, EF):
"""
正孔の占有確率（1 - fe）を計算する。
フェルミ・ディラック分布関数`fe`を用いて正孔の占有確率を計算する。
:param E (float): エネルギー (eV)。
:param T (float): 温度 (K)。
:param EF (float): フェルミ準位 (eV)。
:returns: float: 正孔の占有確率。
"""
return 1.0 - fe(E, T, EF)
[ドキュメント]
def DOSfe(E, T, EF):
"""
DOSとフェルミ・ディラック分布関数の積を計算する。
`DOS(E)`と`fe(E, T, EF)`の積を返す。これは電子密度計算の被積分関数となる。
:param E (float): エネルギー (eV)。
:param T (float): 温度 (K)。
:param EF (float): フェルミ準位 (eV)。
:returns: float: `DOS(E) * fe(E, T, EF)` の値。
"""
return DOS(E) * fe(E, T, EF)
[ドキュメント]
def DOSfh(E, T, EF):
"""
DOSと正孔の占有確率の積を計算する。
`DOS(E)`と`fh(E, T, EF)`の積を返す。これは正孔密度計算の被積分関数となる。
:param E (float): エネルギー (eV)。
:param T (float): 温度 (K)。
:param EF (float): フェルミ準位 (eV)。
:returns: float: `DOS(E) * fh(E, T, EF)` の値。
"""
return DOS(E) * fh(E, T, EF)
[ドキュメント]
def integrate(func, E0, E1, h):
"""
台形公式を用いて数値積分を行う。
指定された関数`func`を`E0`から`E1`まで、ステップ幅`h`で台形公式により数値積分する。
:param func (callable): 積分対象の関数。
:param E0 (float): 積分範囲の下限。
:param E1 (float): 積分範囲の上限。
:param h (float): 積分ステップ幅。
:returns: list[float, float]: 積分結果とダミーのエラー値 (-1.0) のリスト。
"""
n = int((E1 - E0) / h + 1.000001)
h = (E1 - E0) / (n - 1)
y = [func(E0 + i * h) for i in range(n)]
S = 0.5 * (y[0] + y[n-1]) + sum(y[1:n-1])
return [h * S, -1.0]
[ドキュメント]
def Ne(T, EF, E0, E1):
"""
指定された範囲における電子密度を計算する。
`DOSfe`関数を`E0`から`E1`まで積分し、電子数を求める。
:param T (float): 温度 (K)。
:param EF (float): フェルミ準位 (eV)。
:param E0 (float): 積分範囲の下限。
:param E1 (float): 積分範囲の上限。
:returns: float: 計算された電子密度。
"""
global Estep
# N = integrate.quad(lambda E: DOSfe(E, T, EF), E0, E1, limit = nmaxdiv, epsrel = eps)
N = integrate(lambda E: DOSfe(E, T, EF), E0, E1, Estep)
if Debug:
print(" Ne integ error: {}".format(N[1]))
return N[0]
[ドキュメント]
def Nh(T, EF, E0, E1):
"""
指定された範囲における正孔密度を計算する。
`DOSfh`関数を`E0`から`E1`まで積分し、正孔数を求める。
:param T (float): 温度 (K)。
:param EF (float): フェルミ準位 (eV)。
:param E0 (float): 積分範囲の下限。
:param E1 (float): 積分範囲の上限。
:returns: float: 計算された正孔密度。
"""
global Estep
# N = integrate.quad(lambda E: DOSfh(E, T, EF), E0, E1, limit = nmaxdiv, epsrel = eps)
N = integrate(lambda E: DOSfh(E, T, EF), E0, E1, Estep)
if Debug:
print(" Nh integ error: {}".format(N[1]))
return N[0]
[ドキュメント]
def NDp(EF, T):
"""
イオン化されたドナー密度を計算する。
ドナー準位`ED`とドナー濃度`ND`、フェルミ・ディラック分布を用いてイオン化率を計算する。
:param EF (float): フェルミ準位 (eV)。
:param T (float): 温度 (K)。
:returns: float: イオン化されたドナー密度。
"""
global ND, ED, kB
return ND * (1.0 - fe(ED, T, EF))
[ドキュメント]
def NAm(EF, T):
"""
イオン化されたアクセプター密度を計算する。
アクセプター準位`EA`とアクセプター濃度`NA`、フェルミ・ディラック分布を用いてイオン化率を計算する。
:param EF (float): フェルミ準位 (eV)。
:param T (float): 温度 (K)。
:returns: float: イオン化されたアクセプター密度。
"""
global NA, EA, kB
n = NA * fe(EA, T, EF)
# print("n=", EF, T, NA, EA, n)
return n
[ドキュメント]
def CalEF(T, EF0, E0, E1, totNe):
"""
電荷中性条件に基づいてフェルミ準位を計算する (現在未使用の関数)。
`scipy.optimize.newton`メソッドを使用して、電子密度`Ne`が`totNe`に等しくなるフェルミ準位`EF`を見つける。
:param T (float): 温度 (K)。
:param EF0 (float): フェルミ準位の初期推測値 (eV)。
:param E0 (float): 積分範囲の下限。
:param E1 (float): 積分範囲の上限。
:param totNe (float): 目標とする全電子密度。
:returns: float: 電荷中性条件を満たすフェルミ準位。
"""
EF = optimize.newton(lambda EF: Ne(T, EF, E0, E1) - totNe, EF0)
return EF
[ドキュメント]
def exec_T():
"""
温度の関数としてキャリア密度とフェルミ準位を計算・プロットする。
指定された温度範囲でフェルミ準位を二分法で繰り返し決定し、
電子密度、正孔密度、イオン化ドナー/アクセプター密度、ホール係数などを計算する。
結果をプロットし、特性温度も算出する。
:returns: int または None: エラー発生時に0を返す。
"""
global mode, T0, EF0
global file, E_raw, dos_raw
global nDOS, Emin, Emax, Estep
global EV, EC
global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
global dEA, NA, dED, ND, EA, ED
xT = []
xTinv = []
yEF = []
yEint = []
yEentropy = []
yEFa = []
yNe = []
yNh = []
yNDp = []
yNAm = []
yRH = []
EFprev = EF0
#初期範囲
EFmin = EF0 - dEbisec
EFmax = EF0 + dEbisec
dEFbisec = 0.1
print ("%5s\t%8s\t%14s" % ("T(K)", "EF(eV)", "Ne(0 K) (check)"))
for iT in range(nT):
T = Tmin + iT * Tstep
print("")
print("iT {:03d}: {} K".format(iT, T))
kBTe = kB * T / e
dE = nrange * kBTe
E0 = EV - dE
E1 = EC + dE
if E0 < Emin or Emax < E1:
print("*** Integration range error: [{}, {}] exceeds DOS E range [{}, {}]".format(E0, E1, Emin, Emax))
exit()
Ne0K = Ne(0.0, EF0, E0, EF0)
# まず、EFmin,EFmaxにおけるΔQを計算し、それぞれが正・負あるいは負・生となっていることを確認する
print("E range:", EF0, E1)
# Nemin = Ne(T, EFmin, EF0, E1)
Nemin = Ne(T, EFmin, EC - Estep, E1)
# Nhmin = Nh(T, EFmin, E0, EF0)
Nhmin = Nh(T, EFmin, E0, EV + Estep)
NAmmin = NAm(EFmin, T)
NDpmin = NDp(EFmin, T)
dQmin = Nemin + NAmmin - Nhmin - NDpmin
# Nemax = Ne(T, EFmax, EF0, E1)
Nemax = Ne(T, EFmax, EC - Estep, E1)
# Nhmax = Nh(T, EFmax, E0, EF0)
Nhmax = Nh(T, EFmax, E0, EV + Estep)
NAmmax = NAm(EFmax, T)
NDpmax = NDp(EFmax, T)
dQmax = Nemax + NAmmax - Nhmax - NDpmax
print(" min (EF,Ne,Nh,NDp,NAm,dQ)=({:6.4f},{:10.4g},{:10.4g},{:10.4g},{:10.4g},{:10.4g})"
.format(EFmin, Nemin, Nhmin, NDpmin, NAmmin, dQmin))
print(" max (EF,Ne,Nh,NDp,NAm,dQ)=({:6.4f},{:10.4g},{:10.4g},{:10.4g},{:10.4g},{:10.4g})"
.format(EFmax, Nemax, Nhmax, NDpmax, NAmmax, dQmax))
# 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))
print(" dQ=", dQmin, dQmax)
print(" dQmin*dQmax=", dQmin * 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(T, EFhalf, EF0, E1)
Neh = Ne(T, EFhalf, EC - Estep, E1)
# Nhh = Nh(T, EFhalf, E0, EF0)
Nhh = Nh(T, EFhalf, E0, EV + Estep)
NAmh = NAm(EFhalf, T)
NDph = NDp(EFhalf, T)
dQhalf = Neh + NAmh - Nhh - NDph
print(" iter {:03d}: half (EF,Ne,Nh,NDp,NAm,dQ)=({:6.4f},{:10.4g},{:10.4g},{:10.4g},{:10.4g},{:10.4g})"
.format(i, EFhalf, Neh, Nhh, NDph, NAmh, dQhalf))
# EFの精度がepsより小さくなったら計算終了
if abs(EFmin - EFhalf) < epsbisec and abs(EFmax - EFhalf) < epsbisec:
# 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
EF = EFhalf
RH = 1.0 / e * (Nhh - Neh) / pow(Nhh + Neh, 2.0)
xT.append(T)
xTinv.append(1000.0 / T)
yEF.append(EF)
yNe.append(Neh)
yNh.append(Nhh)
yNDp.append(NDph)
yNAm.append(NAmh)
yRH.append(RH)
print(" ***(T,EF,Ne,Nh,NDp,NAm,RH)=({:6.4f},{:6.4g},{:10.4g},{:10.4g},{:10.4g},{:10.4g},{:10.4g})"
.format(T, EF, Neh, Nhh, NDph, NAmh, RH))
EFprev = EF
EFmin = EF - dEFbisec
EFmax = EF + dEFbisec
print ("%5.0f\t%12.6g" % (T, EF))
print("")
yEaNe = []
yEaNeTh = []
yEaNh = []
yEaNhTh = []
print ("%8s %8s %8s %8s %8s %8s %12s %10s %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"))
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])
EaNe = (log(yNe[i2]) - log(yNe[i1])) * dT # meV
EaNeTh = (log(yNe[i2]/Th2) - log(yNe[i1]/Th1)) * dT # meV
EaNh = (log(yNh[i2]) - log(yNh[i1])) * dT # meV
EaNhTh = (log(yNh[i2]/Th2) - log(yNh[i1]/Th1)) * 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("%5.3f %8.3f %8.3g %8.3g %8.3g %8.3g " % (xT[i], yEF[i], yEaNe[i], yEaNeTh[i], yEaNh[i], yEaNhTh[i])
+ "%12g %12g %12g %12g %12g" % (yNe[i], yNh[i], yNDp[i], yNAm[i], yRH[i]))
#=============================
# Plot graphs
#=============================
fig = plt.figure(figsize = (12, 8))
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 = 'EF', color = 'black')
ax1.set_xlabel("T (K)", fontsize = fontsize)
ax1.set_ylabel("E$_F$ (eV)", fontsize = fontsize)
ax1.legend(fontsize = legend_fontsize)
ax2.plot(E_raw, dos_raw, label = 'DOS', color = 'black')
ax2.plot([EF0, EF0], [min(dos_raw), max(dos_raw)], color = 'red', linestyle = 'dashed')
ax2.set_xlabel("E (eV)", fontsize = fontsize)
ax2.set_ylabel("DOS (states/cm$^3$)", fontsize = fontsize)
ax2.legend(fontsize = legend_fontsize)
ax3.plot(xT, yRH, label = 'RH', color = 'black', marker = 'o')
ax3.set_xlabel("T (K)", fontsize = fontsize)
ax3.set_ylabel("R$_H$ (C$^{-1}$cm$^{-3}$)", fontsize = fontsize)
ax3.legend(fontsize = legend_fontsize)
ax4.plot(xTinv, yNe, label = 'Ne', color = 'black', marker = 'o')
ax4.plot(xTinv, yNh, label = 'Nh', color = 'red', marker = 'o')
ax4.plot(xTinv, yNDp, label = 'ND$^+$', color = 'blue', marker = 'x')
ax4.plot(xTinv, yNAm, label = 'NA$^-$', color = 'green', marker = 'x')
ax4.set_yscale('log')
# ax4.plot([0, 0], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax4.set_xlabel("1000/T (K$^{-1}$)", fontsize = fontsize)
ax4.set_ylabel("N (cm$^{-3}$)", fontsize = fontsize)
ax4.legend(fontsize = legend_fontsize)
ax5.plot(xTinv, yEaNe, label = 'Ea(Ne)', color = 'black', marker = 'o')
ax5.plot(xTinv, yEaNeTh, label = 'Ea(Ne*T$^{-0.5}$)', color = 'black', marker = 'o')
ax5.plot(xTinv, yEaNh, label = 'Ea(Nh)', color = 'red', marker = 'o')
ax5.plot(xTinv, yEaNhTh, label = 'Ea(Nh*T$^{-0.5})', color = 'red', marker = 'o')
ax5.set_xlabel("1000/T (K$^{-1}$)", fontsize = fontsize)
ax5.set_ylabel("Ea (meV)", fontsize = fontsize)
ax5.legend(fontsize = legend_fontsize)
ax1.tick_params(labelsize = fontsize)
ax2.tick_params(labelsize = fontsize)
ax3.tick_params(labelsize = fontsize)
ax4.tick_params(labelsize = fontsize)
ax5.tick_params(labelsize = fontsize)
# Rearange the graph axes so that they are not overlapped
plt.tight_layout()
# plt.pause(0.1)
print("Close the graph window to terminate this program")
plt.show()
print("")
usage()
print("")
# print("Press ENTER to exit>>", end = '')
# input()
[ドキュメント]
def exec_EF():
"""
フェルミ準位の関数としてキャリア密度などを計算・プロットする。
指定されたフェルミ準位範囲で電子密度、正孔密度、イオン化ドナー/アクセプター密度、
ホール係数などを計算する。結果をプロットし、有効状態密度や特性温度も算出する。
:returns: None
"""
global mode, T0, EF0
global file, E_raw, dos_raw
global nDOS, Emin, Emax, Estep
global EV, EC
global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
global dEA, NA, dED, ND, EA, ED
kBTe = kB * T0 / e
dE = nrange * kBTe
E0 = EV - dE
E1 = EC + dE
xEF = []
yNe = []
yNh = []
yNDp = []
yNAm = []
yRH = []
EFprev = EF0
print("at {} K".format(T0))
print ("%5s\t%8s\t%14s\t%14s\t%14s\t%14s" % ("EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH"))
for iEF in range(nEF):
EF = EFmin + iEF * EFstep
dE = nrange * kBTe
# Ne1 = Ne(T0, EF, EF0, E1)
Ne1 = Ne(T0, EF, EC - Estep, E1)
# Nh1 = Nh(T0, EF, E0, EF0)
Nh1 = Nh(T0, EF, E0, EV + Estep)
NDp1 = NDp(EF, T0)
NAm1 = NAm(EF, T0)
RH = 1.0 / e * (Nh1 - Ne1) / pow(Nh1 + Ne1, 2.0)
xEF.append(EF)
yNe.append(Ne1)
yNh.append(Nh1)
yNDp.append(NDp1)
yNAm.append(NAm1)
yRH.append(RH)
print ("%10.4f\t%12.6g\t%12.6g\t%12.6g\t%12.6g\t%12.6g" % (EF, Ne1, Nh1, NDp1, NAm1, RH))
yT0Ne = []
yT0Nh = []
print("")
print ("%5s\t%8s\t%12s\t%12s\t%12s\t%12s\t%12s\t%12s" % ("EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH", "T0(Ne)", "T0(Nh)"))
for i in range(nEF):
if i == 0:
i1 = i
i2 = i+1
elif i == nEF - 1:
i1 = i-1
i2 = i
else:
i1 = i-1
i2 = i+1
dEF = xEF[i2] - xEF[i1]
dlnNedEF = (log(yNe[i2]) - log(yNe[i1])) / dEF
dlnNhdEF = (log(yNh[i2]) - log(yNh[i1])) / dEF
T0Ne = 1.0 / kB / dlnNedEF * e
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" % (EF, Ne1, Nh1, NDp1, NAm1, RH, T0Ne, T0Nh))
print("")
Emid = (EV + EC) / 2.0
print("Effective density of states at the mid gap {} eV:".format(Emid))
fNe = interp1d(xEF, yNe, kind = 'cubic')
fNh = interp1d(xEF, yNh, kind = 'cubic')
dlnNedEF = (log(fNe(Emid + Estep)) - log(fNe(Emid - Estep))) / 2.0 / Estep
dlnNhdEF = (log(fNh(Emid + Estep)) - log(fNh(Emid - Estep))) / 2.0 / Estep
T0Ne = 1.0 / kB / dlnNedEF * e
T0Nh = -1.0 / kB / dlnNhdEF * e
NC = fNe(Emid) * exp((EC - Emid) * e / kB / T0Ne)
NV = fNh(Emid) * exp((Emid - EV) * e / kB / T0Nh)
print(" NC = {} cm-3 (T0 = {} K)".format(NC, T0Nh))
print(" NV = {} cm-3 (T0 = {} K)".format(NV, T0Nh))
xEFapprox = [EFmin + i * EFstep for i in range(nEF)]
yNeapprox = [NC * exp(-(EC - xEFapprox[i]) * e / kB / T0) for i in range(nEF)]
yNhapprox = [NV * exp(-(xEFapprox[i] - EV) * e / kB / T0) for i in range(nEF)]
# calculate DC(E)*fe(E), DV(E)*(1-fe(E))
print("")
print("Calculate fe(E), fh(E)=1-fe(E), DC(E)*fe(E), DV(E)*fh(E) at EF={} eV, T={} K".format(EF0, T0))
yfe = []
yfh = []
yDOSfe = []
yDOSfh = []
for i in range(len(E_raw)):
E = E_raw[i]
yfe.append(fe(E, T0, EF0))
yfh.append(fh(E, T0, EF0))
if E >= EC:
yDOSfe.append(DOSfe(E, T0, EF0))
else:
yDOSfe.append(0.0)
if E <= EV:
yDOSfh.append(DOSfh(E, T0, EF0))
else:
yDOSfh.append(0.0)
#=============================
# Plot graphs
#=============================
fig = plt.figure(figsize = (12, 8))
ax2 = fig.add_subplot(2, 2, 1)
ax2r = ax2.twinx()
ax3 = fig.add_subplot(2, 2, 2)
ax4 = fig.add_subplot(2, 2, 3)
ax1 = fig.add_subplot(2, 2, 4)
ax2.plot(E_raw, dos_raw, label = 'DOS', color = 'black')
ax2.plot(E_raw, yDOSfe, label = 'D(E)fe', color = 'blue', linewidth = 0.5)
ax2.plot(E_raw, yDOSfh, label = 'D(E)fh', color = 'red', linewidth = 0.5)
ax2.plot([EF0, EF0], [min(dos_raw), max(dos_raw)], color = 'red', linestyle = 'dashed')
ax2r.plot(E_raw, yfe, label = 'fe', color = 'blue', linewidth = 0.5, linestyle = 'dashed')
ax2r.plot(E_raw, yfh, label = 'fh', color = 'red', linewidth = 0.5, linestyle = 'dashed')
ax2.set_ylim([0.0, max(dos_raw) * 1.1])
ax2r.set_ylim([0.0, 2.0])
ax2.set_xlabel("E (eV)", fontsize = fontsize)
ax2.set_ylabel("DOS (states/cm$^3$)", fontsize = fontsize)
ax2r.set_ylabel("fe, fh", fontsize = fontsize)
# ax2.legend()
h2l, l2l = ax2.get_legend_handles_labels()
h2r, l2r = ax2r.get_legend_handles_labels()
ax2.legend(h2l + h2r, l2l + l2r, fontsize = legend_fontsize)
ax3.plot(xEF, yNe, label = 'Ne', color = 'black') #, marker = 'o')
ax3.plot(xEF, yNh, label = 'Nh', color = 'red') #, marker = 'x')
ax3.plot(xEF, yNDp, label = 'ND$^+$', color = 'blue') #, marker = 'x')
ax3.plot(xEF, yNAm, label = 'NA$^-$', color = 'green') #, marker = 'x')
ax3.plot(xEFapprox, yNeapprox, label = 'Ne(Boltz)', color = 'gray', linewidth = 0.5) #, marker = 'x')
ax3.plot(xEFapprox, yNhapprox, label = 'Nh(Boltz)', color = 'purple', linewidth = 0.5) #, marker = 'x')
ax3.plot([EV, EV], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax3.plot([EC, EC], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax3.set_yscale('log')
ax3.set_xlabel("EF (eV)", fontsize = fontsize)
ax3.set_ylabel("N (cm$^{-3}$)", fontsize = fontsize)
ax3.legend(fontsize = legend_fontsize)
ax4.plot(xEF, yRH, label = 'R$_H$', color = 'black', marker = 'o')
ax4.set_xlabel("EF (eV)", fontsize = fontsize)
ax4.set_ylabel("R$_H$ (C$^{-1}$cm$^{-3}$)", fontsize = fontsize)
ax4.legend(fontsize = legend_fontsize)
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([EV, EV], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax1.plot([EC, EC], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax1.set_ylim([0.0, 2.0*T0])
ax1.set_xlabel("EF (eV)", fontsize = fontsize)
ax1.set_ylabel("T$_0$ (K)", fontsize = fontsize)
ax1.legend(fontsize = legend_fontsize)
ax1.tick_params(labelsize = fontsize)
ax2.tick_params(labelsize = fontsize)
ax2r.tick_params(labelsize = fontsize)
ax3.tick_params(labelsize = fontsize)
ax4.tick_params(labelsize = fontsize)
# Rearange the graph axes so that they are not overlapped
plt.tight_layout()
# plt.pause(0.1)
print("Close the graph window to terminate this program")
plt.show()
print("")
usage()
print("")
# print("Press ENTER to exit>>", end = '')
# input()
[ドキュメント]
def main():
"""
スクリプトのメイン実行関数。
コマンドライン引数を処理し、DOSファイルを読み込み、バンド端やドーピング準位を設定する。
その後、`mode`変数に応じて`exec_T`または`exec_EF`を呼び出し、
キャリア計算と結果のプロットを実行する。
:returns: None
"""
global mode, T0, EF0
global Vcell, file, E_raw, dos_raw
global nDOS, Emin, Emax, Estep
global EV, EC
global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
global dEA, NA, dED, ND, EA, ED
updatevars()
print("mode: ", mode)
print("DOS configuration")
print(" Read from [{}]".format(file))
E_raw, dos_raw = np.genfromtxt(file, skip_header=1, usecols=(0,1), unpack=True)
nDOS = len(E_raw)
for i in range(len(E_raw)):
dos_raw[i] *= 1.0 / (Vcell * 1.0e-24)
Emin = E_raw[0]
Emax = E_raw[nDOS-1]
Estep = E_raw[1] - E_raw[0]
print(" E range: {} - {}, {} eV step".format(Emin, Emax, Estep))
print("")
EV, EC = FindBandEdges(E_raw, dos_raw, EF0, Egth)
Eg = EC - EV
print(" Band edge: EV={} EC={} Eg={} eV".format(EV, EC, Eg))
EFmin = EV + dEFmin
EFmax = EC + dEFmax
EFstep = (EFmax - EFmin) / (nEF - 1)
ED = EC + dED
EA = EV + dEA
print(" Donor level : {} eV, {} cm^-3".format(ED, ND))
print(" Acceptor level: {} eV, {} cm^-3".format(EA, NA))
print("")
print("EF dependence")
print(" EF range: {} - {}, {} eV step".format(EFmin, EFmax, EFstep))
print("")
print("T dependence")
print(" Temperature range: {} - {}, {} K step".format(Tmin, Tmax, Tstep))
print("")
print("Integration configuration")
print(" Integration E range: {} - {} eV, or EF+-{}*kBT eV".format(Einteg0, Einteg1, nrange))
print(" quad() parameters")
print(" releps=", eps)
print(" nmaxdiv=", nmaxdiv)
print("")
print("EF at 0 K = ", EF0, " eV")
Ne0K = Ne(0.0, EF0, -5.0, EF0)
# Ne0K = Ne(0.0, EF0, Einteg0, EF0)
print(" Ne at 0 K = ", Ne0K, " /u.c.")
print("")
if mode == 'T':
exec_T()
elif mode == 'EF':
exec_EF()
else:
print("")
print("Error: Invalid mode [{}]".format(mode))
if __name__ == '__main__':
main()