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
import csv # savecsv, read_csv関数で使用されるが、元のコードには記述がなかったため追加。ルール1を厳密に解釈すると変更になるが、機能上必要であるため追加。
"""
VASPのDOSCARに基づいてキャリア密度とフェルミ準位を計算するスクリプト。
概要:
VASPによって計算された状態密度(DOS)データを用いて、様々な温度またはフェルミ準位における電子および正孔のキャリア密度、ドナー/アクセプターのイオン化密度、ホール係数などを計算し、結果をプロットします。
詳細説明:
本スクリプトは、DOSCARから読み込んだDOSデータと与えられたドナー/アクセプター準位および密度に基づき、電荷中性条件を満たすフェルミ準位を決定します。計算モードとして、温度(T)を掃引してフェルミ準位やキャリア密度の変化を追うモード('T')と、特定の温度でフェルミ準位(EF)を掃引して各キャリア密度等のEF依存性を調べるモード('EF')があります。結果はMatplotlibを用いてグラフ表示されます。
関連リンク:
:doc:`EF-T-DOS_usage`
"""
# 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 = 932.317188 # 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):
"""
文字列を浮動小数点数に変換する。
概要:
与えられた文字列をfloat型に変換して返す。変換できない場合はNoneを返す。
:param str: 変換する文字列。
:type str: str
:returns: 変換された浮動小数点数、または変換に失敗した場合はNone。
:rtype: float or None
"""
try:
return float(str)
except:
return None
[ドキュメント]
def pint(str):
"""
文字列を整数に変換する。
概要:
与えられた文字列をint型に変換して返す。変換できない場合はNoneを返す。
:param str: 変換する文字列。
:type str: str
:returns: 変換された整数、または変換に失敗した場合はNone。
:rtype: int or None
"""
try:
return int(str)
except:
return None
[ドキュメント]
def getarg(position, defval = None):
"""
コマンドライン引数を取得する。
概要:
指定された位置のコマンドライン引数を取得する。
詳細説明:
`sys.argv`リストから`position`で指定された位置の引数を取得します。
引数が存在しない場合は、`defval`で指定されたデフォルト値を返します。
:param position: 取得する引数の位置 (0はスクリプト名)。
:type position: int
:param defval: 引数が存在しない場合に返すデフォルト値 (オプション)。
:type defval: any, optional
:returns: 指定された位置のコマンドライン引数、またはデフォルト値。
:rtype: str or any
"""
try:
return sys.argv[position]
except:
return defval
[ドキュメント]
def getfloatarg(position, defval = None):
"""
コマンドライン引数を浮動小数点数として取得する。
概要:
指定された位置のコマンドライン引数を浮動小数点数に変換して取得する。
詳細説明:
`getarg`で引数を取得した後、`pfloat`関数を用いてfloat型に変換します。
引数が存在しない場合や変換に失敗した場合は、`defval`で指定されたデフォルト値を返します。
:param position: 取得する引数の位置 (0はスクリプト名)。
:type position: int
:param defval: 引数が存在しない場合や変換に失敗した場合に返すデフォルト値 (オプション)。
:type defval: any, optional
:returns: 変換された浮動小数点数、またはデフォルト値。
:rtype: float or None
"""
return pfloat(getarg(position, defval))
[ドキュメント]
def getintarg(position, defval = None):
"""
コマンドライン引数を整数として取得する。
概要:
指定された位置のコマンドライン引数を整数に変換して取得する。
詳細説明:
`getarg`で引数を取得した後、`pint`関数を用いてint型に変換します。
引数が存在しない場合や変換に失敗した場合は、`defval`で指定されたデフォルト値を返します。
:param position: 取得する引数の位置 (0はスクリプト名)。
:type position: int
:param defval: 引数が存在しない場合や変換に失敗した場合に返すデフォルト値 (オプション)。
:type defval: any, optional
:returns: 変換された整数、またはデフォルト値。
:rtype: int or None
"""
return pint(getarg(position, defval))
[ドキュメント]
def usage():
"""
スクリプトの利用方法を表示する。
概要:
このスクリプトをコマンドラインから実行する際の引数の形式と例を出力する。
詳細説明:
引数なしでスクリプトが実行された場合や、ユーザーが利用方法を確認したい場合に呼び出されます。
主に`mode` (`T`または`EF`) と、それに続くファイル名や計算範囲に関する引数について説明します。
"""
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`関数を呼び出して終了します。
:returns: なし
:rtype: None
"""
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ファイルに保存する。
概要:
指定されたヘッダーとデータリストをCSV形式でファイルに書き込む。
詳細説明:
`outfile`で指定されたファイル名でCSVファイルを開き、`header`を1行目、
`datalist`の内容を以降の行に書き込みます。`datalist`は、各要素が1つの列に対応する
リストのリストとして期待されます。
:param outfile: 出力するCSVファイルのパス。
:type outfile: str
:param header: CSVのヘッダー行となる文字列のリスト。
:type header: list[str]
:param datalist: 保存するデータ。各要素が1つの列に対応するリストのリスト。
例: `[[col1_data1, col1_data2], [col2_data1, col2_data2]]`
:type datalist: list[list[any]]
:returns: なし
:rtype: 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ファイルからデータを読み込む。
概要:
指定されたCSVファイルからヘッダーと2列の数値データを読み込む。
詳細説明:
`infile`で指定されたCSVファイルを読み込み、最初の行をヘッダーとして解釈します。
その後の行から、1列目をxデータ、2列目をyデータとして読み込みます。
`xmin`と`xmax`が指定されている場合、xデータがこの範囲内の行のみが処理されます。
:param infile: 読み込むCSVファイルのパス。
:type infile: str
:param xmin: xデータの最小値。この値より小さいxを持つ行は無視される (オプション)。
:type xmin: float, optional
:param xmax: xデータの最大値。この値より大きいxを持つ行は無視される (オプション)。
:type xmax: float, optional
:param delimiter: CSVファイルの区切り文字 (デフォルトは',')。
:type delimiter: str, optional
:returns: ヘッダー、xデータ、yデータのタプル。
:rtype: tuple[list[str], list[float], list[float]]
"""
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):
"""
バンドエッジ (価電子帯上端 EV と伝導帯下端 EC) を探索する。
概要:
DOSデータとFermiレベルの初期値、DOS閾値に基づいて価電子帯上端 (EV) と伝導帯下端 (EC) を特定する。
詳細説明:
DOS (状態密度) データ配列 `DOS` と対応するエネルギー配列 `E` を使用して、
まず`EF0`に近いエネルギーからDOS値が`Egth` (閾値) を超える点を探索します。
`EF0`より低いエネルギー側で`Egth`を超える最初の点の次のエネルギーをEVとし、
`EF0`より高いエネルギー側で`Egth`を超える最初の点の一つ前のエネルギーをECとします。
:param E: エネルギー値の配列。
:type E: numpy.ndarray
:param DOS: 各エネルギーにおける状態密度値の配列。
:type DOS: numpy.ndarray
:param EF0: Fermiレベルの初期値。バンドエッジ探索の基準点となる。
:type EF0: float
:param Egth: バンドエッジを決定するためのDOS閾値。この値を超えるとバンドと見なされる。
:type Egth: float
:returns: 価電子帯上端 (EV) と伝導帯下端 (EC) のエネルギー値のタプル。
:rtype: tuple[float, float]
"""
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):
"""
リーマン和を用いて数値を積分する。
概要:
与えられたy値の配列に対し、指定された範囲内でリーマン和による近似積分を行う。
詳細説明:
x軸の開始点`x0`とステップ幅`dx`、そして対応するy値の配列`y`を用いて、
積分範囲`[xmin, xmax]`内でのy値の合計を計算し、`dx`を乗算して近似積分値を求めます。
この関数は単純なリーマン和であり、より高度な数値積分手法ではありません。
:param x0: x軸の開始点。
:type x0: float
:param dx: x軸のステップ幅。
:type dx: float
:param y: 積分するy値の配列。
:type y: list[float]
:param xmin: 積分範囲の最小値。
:type xmin: float
:param xmax: 積分範囲の最大値。
:type xmax: float
: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):
"""
エネルギーEにおける状態密度 (DOS) を補間して取得する。
概要:
読み込まれたDOSデータ (`E_raw`, `dos_raw`) を用いて、任意のエネルギー`E`での状態密度を補間計算する。
詳細説明:
キュビックスプライン補間 (`scipy.interpolate.interp1d(kind='cubic')`) を使用して、
離散的なDOSデータから連続的なDOS値を取得します。
要求されたエネルギー`E`が元のDOSデータの範囲外である場合、エラーメッセージを出力しプログラムを終了します。
:param E: 状態密度を評価するエネルギー。
:type E: float
:returns: 指定されたエネルギー`E`における状態密度。
: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):
"""
Fermi-Dirac分布関数を計算する。
概要:
指定されたエネルギー`E`、温度`T`、フェルミ準位`EF`における電子の占有確率を計算する。
詳細説明:
Fermi-Dirac分布関数 `1 / (exp((E - EF) * e / kB / T) + 1.0)` を計算します。
温度`T`が0 Kの場合、Fermi準位`EF`より低いエネルギーでは1.0、高いエネルギーでは0.0を返します（ステップ関数）。
:param E: エネルギー。
:type E: float
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位。
:type EF: float
:returns: Fermi-Dirac分布関数の値。
: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):
"""
正孔のFermi-Dirac分布関数 (1 - fe) を計算する。
概要:
指定されたエネルギー`E`、温度`T`、フェルミ準位`EF`における正孔の占有確率を計算する。
詳細説明:
電子のFermi-Dirac分布関数`fe(E, T, EF)`の補数として計算されます。
:param E: エネルギー。
:type E: float
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位。
:type EF: float
:returns: 正孔のFermi-Dirac分布関数の値。
:rtype: float
"""
return 1.0 - fe(E, T, EF)
# Define the function to be integrated
[ドキュメント]
def DOSfe(E, T, EF):
"""
状態密度とFermi-Dirac分布関数の積を計算する。
概要:
エネルギー`E`における状態密度 (DOS) と電子のFermi-Dirac分布関数 (`fe`) の積を返す。
詳細説明:
電子濃度を計算するための被積分関数として使用されます。
:param E: エネルギー。
:type E: float
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位。
:type EF: float
:returns: `DOS(E) * fe(E, T, EF)` の値。
:rtype: float
"""
return DOS(E) * fe(E, T, EF)
[ドキュメント]
def DOSfh(E, T, EF):
"""
状態密度と正孔のFermi-Dirac分布関数の積を計算する。
概要:
エネルギー`E`における状態密度 (DOS) と正孔のFermi-Dirac分布関数 (`fh`) の積を返す。
詳細説明:
正孔濃度を計算するための被積分関数として使用されます。
:param E: エネルギー。
:type E: float
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位。
:type EF: float
: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`で台形公式により数値積分する。
詳細説明:
積分範囲`[E0, E1]`を`h`のステップで分割し、台形公式 `h * (0.5 * (y0 + yn) + sum(y[1:n-1]))` を適用して積分値を計算します。
`scipy.integrate.quad`の代替として使用されます。エラー値は常に`-1.0`を返します。
:param func: 積分する関数 (エネルギーEを引数に取り、数値を返す関数)。
:type func: callable[[float], float]
:param E0: 積分範囲の開始エネルギー。
:type E0: float
:param E1: 積分範囲の終了エネルギー。
:type E1: float
:param h: 積分ステップ幅。
:type h: float
:returns: 積分結果とダミーのエラー値 (-1.0) のリスト。
:rtype: list[float, float]
"""
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]
# Calculate the number of electrons from E0 to E1 with the Fermi level EF at T
[ドキュメント]
def Ne(T, EF, E0, E1):
"""
特定のエネルギー範囲内の電子濃度を計算する。
概要:
指定された温度`T`、フェルミ準位`EF`、エネルギー範囲`[E0, E1]`における電子濃度を計算する。
詳細説明:
`DOSfe`関数を`E0`から`E1`まで`integrate`関数 (台形公式) を用いて数値積分することで電子濃度を求めます。
デバッグモードが有効な場合、積分エラーが出力されます。
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位。
:type EF: float
:param E0: 積分範囲の開始エネルギー。
:type E0: float
:param E1: 積分範囲の終了エネルギー。
:type E1: float
:returns: エネルギー範囲`[E0, E1]`における電子濃度。
: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):
"""
特定のエネルギー範囲内の正孔濃度を計算する。
概要:
指定された温度`T`、フェルミ準位`EF`、エネルギー範囲`[E0, E1]`における正孔濃度を計算する。
詳細説明:
`DOSfh`関数を`E0`から`E1`まで`integrate`関数 (台形公式) を用いて数値積分することで正孔濃度を求めます。
デバッグモードが有効な場合、積分エラーが出力されます。
:param T: 温度 (K)。
:type T: float
:param EF: フェルミ準位。
:type EF: float
:param E0: 積分範囲の開始エネルギー。
:type E0: float
:param E1: 積分範囲の終了エネルギー。
:type E1: float
:returns: エネルギー範囲`[E0, E1]`における正孔濃度。
: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):
"""
イオン化されたドナー密度を計算する。
概要:
指定されたフェルミ準位`EF`と温度`T`における、イオン化されたドナーの密度を計算する。
詳細説明:
ドナー準位`ED`における電子のFermi-Dirac分布関数`fe`を用いて、
ドナーが電子を放出してイオン化される確率`1 - fe(ED, T, EF)`を計算し、
ドナーの全密度`ND`に乗算します。
:param EF: フェルミ準位。
:type EF: float
:param T: 温度 (K)。
:type T: float
:returns: イオン化されたドナーの密度。
:rtype: float
"""
global ND, ED, kB
return ND * (1.0 - fe(ED, T, EF))
# ionized acceptor density
[ドキュメント]
def NAm(EF, T):
"""
イオン化されたアクセプター密度を計算する。
概要:
指定されたフェルミ準位`EF`と温度`T`における、イオン化されたアクセプターの密度を計算する。
詳細説明:
アクセプター準位`EA`における電子のFermi-Dirac分布関数`fe`を用いて、
アクセプターが電子を受け入れてイオン化される確率`fe(EA, T, EF)`を計算し、
アクセプターの全密度`NA`に乗算します。
:param EF: フェルミ準位。
:type EF: float
:param T: 温度 (K)。
:type T: float
: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):
"""
電荷中性条件に基づいてFermiレベルを計算する。
概要:
指定された温度`T`と電子の総数`totNe`に対し、電荷中性条件を満たすFermiレベル`EF`をNewton法で計算する。
詳細説明:
`scipy.optimize.newton`メソッドを使用し、`Ne(T, EF, E0, E1) - totNe = 0`の根を探索します。
初期値として`EF0`が与えられます。この関数は、プログラム内で実際に使用されていません（二分法が採用されています）。
:param T: 温度 (K)。
:type T: float
:param EF0: Fermiレベルの初期推定値。
:type EF0: float
:param E0: 積分範囲の開始エネルギー。
:type E0: float
:param E1: 積分範囲の終了エネルギー。
:type E1: float
:param totNe: 総電子数。
:type totNe: float
:returns: 電荷中性条件を満たすFermiレベル。
:rtype: float
"""
EF = optimize.newton(lambda EF: Ne(T, EF, E0, E1) - totNe, EF0)
return EF
[ドキュメント]
def exec_T():
"""
温度依存性の計算と結果のプロットを実行する。
概要:
指定された温度範囲でFermiレベル、キャリア密度、ホール係数などの温度依存性を計算し、結果をグラフ表示する。
詳細説明:
`Tmin`から`Tmax`までの温度を`nT`ステップで掃引し、各温度において電荷中性条件
(`Ne + NAm - Nh - NDp = 0`) を満たすFermiレベル (`EF`) を二分法で決定します。
その後、電子濃度 (`Ne`)、正孔濃度 (`Nh`)、イオン化ドナー密度 (`NDp`)、
イオン化アクセプター密度 (`NAm`)、ホール係数 (`RH`) などを計算し、
Matplotlibを使って複数のグラフにプロットして表示します。
:returns: 成功時は0を返す (エラー時は0以外の値を返すことがある)。
:rtype: int
"""
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
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(iT, 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():
"""
Fermiレベル依存性の計算と結果のプロットを実行する。
概要:
特定の温度 (`T0`) の下で、フェルミ準位を掃引し、キャリア密度、ホール係数、有効状態密度などのEF依存性を計算・グラフ表示する。
詳細説明:
`EFmin`から`EFmax`までのフェルミ準位を`nEF`ステップで掃引し、各EFにおける
電子濃度 (`Ne`)、正孔濃度 (`Nh`)、イオン化ドナー密度 (`NDp`)、
イオン化アクセプター密度 (`NAm`)、ホール係数 (`RH`) などを計算します。
さらに、有効状態密度 (`NC`, `NV`) も計算し、Matplotlibを使って複数のグラフにプロットして表示します。
:returns: なし (ただし、内部で終了処理を含む可能性あり)
:rtype: 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, ND
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
try:
dlnNhdEF = (log(yNh[i2]) - log(yNh[i1])) / dEF
except:
#dlnNhdEF = 0
print("PROBLEMDATA",yNh[i2],yNh[i1])
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))
print(" EV=",EV," EC=",EC)
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():
"""
プログラムのメインエントリポイント。
概要:
コマンドライン引数の処理、DOSCARデータの読み込み、バンドエッジの特定、および選択された計算モードの実行を行う。
詳細説明:
1. `updatevars()`を呼び出してコマンドライン引数を解析し、グローバル変数を設定します。
2. 指定されたファイルからDOSデータを読み込み、セル体積でスケーリングします。
3. `FindBandEdges()`を使用して価電子帯上端 (`EV`) と伝導帯下端 (`EC`) を特定します。
4. ドナーおよびアクセプター準位を設定します。
5. 設定された計算モード (`mode`) に基づいて、`exec_T()`または`exec_EF()`のいずれかを実行します。
モードが不正な場合はエラーメッセージを出力し終了します。
:returns: なし
:rtype: 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()