EF-T-DOS.py 技術ドキュメント
プログラムの動作
EF-T-DOS.py は、VASPのDOSCARファイルから変換された状態密度 (DOS) データに基づいて、半導体材料のキャリア密度、フェルミ準位、ホール係数を計算・解析するためのPythonスクリプトです。このプログラムは、以下の2つの動作モードを提供します。
温度スキャンモード (
Tモード)指定された温度範囲で、フェルミ準位、電子濃度 (\(N_e\))、正孔濃度 (\(N_h\))、イオン化ドナー濃度 (\(N_D^+\))、イオン化アクセプター濃度 (\(N_A^-\))、ホール係数 (\(R_H\)) の温度依存性を計算します。
活性化エネルギー (\(E_a\)) も温度の逆数に対して計算され、キャリア輸送メカニズムの分析に役立ちます。
結果はグラフとしてプロットされ、標準出力にも詳細なデータが出力されます。
フェルミ準位スキャンモード (
EFモード)固定された基準温度 (\(T_0\)) で、フェルミ準位を一定範囲で変化させながら、上記のキャリア関連量 (\(N_e, N_h, N_D^+, N_A^-, R_H\)) を計算します。
伝導帯と価電子帯の有効状態密度 (\(N_C, N_V\)) を推定します。
フェルミ分布関数とDOSの積、および各キャリア密度のフェルミ準位依存性などがグラフとしてプロットされ、標準出力にも詳細なデータが出力されます。
このプログラムは、第一原理計算によって得られた状態密度データを用いて、ドーピングされた半導体の電気的特性を理論的に予測し、実験結果との比較や材料設計の指針を得ることを目的としています。
原理
本プログラムは、半導体の電子構造と統計力学に基づいて、キャリア濃度とフェルミ準位を計算します。
物理定数
pi: 円周率h: プランク定数 (\(6.6260755 \times 10^{-34}\) Js)hbar: ディラック定数 (\(1.05459 \times 10^{-34}\) Js)c: 光速 (\(2.99792458 \times 10^8\) m/s)e: 電気素量 (\(1.60218 \times 10^{-19}\) C)kB: ボルツマン定数 (\(1.380658 \times 10^{-23}\) J/K)me: 電子質量 (\(9.1093897 \times 10^{-31}\) kg)
状態密度 (DOS)
入力ファイルから読み込まれる状態密度データ \(DOS(E)\) は、線形補間（scipy.interpolate.interp1d の cubic オプションを使用）によって連続関数として扱われます。このDOSは、単位体積あたりの状態数 (states/eV/cm\(^3\)) に変換されて使用されます。
バンド端の決定
価電子帯上端 \(E_V\) および伝導帯下端 \(E_C\) は、FindBandEdges 関数によって決定されます。この関数は、指定された状態密度閾値 Egth (デフォルト 0.01) を超えるDOSの立ち上がり/立ち下がりのエネルギーをバンド端として識別します。
フェルミ・ディラック分布関数
電子の占有確率を記述するフェルミ・ディラック分布関数 \(f_e(E, T, E_F)\) は以下の式で計算されます。 $\(f_e(E, T, E_F) = \frac{1}{\exp\left(\frac{E - E_F}{k_B T / e_0}\right) + 1}\)\( ここで、\)E\( はエネルギー、\)T\( は絶対温度、\)E_F\( はフェルミ準位です。\)E\(, \)E_F\( は \)eV\( 単位、\)k_B T / e_0\( は \)eV\( 単位の熱エネルギー（ここで \)e_0\( は電気素量の値 \)1.60218 \times 10^{-19}$ をジュールから電子ボルトへの換算係数として使用）。
正孔の占有確率を記述する関数 \(f_h(E, T, E_F)\) は以下の式で計算されます。 $\(f_h(E, T, E_F) = 1 - f_e(E, T, E_F)\)$
キャリア密度
電子濃度 \(N_e\) および正孔濃度 \(N_h\) は、状態密度 \(DOS(E)\) とフェルミ・ディラック分布関数を掛け合わせたものをエネルギーで積分することで求められます。 $\(N_e = \int_{E_{min}}^{E_{max}} DOS(E) \cdot f_e(E, T, E_F) dE\)\( \)\(N_h = \int_{E_{min}}^{E_{max}} DOS(E) \cdot f_h(E, T, E_F) dE\)\( 積分は、プログラム内で実装された台形積分法によって数値的に実行されます。積分範囲は、価電子帯上端 \)E_V\( および伝導帯下端 \)E_C\( を中心に、\)k_B T\( の数倍の範囲 (`nrange` * \)k_B T$) をカバーするように設定されます。
ドーパント濃度
イオン化ドナー密度 (\(N_D^+\)): ドナー準位 \(E_D\) は伝導帯 \(E_C\) から
dED離れた位置にあります。ドナーは電子を放出してイオン化します。 $\(N_D^+ = N_D \left(1 - f_e(E_D, T, E_F)\right)\)$イオン化アクセプター密度 (\(N_A^-\)): アクセプター準位 \(E_A\) は価電子帯 \(E_V\) から
dEA離れた位置にあります。アクセプターは電子を受け入れてイオン化します。 $\(N_A^- = N_A \cdot f_e(E_A, T, E_F)\)$
電荷中性条件とフェルミ準位の決定
フェルミ準位 \(E_F\) は、材料全体の電荷中性条件を満たすように決定されます。 $\(N_e + N_A^- - N_h - N_D^+ = 0\)\( `T` モードでは、この非線形方程式を解くために二分法が用いられます。`EF` モードでは、\)E_F$ を直接スキャンして各量を計算します。
ホール係数
ホール係数 \(R_H\) は、キャリア濃度から以下の式で計算されます。 $\(R_H = \frac{1}{e_0} \frac{N_h - N_e}{(N_h + N_e)^2}\)\( ここで \)e_0$ は電気素量です。
活性化エネルギー (Tモードのみ)
キャリア濃度 \(N\) (電子または正孔) の温度依存性から活性化エネルギー \(E_a\) が計算されます。活性化エネルギーは、\(\ln N\) 対 \(1/T\) のアレニウスプロットの勾配から得られます。 $\(E_a = -k_B \frac{d(\ln N)}{d(1/T)}\)\( また、\)N/\sqrt{T}$ に対する活性化エネルギーも計算されます。
有効状態密度 (EFモードのみ)
ボルツマン近似が適用できるバンド端付近では、キャリア濃度は以下のように近似されます。 $\(N_e \approx N_C \exp\left(-\frac{E_C - E_F}{k_B T / e_0}\right)\)\( \)\(N_h \approx N_V \exp\left(-\frac{E_F - E_V}{k_B T / e_0}\right)\)\( ここで \)N_C\( と \)N_V$ はそれぞれ伝導帯および価電子帯の有効状態密度です。本プログラムでは、フェルミ準位スキャン結果からこれらの有効状態密度を推定します。
必要な非標準ライブラリとインストール方法
このプログラムは以下のPython非標準ライブラリを使用します。
numpyscipymatplotlib
これらのライブラリは、Pythonのパッケージマネージャーである pip を使ってインストールできます。
pip install numpy scipy matplotlib
必要な入力ファイル
プログラムは、VASPのDOSCARから変換された状態密度データファイル DOS.dat を入力として必要とします。
ファイル名:
DOS.dat(デフォルト。コマンドライン引数で変更可能)形式: CSVまたはスペース区切りのテキストファイル。
1行目はヘッダー行としてスキップされます。
1列目はエネルギー (eV)
2列目は状態密度 (states/eV) （単位胞あたりの状態密度を想定）
DOS.dat の例:
# Energy(eV) DOS(states/eV)
-5.000 0.000000e+00
-4.990 1.234567e+18
-4.980 2.345678e+18
...
その他の設定:
Vcell: 単位胞の体積 (A\(^3\))。プログラム内でハードコードされており、DOS.datの状態密度をstates/eV/cm^3に変換するために使用されます。デフォルト値は212.309 A^3です。
生成される出力ファイル
本プログラムは、計算結果をファイルとして直接保存しません。
グラフィック出力:
matplotlibを使用して、計算されたフェルミ準位、キャリア密度、ホール係数、活性化エネルギーなどのグラフが画面に表示されます。ユーザーは表示されたグラフウィンドウから手動で画像を保存できます。標準出力: 計算の進行状況、設定値、および各計算ステップでのフェルミ準位、キャリア密度、ホール係数などの詳細なデータがターミナルに表示されます。
コマンドラインでの使用例 (Usage)
基本的なコマンドラインでの実行形式は以下の通りです。
python EF-T-DOS.py mode (file args)
モードオプション:
温度スキャンモード (
T)python EF-T-DOS.py T file Tmin Tmax nT
T: 温度スキャンモードを指定file: 状態密度データファイル名 (例:DOS.dat)Tmin: 最小温度 (K)Tmax: 最大温度 (K)nT: 温度ステップ数
フェルミ準位スキャンモード (
EF)python EF-T-DOS.py EF file nEF
EF: フェルミ準位スキャンモードを指定file: 状態密度データファイル名 (例:DOS.dat)nEF: フェルミ準位のステップ数
コマンドラインでの具体的な使用例
例1: 温度スキャンモードでキャリア密度とフェルミ準位の温度依存性を計算
この例では、DOS.dat ファイルを使用し、温度300Kから600Kまで11ステップ（50K間隔）でフェルミ準位、キャリア密度、ホール係数、活性化エネルギーを計算し、グラフを表示します。
python EF-T-DOS.py T DOS.dat 300 600 11
実行結果の抜粋:
mode: T
DOS configuration
Read from [DOS.dat]
E range: -5.0 - 5.0, 0.01 eV step
Band edge: EV=0.04 EC=0.55 Eg=0.51 eV
Donor level : 0.27 eV, 5.0e+16 cm^-3
Acceptor level: 0.36 eV, 8.0e+16 cm^-3
EF dependence
EF range: -0.16 - 0.75, 0.0185 eV step
T dependence
Temperature range: 300 - 600, 30 K step
Integration configuration
Integration E range: -3.0 - 3.0 eV, or EF+-6.0*kBT eV
quad() parameters
releps= 1.0e-05
nmaxdiv= 150
EF at 0 K = 0.4 eV
Ne at 0 K = 0.0 /u.c.
T(K) EF(eV) Ne(0 K) (check)
iT 000: 300 K
E range: 0.4 0.55
min (EF,Ne,Nh,NDp,NAm,dQ)=(-0.1600, 1.1077e+00, 1.0664e+19, 4.9962e+16, 8.0000e+16,-1.0506e+19)
max (EF,Ne,Nh,NDp,NAm,dQ)=( 0.7500, 2.0838e+19, 2.5835e+00, 1.2185e+15, 8.0000e+16, 2.0759e+19)
dQ= -1.0506041639353986e+19 2.0759247385929827e+19
dQmin*dQmax= -2.1818296715694726e+38
iter 000: half (EF,Ne,Nh,NDp,NAm,dQ)=(0.2950, 1.4014e+18, 5.7601e+16, 4.0931e+16, 7.9995e+16, 1.3435e+18)
... (中略) ...
iter 025: half (EF,Ne,Nh,NDp,NAm,dQ)=(0.2885, 1.2588e+18, 6.2201e+16, 4.2230e+16, 7.9996e+16, 1.1889e+18)
***(T,EF,Ne,Nh,NDp,NAm,RH)=(300.0000,0.2885,1.2588e+18,6.2201e+16,4.2230e+16,7.9996e+16,5.3444e+02)
300 0.288529
... (各温度での結果が続き、最終的にグラフが表示される) ...
例2: フェルミ準位スキャンモードでキャリア密度や有効状態密度のEF依存性を計算
この例では、DOS.dat ファイルを使用し、デフォルト温度300Kでフェルミ準位を50ステップ変化させながら、キャリア密度、ホール係数、有効状態密度を計算し、グラフを表示します。
python EF-T-DOS.py EF DOS.dat 50
実行結果の抜粋:
mode: EF
DOS configuration
Read from [DOS.dat]
E range: -5.0 - 5.0, 0.01 eV step
Band edge: EV=0.04 EC=0.55 Eg=0.51 eV
Donor level : 0.27 eV, 5.0e+16 cm^-3
Acceptor level: 0.36 eV, 8.0e+16 cm^-3
EF dependence
EF range: -0.16 - 0.75, 0.0185 eV step
T dependence
Temperature range: 300 - 600, 30 K step
Integration configuration
Integration E range: -3.0 - 3.0 eV, or EF+-6.0*kBT eV
quad() parameters
releps= 1.0e-05
nmaxdiv= 150
EF at 0 K = 0.4 eV
Ne at 0 K = 0.0 /u.c.
at 300 K
EF(eV) Ne Nh ND+ NA- RH
-0.1600 1.1077e+00 1.0664e+19 4.9962e+16 8.0000e+16 -1.0506e+01
-0.1414 2.3551e+00 1.0664e+19 4.9962e+16 8.0000e+16 -1.0506e+01
-0.1229 4.9818e+00 1.0664e+19 4.9962e+16 8.0000e+16 -1.0506e+01
... (中略) ...
0.7315 1.9069e+19 3.1508e+00 1.2985e+15 8.0000e+16 -2.6074e-09
0.7500 2.0838e+19 2.5835e+00 1.2185e+15 8.0000e+16 -2.3175e-09
EF(eV) Ne Nh ND+ NA- RH T0(Ne) T0(Nh)
-0.1600 1.1077e+00 1.0664e+19 4.9962e+16 8.0000e+16 -1.0506e+01 3.4862e+02 1.9960e+02
-0.1414 2.3551e+00 1.0664e+19 4.9962e+16 8.0000e+16 -1.0506e+01 3.4862e+02 1.9960e+02
... (中略) ...
Effective density of states at the mid gap 0.295 eV:
NC = 3.398031313364239e+18 cm-3 (T0 = 199.60012291583562 K)
NV = 3.398031313364239e+18 cm-3 (T0 = 199.60012291583562 K)
Calculate fe(E), fh(E)=1-fe(E), DC(E)*fe(E), DV(E)*fh(E) at EF=0.4 eV, T=300 K
... (最終的にグラフが表示される) ...