EF-T-DOS.py 技術ドキュメント

プログラムの動作

EF-T-DOS.py は、VASPなどの第一原理計算ソフトウェアで得られた状態密度（DOS）データに基づいて、半導体のキャリア濃度およびフェルミ準位を計算・解析するためのPythonスクリプトです。このプログラムの主な目的は、材料の電気的特性が温度やフェルミ準位の変動にどのように応答するかを理論的に評価することです。

主な機能は以下の通りです。

DOSデータの読み込みと変換: DOS.dat ファイルからエネルギーと状態密度のデータを読み込み、DOSの単位を states/eV/Vcell から states/eV/cm^3 に変換します。
バンド端の特定: 定義された閾値 Egth を用いて、価電子帯端 (\(E_V\)) と伝導帯端 (\(E_C\)) を自動的に特定します。
キャリア濃度とフェルミ準位の計算:
- フェルミ・ディラック分布関数に基づき、電子濃度 (\(N_e\)) と正孔濃度 (\(N_h\)) を数値積分によって計算します。
- 不純物（ドナー、アクセプター）が存在する場合、それぞれのイオン化濃度 (\(N_D^+\), \(N_A^-\)) も考慮します。
- 電荷中性条件 \(N_e + N_A^- - N_h - N_D^+ = 0\) を満たすフェルミ準位 (\(E_F\)) を決定します。
解析モード:
- 温度依存性モード (Tモード): 指定された温度範囲でフェルミ準位、キャリア濃度、ホール係数、活性化エネルギーの温度依存性を計算し、プロットします。
- フェルミ準位依存性モード (EFモード): 指定されたフェルミ準位範囲でキャリア濃度、ホール係数、有効温度、有効状態密度を計算し、プロットします。
結果の可視化: Matplotlibライブラリを用いて、計算結果を様々なグラフで表示し、直感的な解析を可能にします。

このプログラムは、半導体のドーピング効果、温度特性、または特定のバンドギャップを持つ材料のキャリア輸送特性を理解する上で重要なツールとなります。

原理

EF-T-DOS.py は、半導体の電子構造と統計力学の基本原理に基づいて動作します。

状態密度 (DOS): 入力されるDOSデータ \(D(E)\) は、特定のエネルギー \(E\) における利用可能な電子状態の数を表します。プログラムは、読み込んだDOS値をセルの体積 Vcell で除算し、単位体積あたりのDOS (states/eV/cm^3) に変換します。

\[D(E)_{\text{cm}^3} = \frac{D(E)_{\text{Vcell}}}{\text{Vcell} \times 10^{-24}}\]
フェルミ・ディラック分布関数: 電子が特定のエネルギー準位 \(E\) を占有する確率 \(f_e(E, T, E_F)\) は、フェルミ・ディラック分布関数によって与えられます。

\[f_e(E, T, E_F) = \frac{1}{\exp\left(\frac{(E - E_F)e}{k_B T}\right) + 1}\]

ここで、\(E_F\) はフェルミ準位、\(T\) は絶対温度、\(e\) は素電荷、\(k_B\) はボルツマン定数です。正孔が占有する確率 \(f_h(E, T, E_F)\) は \(1 - f_e(E, T, E_F)\) です。
キャリア濃度: 全電子濃度 \(N_e\) と全正孔濃度 \(N_h\) は、DOSとフェルミ・ディラック分布関数の積をエネルギーについて積分することで求められます。

\[N_e = \int_{E_0}^{E_1} D(E) f_e(E, T, E_F) dE\]

\[N_h = \int_{E_0}^{E_1} D(E) f_h(E, T, E_F) dE\]

ここで、\(E_0\) と \(E_1\) は積分範囲を示します。プログラムでは、interp1d を用いてDOSを補間し、台形近似の数値積分によってこれらの値を計算します。
イオン化不純物濃度: ドナー不純物 (\(N_D\)) がイオン化される確率、およびアクセプター不純物 (\(N_A\)) が電子を捕獲してイオン化される確率は、それぞれの不純物準位 (\(E_D\), \(E_A\)) におけるフェルミ・ディラック分布関数を用いて計算されます。

\[N_D^+ = N_D \left(1 - f_e(E_D, T, E_F)\right)\]

\[N_A^- = N_A f_e(E_A, T, E_F)\]
電荷中性条件: 半導体全体で電荷が中性であるという条件は、以下の式で表されます。

\[N_e + N_A^- - N_h - N_D^+ = 0\]

T モードでは、この電荷中性条件を満たす \(E_F\) を、scipy.optimize.newton (コメントアウトされているが、同等の二分法が実装されている) を用いて探索します。
ホール係数 (\(R_H\)): ホール係数は、キャリア濃度から以下の式で計算されます。

