"""
VASPのDOSCARファイルに基づいてキャリア密度とフェルミ準位を計算するスクリプト。
概要:
VASPが出力するDOSCARファイルから得られた状態密度(DOS)データを用いて、
半導体のフェルミ準位(EF)とキャリア密度(電子、正孔、イオン化ドーパント)の温度依存性またはEF依存性を計算し、
結果をグラフで可視化します。
詳細説明:
このスクリプトは、VASPのDOSCARから読み取った全状態密度(DOS)と、
INCAR/POSCAR/CONTCARから取得したセル情報を用いて、半導体中の電荷キャリア濃度を計算します。
主に以下の二つのモードで動作します。
1. **温度依存性 ('T'モード)**:
指定された温度範囲で、電荷中性条件を満たすフェルミ準位を数値的に探索します。
各温度における電子濃度、正孔濃度、イオン化ドナー・アクセプター濃度、ホール係数、
および見かけの活性化エネルギーを計算し、結果をプロットします。
フェルミ準位の探索には、バイセクション法またはニュートン法が用いられます。
2. **フェルミ準位依存性 ('EF'モード)**:
指定された固定温度において、フェルミ準位を変化させた場合のキャリア濃度、ホール係数、
および熱平衡時の有効温度 (T0) を計算します。
これにより、伝導帯と価電子帯の有効状態密度を推定することも可能です。
結果は、DOSプロット、キャリア濃度対EFプロット、ホール係数対EFプロット、T0対EFプロットとして表示されます。
スクリプトはVASPの出力ファイルを自動的に検索し、DOSデータを体積あたりに変換します。
また、ユーザーはドナーおよびアクセプターの準位と濃度を設定して、計算に含めることができます。
関連リンク:
:doc:`vasp_ef_usage`
"""
import os
import sys
from math import exp, sqrt
import numpy as np
from math import log, exp
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
from tklib.tkfile import tkFile
import tklib.tkre as tkre
from tklib.tkutils import IsDir, IsFile, SplitFilePath
from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
from tklib.tksci.tksci import Reduce01, Round
from tklib.tksci.tkmatrix import make_matrix1, make_matrix2, make_matrix3
from tklib.tkcrystal.tkcif import tkCIF, tkCIFData
from tklib.tkcrystal.tkcrystal import tkCrystal
from tklib.tkcrystal.tkvasp import tkVASP
from tklib.tksci import tkequation
import tklib.tkcsv
# 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";
#==================================
# global variables
#==================================
Debug = 0
#mode: 'EF', 'T'
mode = 'T'
# Sample information (DOS.data is converted to states/eV/Vcell)
# read from DOSCAR
CAR_path = '.'
DOSCAR_path = "TotalDOS-SnSe.dat"
#DOSCAR_path = "DOS.dat"
UseEF0 = 0
T0 = 300 # K
EF0 = 0.2 # EF at charge neutral
Emid = EF0 # E in bandgap, used to search EC and EV
Egth = 0.01 # Critera DOS value to search band edges
# Other levels
dEA = 0.3 # Acceptor level measured from EV, eV
NA = 0.0 # 1.0e14
dED = -0.28 # Donor level measured from EC, eV
ND = 0.0 # 0.5e17
EA = None # Acceptor states are not considered if EA = None
ED = None # Donor states are not considered if ED = None
# EF range for mode = 'EF'
dEFmin = -0.2 # measured from EV, eV
dEFmax = 0.2 # measured from EC, eV
nEF = 50
# Temperature range for mode = 'T'
Tmin = 300 # K
Tmax = 600
nT = 11
Tstep = None
# 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_bisec = 1.0e-6
nmaxiter_bisection = 200
# Newton method parmeters
dump_newton = 0.5
eps_newton = 1.0e-6
h_newton = 1.0e-6
nmaxiter_newton = 30
iprintinterval = 1
# Global variables
Ne0K = None
Vcell = None # A^3
E_raw = None
dos_raw = None
nDOS = None
Emin = None
Emax = None
Estep = None
EV = None
EC = None
EFmin = None
EFmax = None
EFstep = None
#=============================
# Treat argments
#=============================
[ドキュメント]
def pfloat(str):
"""
概要: 文字列を浮動小数点数に変換します。
詳細説明:
与えられた文字列を浮動小数点数に変換します。
変換に失敗した場合はNoneを返します。
Parameters:
:param str: 変換する文字列。
:type str: str
Returns:
:returns: 変換された浮動小数点数、またはNone。
:rtype: float or None
"""
try:
return float(str)
except:
return None
[ドキュメント]
def pint(str):
"""
概要: 文字列を整数に変換します。
詳細説明:
与えられた文字列を整数に変換します。
変換に失敗した場合はNoneを返します。
Parameters:
:param str: 変換する文字列。
:type str: str
Returns:
:returns: 変換された整数、またはNone。
:rtype: int or None
"""
try:
return int(str)
except:
return None
[ドキュメント]
def getarg(position, defval = None):
"""
概要: コマンドライン引数を取得します。
詳細説明:
sys.argvから指定された位置のコマンドライン引数を取得します。
引数が存在しない場合はデフォルト値を返します。
Parameters:
:param position: 取得する引数の位置（0はスクリプト名）。
:type position: int
:param defval: 引数が見つからなかった場合に返すデフォルト値。デフォルトはNone。
:type defval: any, optional
Returns:
:returns: 指定された位置の引数、またはデフォルト値。
:rtype: str or any
"""
try:
return sys.argv[position]
except:
return defval
[ドキュメント]
def getfloatarg(position, defval = None):
"""
概要: コマンドライン引数を浮動小数点数として取得します。
詳細説明:
getargで取得した文字列をpfloat関数で浮動小数点数に変換して返します。
Parameters:
:param position: 取得する引数の位置。
:type position: int
:param defval: 引数が見つからなかった場合に返すデフォルト値。デフォルトはNone。
:type defval: float, optional
Returns:
:returns: 変換された浮動小数点数、またはNone。
:rtype: float or None
"""
return pfloat(getarg(position, defval))
[ドキュメント]
def getintarg(position, defval = None):
"""
概要: コマンドライン引数を整数として取得します。
詳細説明:
getargで取得した文字列をpint関数で整数に変換して返します。
Parameters:
:param position: 取得する引数の位置。
:type position: int
:param defval: 引数が見つからなかった場合に返すデフォルト値。デフォルトはNone。
:type defval: int, optional
Returns:
:returns: 変換された整数、またはNone。
:rtype: int or None
"""
return pint(getarg(position, defval))
[ドキュメント]
def usage():
"""
概要: スクリプトのコマンドライン引数の使い方を表示します。
詳細説明:
この関数は、プログラムの正しい実行方法と、各モード（'T'または'EF'）での
引数の例を標準出力に表示します。
Parameters:
なし
Returns:
なし
"""
global mode
global CAR_path, DOSCAR_path, Vcell
global E_raw, dos_raw, nDOS, Emin, Emax, Estep
global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
global T0, EF0
if CAR_path == '':
CAR_path = '.'
print("")
print("Usage: Variables in () are optional")
print(" python {} mode (file args)".format(sys.argv[0], ))
print(" (i) python {} T CAR_path DOSCAR_path Tmin Tmax nT)".format(sys.argv[0]))
print(" ex: python {} T {} {} {} {} {}".format(sys.argv[0], CAR_path, DOSCAR_path, Tmin, Tmax, nT))
print(" (ii) python {} EF CAR_path DOSCAR_path nEF EF0 T0)".format(sys.argv[0]))
print(" ex: python {} EF {} {} {} {} {}".format(sys.argv[0], CAR_path, DOSCAR_path, nEF, EF0, T0))
[ドキュメント]
def updatevars():
"""
概要: コマンドライン引数に基づいてグローバル変数を更新します。
詳細説明:
sys.argvを解析し、モードに応じたグローバル変数（CAR_path, DOSCAR_path,
Tmin, Tmax, nT, EF0, T0, nEFなど）の値を設定します。
Parameters:
なし
Returns:
なし
"""
global mode
global CAR_path, DOSCAR_path, Vcell
global E_raw, dos_raw, nDOS, Emin, Emax, Estep
global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
global T0, EF0
global Tmin, Tmax, nT, Tstep
argv = sys.argv
if len(argv) == 1:
usage()
exit()
mode = getarg( 1, mode)
CAR_path = getarg( 2, CAR_path)
DOSCAR_path = getarg( 3, DOSCAR_path)
if mode == 'EF':
nEF = getintarg(4, nEF)
if mode == 'T':
Tmin = getintarg (4, Tmin)
Tmax = getintarg (5, Tmax)
nT = getintarg (6, nT)
EF0 = getfloatarg(7, EF0)
T0 = getfloatarg(8, T0)
Tstep = (Tmax - Tmin) / (nT - 1)
#=============================
# other functions
#=============================
[ドキュメント]
def savecsv(outfile, header, datalist):
"""
概要: データをCSVファイルに保存します。
詳細説明:
指定されたファイルパスにヘッダー行とデータリストを書き込みます。
datalistは各カラムのデータを含むリストのリストとして与えられます。
ファイル書き込みに失敗した場合はエラーメッセージを表示します。
Parameters:
:param outfile: 出力ファイル名。
:type outfile: str
:param header: CSVファイルのヘッダー行。
:type header: list of str
:param datalist: 書き込むデータ。各内部リストは1つのカラムに対応。
:type datalist: list of list
Returns:
なし
"""
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ファイルからデータを読み込みます。
詳細説明:
指定されたCSVファイルからヘッダーと2つのカラム（通常x, y）のデータを読み込みます。
xminとxmaxが指定されている場合、x値がその範囲内のデータのみを抽出します。
ファイル読み込みに失敗した場合はエラーメッセージを表示し、プログラムを終了します。
Parameters:
:param infile: 入力ファイル名。
:type infile: str
:param xmin: 読み込むx値の最小値。Noneの場合、制限なし。
:type xmin: float, optional
:param xmax: 読み込むx値の最大値。Noneの場合、制限なし。
:type xmax: float, optional
:param delimiter: CSVの区切り文字。デフォルトは','。
:type delimiter: str, optional
Returns:
:returns: (header, x_data, y_data)のタプル。
header (list of str): CSVファイルのヘッダー。
x_data (list of float): 最初のカラムのデータ。
y_data (list of float): 2番目のカラムのデータ。
:rtype: tuple
"""
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, Emidgap, DOSth):
"""
概要: バンドエッジ (価電子帯上端 EV, 伝導帯下端 EC) を探索します。
詳細説明:
状態密度(DOS)データと指定されたバンドギャップ中間エネルギーの初期推定値、
およびDOSのしきい値を用いて、価電子帯上端 (EV) と伝導帯下端 (EC) を特定します。
EmidgapからDOSがしきい値DOSthを下回る領域を探索し、そこからバンドエッジを決定します。
Parameters:
:param E: エネルギーの配列。
:type E: list of float
:param DOS: DOSの配列。
:type DOS: list of float
:param Emidgap: バンドギャップ中間エネルギーの初期推定値 (eV)。
:type Emidgap: float
:param DOSth: DOSのしきい値。これを超えるDOS値を持つエネルギーをバンドエッジとみなします。
:type DOSth: float
Returns:
:returns: (EV, EC)のタプル。
EV (float): 価電子帯上端のエネルギー (eV)。
EC (float): 伝導帯下端のエネルギー (eV)。
:rtype: tuple
"""
EV = 1.0e30
EC = 1.0e30
i0 = -1
for i in range(len(E)):
if E[i] >= Emidgap:
i0 = i
break
for i in range(i0, 0, -1):
if DOS[i] > DOSth:
EV = E[i+1]
break
for i in range(i0, len(E), +1):
if DOS[i] > DOSth:
EC = E[i-1]
break
return EV, EC
[ドキュメント]
def rieman(x0, dx, y, xmin, xmax):
"""
概要: リーマン和を用いて積分を近似計算します。
詳細説明:
与えられたデータ配列yとステップサイズdx、開始点x0を用いて、
指定された範囲[xmin, xmax]におけるリーマン和を計算します。
この関数は現在、スクリプト内で使用されていません。
Parameters:
:param x0: データの開始x座標。
:type x0: float
:param dx: x方向のステップサイズ。
:type dx: float
:param y: y値の配列。
:type y: list of float
:param xmin: 積分範囲の下限。
:type xmin: float
:param xmax: 積分範囲の上限。
:type xmax: float
Returns:
:returns: 計算されたリーマン和。
:rtype: 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`で
補間されたDOS関数を定義し、与えられたエネルギー`E`でのDOS値を返します。
エネルギーがDOSデータ範囲外の場合、エラーを発生させプログラムを終了します。
Parameters:
:param E: エネルギー (eV)。
:type E: float
Returns:
:returns: 指定されたエネルギーにおけるDOS (states/eV/Vcell)。
:rtype: float
"""
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)
# Fermi-Dirac distribution function
[ドキュメント]
def fe(E, T, EF):
"""
概要: フェルミ・ディラック分布関数を計算します (電子の場合)。
詳細説明:
指定されたエネルギー`E`、温度`T`、フェルミ準位`EF`における電子の占有確率を計算します。
温度`T`が0の場合、EFより低いエネルギーでは1、高いエネルギーでは0を返します。
Parameters:
:param E: エネルギー (eV)。
:type E: float
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位 (eV)。
:type EF: float
Returns:
:returns: フェルミ・ディラック分布関数の値 (0から1)。
:rtype: 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):
"""
概要: フェルミ・ディラック分布関数の穴の占有確率を計算します (正孔の場合)。
詳細説明:
電子の場合のフェルミ・ディラック分布関数`fe`を用いて、正孔の占有確率 (1 - fe) を計算します。
Parameters:
:param E: エネルギー (eV)。
:type E: float
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位 (eV)。
:type EF: float
Returns:
:returns: 正孔の占有確率 (0から1)。
:rtype: float
"""
return 1.0 - fe(E, T, EF)
# Define the function to be integrated
[ドキュメント]
def DOSfe(E, T, EF):
"""
概要: DOSと電子のフェルミ・ディラック分布関数の積を計算します。
詳細説明:
積分計算のために、エネルギー`E`におけるDOSと電子の占有確率`fe`の積を返します。
Parameters:
:param E: エネルギー (eV)。
:type E: float
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位 (eV)。
:type EF: float
Returns:
:returns: DOS(E) * fe(E, T, EF) の値。
:rtype: float
"""
return DOS(E) * fe(E, T, EF)
[ドキュメント]
def DOSfh(E, T, EF):
"""
概要: DOSと正孔のフェルミ・ディラック分布関数の積を計算します。
詳細説明:
積分計算のために、エネルギー`E`におけるDOSと正孔の占有確率`fh`の積を返します。
Parameters:
:param E: エネルギー (eV)。
:type E: float
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位 (eV)。
:type EF: float
Returns:
:returns: DOS(E) * fh(E, T, EF) の値。
:rtype: float
"""
return DOS(E) * fh(E, T, EF)
[ドキュメント]
def integrate(func, E0, E1, h):
"""
概要: 台形公式を用いて関数を数値積分します。
詳細説明:
指定された関数`func`をエネルギー範囲`[E0, E1]`でステップサイズ`h`を用いて台形公式で数値積分します。
戻り値の2つ目の要素は常に-1.0で、エラー情報は含まれません。
Parameters:
:param func: 積分する関数。引数としてエネルギー (float) を取ります。
:type func: callable
:param E0: 積分範囲の下限 (eV)。
:type E0: float
:param E1: 積分範囲の上限 (eV)。
:type E1: float
:param h: 積分ステップサイズ (eV)。
:type h: float
Returns:
:returns: [積分値, エラーコード]のリスト。エラーコードは常に-1.0。
:rtype: list
"""
n = int((E1 - E0) / h + 1.000001)
h = (E1 - E0) / n
y = [func(E0 + i * h) for i in range(n+1)]
S = 0.5 * (y[0] + y[n]) + sum(y[1:n])
return [h * S, -1.0]
# Calculate the number of electrons from E0 to E1 with the Fermi level EF at T
[ドキュメント]
def Ne(T, EF, E0, E1):
"""
概要: 指定されたエネルギー範囲における電子密度を計算します。
詳細説明:
`DOSfe`関数を`integrate`関数で数値積分し、温度`T`、フェルミ準位`EF`、
エネルギー範囲`[E0, E1]`における電子の総数を計算します。
Parameters:
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位 (eV)。
:type EF: float
:param E0: 積分範囲の下限 (eV)。
:type E0: float
:param E1: 積分範囲の上限 (eV)。
:type E1: float
Returns:
:returns: 計算された電子密度。
:rtype: 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`関数を`integrate`関数で数値積分し、温度`T`、フェルミ準位`EF`、
エネルギー範囲`[E0, E1]`における正孔の総数を計算します。
Parameters:
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位 (eV)。
:type EF: float
:param E0: 積分範囲の下限 (eV)。
:type E0: float
:param E1: 積分範囲の上限 (eV)。
:type E1: float
Returns:
:returns: 計算された正孔密度。
:rtype: 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]
# ionized donor density
[ドキュメント]
def NDp(EF, T):
"""
概要: イオン化ドナー密度を計算します。
詳細説明:
全ドナー濃度`ND`とドナー準位`ED`、温度`T`、フェルミ準位`EF`を用いて、
イオン化されたドナーの密度を計算します。
ドナー準位は電子によって占有されていない場合にイオン化されます (`1 - fe(ED, T, EF)`)。
Parameters:
:param EF: フェルミ準位 (eV)。
:type EF: float
:param T: 温度 (K)。
:type T: float
Returns:
:returns: イオン化ドナー密度。
:rtype: float
"""
global ND, ED, kB
return ND * (1.0 - fe(ED, T, EF))
# ionized acceptor density
[ドキュメント]
def NAm(EF, T):
"""
概要: イオン化アクセプター密度を計算します。
詳細説明:
全アクセプター濃度`NA`とアクセプター準位`EA`、温度`T`、フェルミ準位`EF`を用いて、
イオン化されたアクセプターの密度を計算します。
アクセプター準位は電子によって占有されている場合にイオン化されます (`fe(EA, T, EF)`)。
Parameters:
:param EF: フェルミ準位 (eV)。
:type EF: float
:param T: 温度 (K)。
:type T: float
Returns:
:returns: イオン化アクセプター密度。
:rtype: float
"""
global NA, EA, kB
n = NA * fe(EA, T, EF)
# print("n=", EF, T, NA, EA, n)
return n
# Calculate EF base on
# the charge neutrality (electron number) condition, Ne(T, EF, dE) - N = 0,
# using the Newton method with the initial value EF0
[ドキュメント]
def CalEF(T, EF0, E0, E1, totNe):
"""
概要: 電荷中性条件に基づいてフェルミ準位を計算します。
詳細説明:
ニュートン法(`scipy.optimize.newton`) を使用して、指定された総電子数`totNe`と
電子密度`Ne`が等しくなるようなフェルミ準位`EF`を探索します。
この関数は現在、スクリプト内で使用されていません。
Parameters:
:param T: 温度 (K)。
:type T: float
:param EF0: フェルミ準位の初期推定値 (eV)。
:type EF0: float
:param E0: 電子密度計算の積分範囲下限 (eV)。
:type E0: float
:param E1: 電子密度計算の積分範囲上限 (eV)。
:type E1: float
:param totNe: 総電子数。
:type totNe: float
Returns:
:returns: 計算されたフェルミ準位 (eV)。
:rtype: float
"""
EF = optimize.newton(lambda EF: Ne(T, EF, E0, E1) - totNe, EF0)
return EF
[ドキュメント]
def cal_densities(T, EF, ECmin, ECmax, EVmin, EVmax):
"""
概要: 電子、正孔、イオン化ドナー、イオン化アクセプターの密度を計算します。
詳細説明:
指定された温度`T`とフェルミ準位`EF`、および伝導帯と価電子帯の積分範囲を用いて、
各種キャリア密度と電荷不均衡を計算します。
電荷不均衡は Ne + NAm - Nh - NDp で定義されます。
Parameters:
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位 (eV)。
:type EF: float
:param ECmin: 伝導帯積分範囲の下限 (eV)。
:type ECmin: float
:param ECmax: 伝導帯積分範囲の上限 (eV)。
:type ECmax: float
:param EVmin: 価電子帯積分範囲の下限 (eV)。
:type EVmin: float
:param EVmax: 価電子帯積分範囲の上限 (eV)。
:type EVmax: float
Returns:
:returns: (Ne1, Nh1, NDp1, NAm1, dN) のタプル。
Ne1 (float): 電子密度。
Nh1 (float): 正孔密度。
NDp1 (float): イオン化ドナー密度。
NAm1 (float): イオン化アクセプター密度。
dN (float): 電荷不均衡 (Ne1 + NAm1 - Nh1 - NDp1)。
:rtype: tuple
"""
Ne1 = Ne(T, EF, ECmin, ECmax)
Nh1 = Nh(T, EF, EVmin, EVmax)
NDp1 = NDp(EF, T)
NAm1 = NAm(EF, T)
return Ne1, Nh1, NDp1, NAm1, Ne1 + NAm1 - Nh1 - NDp1
[ドキュメント]
def dQ(T, EF, ECmin, ECmax, EVmin, EVmax, N0):
"""
概要: 電荷中性からの逸脱を計算します。
詳細説明:
`cal_densities`関数を用いて計算された電荷不均衡から、
基準となる過剰電子密度`N0`を差し引いた値を返します。
これはフェルミ準位探索の目的関数として使用されます。
Parameters:
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位 (eV)。
:type EF: float
:param ECmin: 伝導帯積分範囲の下限 (eV)。
:type ECmin: float
:param ECmax: 伝導帯積分範囲の上限 (eV)。
:type ECmax: float
:param EVmin: 価電子帯積分範囲の下限 (eV)。
:type EVmin: float
:param EVmax: 価電子帯積分範囲の上限 (eV)。
:type EVmax: float
:param N0: 基準となる過剰電子密度。
:type N0: float
Returns:
:returns: 電荷不均衡から基準値を差し引いた値。
:rtype: float
"""
Ne1, Nh1, NDp1, NAm1, dN = cal_densities(T, EF, ECmin, ECmax, EVmin, EVmax)
return dN - N0
[ドキュメント]
def diff(h, T, EF, ECmin, ECmax, EVmin, EVmax, N0):
"""
概要: `dQ`関数の数値微分を計算します。
詳細説明:
中心差分法を用いて、`dQ`関数をフェルミ準位`EF`に関して数値微分します。
これはフェルミ準位探索（ニュートン法）のヤコビアンとして使用されます。
Parameters:
:param h: 微分に使用するステップサイズ。
:type h: float
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位 (eV)。
:type EF: float
:param ECmin: 伝導帯積分範囲の下限 (eV)。
:type ECmin: float
:param ECmax: 伝導帯積分範囲の上限 (eV)。
:type ECmax: float
:param EVmin: 価電子帯積分範囲の下限 (eV)。
:type EVmin: float
:param EVmax: 価電子帯積分範囲の上限 (eV)。
:type EVmax: float
:param N0: 基準となる過剰電子密度。
:type N0: float
Returns:
:returns: dQ関数のEFに関する数値微分値。
:rtype: float
"""
return (dQ(T, EF + h, ECmin, ECmax, EVmin, EVmax, N0) - dQ(T, EF - h, ECmin, ECmax, EVmin, EVmax, N0)) / 2.0 / h
[ドキュメント]
def callbackfunc(obj):
"""
概要: 最適化アルゴリズムのイテレーション中に呼び出されるコールバック関数です。
詳細説明:
`tkequation.Equation`クラスのソルバーがイテレーションごとに現在の状態
（イテレーション回数、現在のx値、探索範囲、ステップサイズなど）を
コンソールに出力します。これにより、収束過程を監視できます。
Parameters:
:param obj: 現在の最適化ソルバーオブジェクト。
:type obj: tkequation.Equation
Returns:
:returns: 常に1を返します（ソルバーの続行を指示）。
:rtype: int
"""
if obj.xa is None:
print("Iter {:5d}: x: {:>16.12f}, dx = {:>10.4g}"
.format(obj.iter, obj.x, obj.dx))
else:
print("Iter {:5d}: x: {:>16.12f} in [{:>16.12f}, {:>16.12f}], dx = {:>10.4g}"
.format(obj.iter, obj.x, obj.xa, obj.xb, obj.dx))
return 1
[ドキュメント]
def exec_T():
"""
概要: 温度変化に伴うフェルミ準位とキャリア密度を計算し、プロットします。
詳細説明:
`Tmin`から`Tmax`までの温度範囲で、電荷中性条件を満たすフェルミ準位を
バイセクション法またはニュートン法で探索します。
各温度におけるフェルミ準位、電子濃度、正孔濃度、イオン化ドナー・アクセプター濃度、
ホール係数、および見かけの活性化エネルギーを計算し、コンソールに出力します。
最後に、Matplotlibを用いて結果をグラフで可視化します。
Parameters:
なし
Returns:
なし
"""
global mode, T0, EF0, Ne0K
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, ND
xT = []
xTinv = []
yEF = []
yEint = []
yEentropy = []
yEFa = []
yNe = []
yNh = []
yNDp = []
yNAm = []
yRH = []
EF = EF0
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
if EV - dE < Emin or Emax < EC + dE:
print("*** Integration range error: [{}, {}] exceeds DOS E range [{}, {}]".format(E0, E1, Emin, Emax))
exit()
# Calculate excess electron density N0 if UseEF0 != 0
Ne0, Nh0, NDp0, NAm0, dQ0 = cal_densities(T0, EF0, EC - Estep, EC + dE, EV - dE, EV + Estep)
if UseEF0:
N0 = Ne0 - Nh0
print("Ne0 for EF={:8.4f} eV at {:8.4f} K = {:12.6g} cm^-3 (={:12.4g} - {:12.4g})".format(EF0, T0, N0, Ne0, Nh0))
else:
N0 = 0
if iT == 0:
# if iT == 0 or iT > 0:
EFmin = EF - dEbisec
EFmax = EF + dEbisec
solver = tkequation.Equation(
debug_explicit = 1,
method = 'bisection',
# method = 'newton',
func = lambda EF: dQ(T0, EF, EC - Estep, EC + dE, EV - dE, EV + Estep, N0),
xa = EFmin,
xb = EFmax,
nmaxiter = nmaxiter_bisection,
eps = eps_bisec,
callback = callbackfunc,
isprint = 0
)
else:
solver = tkequation.Equation(
debug_explicit = 1,
method = 'newton',
func = lambda EF: dQ(T0, EF, EC - Estep, EC + dE, EV - dE, EV + Estep, N0),
diff1func = lambda EF: diff(h_newton, T0, EF, EC - Estep, EC + dE, EV - dE, EV + Estep, N0),
x0 = EF,
dump = dump_newton,
nmaxiter = nmaxiter_newton,
eps = eps_newton,
callback = callbackfunc,
isprint = 0
)
x = solver.solve()
if solver.iter >= 0:
print(" Success: Convergence reached at iter = {}, x = {:10.6g}, f = {:8.4g}"
.format(solver.iter, solver.x, solver.f))
else:
print(" Failed: Convergence did not reach")
return 0
EF = solver.x
dQh = solver.f
Neh, Nhh, NDph, NAmh, dQh = cal_densities(T0, EF, EC - Estep, EC + dE, EV - dE, EV + Estep)
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))
print ("%5.0f\t%12.6g" % (T, EF))
print("")
yEaNe = []
yEaNeTh = []
yEaNh = []
yEaNhTh = []
print ("%8s %10s %8s %8s %8s %8s %10s %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]+1.0e-10) - log(yNe[i1]+1.0e-10)) * dT # meV
EaNeTh = (log(yNe[i2]/Th2+1.0e-10) - log(yNe[i1]/Th1+1.0e-10)) * dT # meV
EaNh = (log(yNh[i2]+1.0e-10) - log(yNh[i1]+1.0e-10)) * dT # meV
EaNhTh = (log(yNh[i2]/Th2+1.0e-10) - log(yNh[i1]/Th1+1.0e-10)) * 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.6f %10.6g %10.6g %10.6g %10.6g " % (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 = '$E_F$', color = 'black')
ax1.set_xlabel("T (K)")
ax1.set_ylabel("E$_F$ (eV)")
ax1.legend()
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)")
ax2.set_ylabel("DOS (states/cm$^3$)")
ax2.legend()
ax3.plot(xT, yRH, label = '$R_H$', 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')
ax4.plot(xTinv, yNh, label = '$N_h$', color = 'red', marker = 'o')
ax4.plot(xTinv, yNDp, label = '$N_D^+$', color = 'blue', marker = 'x')
ax4.plot(xTinv, yNAm, label = '$N_A^-$', 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}$)")
ax4.set_ylabel("$N$ (cm$^{-3}$)")
ax4.legend()
ax5.plot(xTinv, yEaNe, label = '$E_a$($N_e$)', color = 'black', marker = 'o')
ax5.plot(xTinv, yEaNeTh, label = '$E_a$($N_eT^{-0.5}$)', color = 'black', marker = 'o')
ax5.plot(xTinv, yEaNh, label = '$E_a$($N_h)', color = 'red', marker = 'o')
ax5.plot(xTinv, yEaNhTh, label = '$E_a$($N_hT^{-0.5})', color = 'red', marker = 'o')
ax5.set_xlabel("1000/$T$ (K$^{-1}$)")
ax5.set_ylabel("$E_a$ (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():
"""
概要: フェルミ準位変化に伴うキャリア密度と関連パラメータを計算し、プロットします。
詳細説明:
指定された固定温度`T0`において、`EFmin`から`EFmax`までのフェルミ準位範囲で、
電子密度、正孔密度、イオン化ドナー・アクセプター密度、ホール係数、
および熱平衡時の有効温度 (T0) を計算します。
さらに有効状態密度を推定します。
計算結果をコンソールに出力し、Matplotlibでグラフとして表示します。
Parameters:
なし
Returns:
なし
"""
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, ND
kBTe = kB * T0 / e
dE = nrange * kBTe
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, Nh1, NDp1, NAm1, dQ = cal_densities(T0, EF, EC - Estep, EC + dE, EV - dE, EV + Estep)
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 ("%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%8s\t%8s" % ("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%8.4g\t%8.4g" % (xEF[i], yNe[i], yNh[i], yNDp[i], yNAm[i], yRH[i], 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 = {:12.6g} cm-3 (T0 = {:12.6g} K)".format(NC, T0Nh))
print(" NV = {:12.6g} cm-3 (T0 = {:12.6g} 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)]
#=============================
# Plot graphs
#=============================
fig = plt.figure(figsize = (8, 8))
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(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)")
ax2.set_ylabel("DOS (states/cm$^3$)")
# ax2.legend()
ax3.plot(xEF, yNe, label = '$N_e$', color = 'black') #, marker = 'o')
ax3.plot(xEF, yNh, label = '$N_h$', color = 'red') #, marker = 'x')
ax3.plot(xEF, yNDp, label = '$N_D^+$', color = 'blue') #, marker = 'x')
ax3.plot(xEF, yNAm, label = '$N_A^-$', color = 'green') #, marker = 'x')
ax3.plot(xEFapprox, yNeapprox, label = '$N_e$(Boltz)', color = 'gray', linewidth = 0.5) #, marker = 'x')
ax3.plot(xEFapprox, yNhapprox, label = '$N_h$(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("$E_F$ (eV)")
ax3.set_ylabel("N (cm$^{-3}$)")
ax3.legend()
ax4.plot(xEF, yRH, label = 'R$_H$', color = 'black', marker = 'o', markersize = 1.5)
ax4.set_xlabel("$E_F$ (eV)")
ax4.set_ylabel("$R_H$ (C$^{-1}$cm$^{-3}$)")
# ax4.legend()
ax1.plot(xEF, yT0Ne, label = '$T_0$($N_e$)', color = 'black') #, marker = 'o')
ax1.plot(xEF, yT0Nh, label = '$T_0$($N_h$)', 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("$E_F$ (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 main():
"""
概要: スクリプトのメインエントリポイントです。
詳細説明:
コマンドライン引数を解析し、VASPファイル（POSCAR, DOSCAR）から
結晶構造情報と状態密度 (DOS) データを読み込みます。
読み込んだ情報に基づいて、バンドエッジ (EV, EC)、ドーパント準位、
および数値積分パラメータを設定します。
その後、グローバル変数`mode`の値に応じて`exec_T`または`exec_EF`関数を呼び出し、
キャリア密度とフェルミ準位の計算を実行します。
Parameters:
なし
Returns:
なし
"""
global mode, T0, EF0, Ne0K
global CAR_path, DOSCAR_path, Vcell
global E_raw, dos_raw, nDOS, Emin, Emax, Estep
global EV, EC
global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
global dEA, NA, dED, ND, EA, ED
global Einteg0
updatevars()
print("")
print("=======================================================")
print(" Calculate EF/T dependence of semiconductor properties")
print("=======================================================")
print("mode: ", mode)
print("*** Read VASP files:")
vasp = tkVASP()
CAR_path = vasp.getdir(CAR_path)
INCAR_path = vasp.get_INCAR(CAR_path)
POSCAR_path = vasp.get_POSCAR(CAR_path)
CONTCAR_path = vasp.get_CONTCAR(CAR_path)
OUTCAR_path = vasp.get_OUTCAR(CAR_path)
print(" CAR dir: ", CAR_path)
print(" INCAR : ", INCAR_path)
print(" POSCAR : ", POSCAR_path)
print(" CONTCAR: ", CONTCAR_path)
print(" OUTCAR : ", OUTCAR_path)
print("")
print("*** Read crystal structure from [{}]".format(POSCAR_path))
cry = vasp.read_poscar(POSCAR_path)
# cry.PrintInf("cell")
a, b, c, alpha, beta, gamm = cry.LatticeParameters()
Vcell = cry.Volume()
print(" cell: {:12.8f} {:12.8f} {:12.8f} A {:10.6f} {:10.6f} {:10.6f}".format(a, b, c, alpha, beta, gamm))
print(" volume: {:12.6f} A^-3".format(Vcell))
print("")
print("*** Read Total DOS from [{}]".format(DOSCAR_path))
E_raw, dos_raw = np.genfromtxt(DOSCAR_path, 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("*** Search EV and EC from the midgap energy {} eV with threshold DOS value of {}".format(Emid, Egth))
EV, EC = FindBandEdges(E_raw, dos_raw, Emid, Egth)
Eg = EC - EV
print(" Band edges: EV={:10.6g} EC={:10.6g} Eg={:10.6g} eV".format(EV, EC, Eg))
EFmin = EV + dEFmin
EFmax = EC + dEFmax
EFstep = (EFmax - EFmin) / (nEF - 1)
print("")
print("*** Additional levels")
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 {} K = {} eV".format(T0, EF0))
Einteg0 = -5.0
Ne0K = Ne(0.0, EF0, Einteg0, EF0)
print(" Ne at {} K from {} eV = {:12.6g} cm^-3.".format(T0, Einteg0, Ne0K))
print("")
if mode == 'T':
exec_T()
elif mode == 'EF':
exec_EF()
else:
terminate("Error: Invalid mode [{}]".format(mode), usage = usage)
if __name__ == '__main__':
main()