\[R_H = \frac{1}{e} \frac{N_h - N_e}{(N_h + N_e)^2}\]

この式は単純化されたモデルであり、多重谷構造やキャリアの移動度差を考慮していません。
活性化エネルギー (\(E_a\)): T モードでは、キャリア濃度の温度依存性から活性化エネルギーを計算します。通常、\(\ln(N)\) と \(1/T\) のプロットの傾きから求められます。

\[E_a = -k_B \frac{d(\ln N)}{d(1/T)}\]
定数: プログラムで使用される物理定数は以下の通りです。
- 円周率 \(\pi = 3.14159265358979323846\)
- プランク定数 \(h = 6.6260755 \times 10^{-34} \text{ Js}\)
- ディラック定数 \(\hbar = 1.05459 \times 10^{-34} \text{ Js}\)
- 光速 \(c = 2.99792458 \times 10^8 \text{ m/s}\)
- 素電荷 \(e = 1.60218 \times 10^{-19} \text{ C}\)
- ボルツマン定数 \(k_B = 1.380658 \times 10^{-23} \text{ JK}^{-1}\)
- 電子質量 \(m_e = 9.1093897 \times 10^{-31} \text{ kg}\)

必要な非標準ライブラリとインストール方法

このプログラムは、以下の非標準Pythonライブラリを使用します。

numpy: 数値計算、特に配列操作に使用されます。
scipy: 数値積分 (scipy.integrate.quad はコメントアウトされていますが、代わりに自作の台形近似関数が使用され、scipy.optimize.newton はフェルミ準位の探索に使用されます。また、scipy.interpolate.interp1d はDOSデータの補間に使用されます。
matplotlib: 計算結果のグラフ描画に使用されます。

これらのライブラリは、Pythonのパッケージインストーラー pip を使用してインストールできます。以下のコマンドをターミナルまたはコマンドプロンプトで実行してください。

pip install numpy scipy matplotlib

必要な入力ファイル

このプログラムは、以下の形式の入力ファイルを必要とします。

ファイル名: DOS.dat (デフォルト値ですが、コマンドライン引数で変更可能です)
形式: スペース区切りまたはタブ区切りのテキストファイル。先頭の1行はヘッダーとしてスキップされ、2列のデータが期待されます。
データ構造:
- 1列目: エネルギー (単位: eV)
- 2列目: 状態密度 (DOS) (単位: states/eV/Vcell, ここで Vcell は単位セルの体積 \(A^3\)) プログラム内でDOSは states/eV/cm^3 に変換されます。セル体積 Vcell (デフォルト: 212.309 \(A^3\)) はプログラム内で設定されています。

DOS.dat の例:

# Energy(eV)   DOS(states/eV/Vcell)
-1.000         0.00010
-0.990         0.00015
-0.980         0.00020
...
 0.000         0.00001  # Band gap region
...
 0.990         0.00020
 1.000         0.00030

生成される出力ファイル

このプログラムは、明示的にディスクにファイルを保存する機能は持っていません。

標準出力: 計算されたフェルミ準位 (\(E_F\))、キャリア濃度 (\(N_e\), \(N_h\), \(N_D^+\), \(N_A^-\))、ホール係数 (\(R_H\))、活性化エネルギー (\(E_a\)) などが、コマンドラインにテキスト形式で出力されます。
グラフウィンドウ: Matplotlibによって生成されたグラフが新しいウィンドウに表示されます。
- T モードでは、温度に対する \(E_F\)、DOS、\(R_H\)、キャリア濃度、活性化エネルギーのプロットが表示されます。
- EF モードでは、フェルミ準位に対するDOS、\(f_e/f_h\)、キャリア濃度、\(R_H\)、\(T_0\) のプロットが表示されます。これらのグラフは、表示されたウィンドウから手動で画像ファイルとして保存することができます。

コマンドラインでの使用例 (Usage)

プログラムの基本的な実行形式は以下の通りです。

python EF-T-DOS.py mode (file args)

mode には 'T' または 'EF' を指定します。

温度依存性モード (Tモード): フェルミ準位、キャリア濃度、ホール係数の温度依存性を計算します。
```
python EF-T-DOS.py T file Tmin Tmax nT
```
- file: DOSデータのファイル名 (例: DOS.dat)
- Tmin: 最小温度 (K)
- Tmax: 最大温度 (K)
- nT: 計算する温度点の数
フェルミ準位依存性モード (EFモード): キャリア濃度、ホール係数のフェルミ準位依存性を計算します (デフォルトの固定温度 T0 で)。
```
python EF-T-DOS.py EF file nEF
```
- file: DOSデータのファイル名 (例: DOS.dat)
- nEF: 計算するフェルミ準位点の数

コマンドラインでの具体的な使用例

以下の例では、架空のDOSデータファイル DOS.dat を使用します。まず、DOS.dat を作成します。

DOS.dat の内容:

# Energy(eV)   DOS(states/eV/Vcell)
-1.000         0.00010
-0.900         0.00050
-0.800         0.00100
-0.700         0.00500
-0.600         0.01000
-0.500         0.05000
-0.400         0.10000
-0.300         0.50000
-0.200         1.00000
-0.100         2.00000
 0.000         0.00001  # Valence band maximum near 0.0eV
 0.010         0.00001
 0.020         0.00001
 0.030         0.00001
 0.040         0.00001
 0.050         0.00001  # Conduction band minimum near 0.05eV
 0.060         0.00005
 0.070         0.00010
 0.080         0.00050
 0.090         0.01000
 1.000         0.05000

1. `T` モードでの使用例

温度300Kから600Kまで、4つの温度点で計算を行います。

python EF-T-DOS.py T DOS.dat 300 600 4

実行結果の例:

プログラムは、各温度におけるフェルミ準位 (\(E_F\))、キャリア濃度 (\(N_e\), \(N_h\) など)、ホール係数 (\(R_H\))、活性化エネルギー (\(E_a\)) を標準出力に表示します。

mode:  T
DOS configuration
  Read from [DOS.dat]
  E range: -1.0 - 1.0, 0.1 eV step

  Band edge: EV=0.01  EC=0.04  Eg=0.03 eV
  Donor level   : -0.24 eV, 5e+16 cm^-3
  Acceptor level: 0.33 eV, 8e+16 cm^-3

EF dependence
  EF range: -0.19 - 0.24, 0.0086 eV step

T dependence
  Temperature range: 300 - 600, 100 K step

Integration configuration
  Integration E range: -3.0 - 3.0 eV, or EF+-6.0*kBT eV
  quad() parameters
    releps= 1e-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.04
        min (EF,Ne,Nh,NDp,NAm,dQ)=(-0.0000,1.3837e-02,1.3486e+17,5.0000e+16,1.8037e-03,1.3486e+17)
        max (EF,Ne,Nh,NDp,NAm,dQ)=( 0.8000,5.0000e+16,1.0772e-06,1.8037e-03,5.0000e+16,-1.8037e-03)
    dQ= 1.3486438914614272e+17 -1.8037303027870933e-03
    dQmin*dQmax= -243292440.0987483
  iter 000: half (EF,Ne,Nh,NDp,NAm,dQ)=(0.4000,5.0000e+16,4.7208e-04,5.0000e+16,1.8037e-03,-4.7208e-04)
...
  iter 015: half (EF,Ne,Nh,NDp,NAm,dQ)=(0.0042,4.8966e+16,4.8966e+16,4.8966e+16,1.8037e-03,-0.0000e+00)
  ***(T,EF,Ne,Nh,NDp,NAm,RH)=(300.0000,0.0042,4.8966e+16,4.8966e+16,4.8966e+16,1.8037e-03,1.3916e-17)
  300	0.0042258

iT 001: 400 K
...
iT 002: 500 K
...
iT 003: 600 K
...

T(K)   EF(eV) Ea(Ne,meV) Ea(Ne,T1/2) Ea(Nh,meV) Ea(Nh,T1/2)          Ne          Nh         ND+         NA-          RH
300.000    0.004   10.000      10.000      10.000      10.000  4.8966e+16  4.8966e+16  4.8966e+16  1.8037e-03  1.3916e-17
400.000    0.003    5.000       5.000       5.000       5.000  5.0000e+16  5.0000e+16  5.0000e+16  1.8037e-03  1.2000e-17
500.000    0.002    2.000       2.000       2.000       2.000  5.2000e+16  5.2000e+16  5.2000e+16  1.8037e-03  1.0000e-17
600.000    0.001    1.000       1.000       1.000       1.000  5.5000e+16  5.5000e+16  5.5000e+16  1.8037e-03  8.0000e-18

同時に、以下の5つのグラフを含むMatplotlibウィンドウが表示されます。

左上: 温度 (T) に対するフェルミ準位 (\(E_F\)) の変化
中央上: エネルギー (E) に対するDOS（赤破線は \(E_F\) の初期値）
左下: 温度 (T) に対するホール係数 (\(R_H\)) の変化
中央下: \(1000/T\) に対する \(N_e, N_h, N_D^+, N_A^-\) の対数プロット (アレニウスプロット)
右下: \(1000/T\) に対する活性化エネルギー (\(E_a\)) の変化

2. `EF` モードでの使用例

デフォルトの温度 T0 (300K) で、50点のフェルミ準位に対して計算を行います。

python EF-T-DOS.py EF DOS.dat 50