vasp_ef.py 技術ドキュメント

プログラムの動作

vasp_ef.py は、VASP (Vienna Ab initio Simulation Package) の計算結果であるDOSCARファイルから、半導体のキャリア密度とフェルミ準位を温度またはフェルミ準位の関数として計算するPythonプログラムです。

主な機能:

VASPのDOSCARファイルからエネルギーと状態密度 (DOS) データを読み込み、単位体積あたりにスケーリングします。
POSCARファイルから結晶のセル体積を読み込みます。
状態密度データから価電子帯の頂上 (EV) と伝導帯の底 (EC) を特定し、バンドギャップ ($E_g$) を計算します。
指定された温度 $T$ とフェルミ準位 $E_F$ における電子、正孔、イオン化ドナー、イオン化アクセプターのキャリア密度をフェルミ・ディラック統計に基づいて計算します。
電荷中性条件 ($N_e + N_A^- - N_h - N_D^+ = 0$) を満たすフェルミ準位を、ニュートン法または二分法を用いて数値的に求めます。
モード 'T': 特定の温度範囲でフェルミ準位の温度依存性、キャリア密度 ($N_e$, $N_h$, $N_D^+$, $N_A^-$) の温度依存性、ホール係数 ($R_H$)、および活性化エネルギー ($E_a$) を計算し、グラフとしてプロットします。
モード 'EF': 特定のフェルミ準位範囲でキャリア密度 ($N_e$, $N_h$, $N_D^+$, $N_A^-$) のフェルミ準位依存性、ホール係数 ($R_H$)、および有効状態密度 ($N_C$, $N_V$) の近似値を計算し、グラフとしてプロットします。
計算結果は標準出力に表示され、matplotlib を使用して視覚的にプロットされます。

解決する課題: VASPから得られる状態密度情報を用いて、半導体の重要な電気的特性（フェルミ準位、キャリア密度、ホール係数など）を理論的に評価し、その温度依存性やフェルミ準位依存性を解析することを目的としています。これにより、材料のドーピング効果や温度変化による伝導特性の変化を予測・理解するのに役立ちます。

原理

このプログラムは、半導体の電子状態とキャリア統計に基づいています。

状態密度 (DOS): VASPのDOSCARファイルから得られた状態密度 $DOS_{raw}(E)$ は、結晶のセル体積 $V_{cell}$ を用いて単位体積あたりの状態密度 $DOS(E)$ に変換されます。 $$DOS(E) = \frac{DOS_{raw}(E)}{V_{cell} \cdot 10^{-24}}$$ ここで $V_{cell}$ はÅ$^3$ 単位であり、$10^{-24}$ はÅ$^3$ から cm$^3$ への変換係数です。エネルギー $E$ に対する $DOS(E)$ は、scipy.interpolate.interp1d を用いて補間されます。
フェルミ・ディラック分布関数: 電子が特定のエネルギー準位 $E$ を占有する確率は、温度 $T$ とフェルミ準位 $E_F$ に依存するフェルミ・ディラック分布関数 $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$ は電気素量、$k_B$ はボルツマン定数です。正孔が占有する確率は $f_h(E, T, E_F) = 1 - f_e(E, T, E_F)$ です。
キャリア密度: 電子濃度 $N_e$ および正孔濃度 $N_h$ は、それぞれ以下のように状態密度とフェルミ・ディラック分布関数をエネルギーについて積分することで求められます。 $$N_e(T, E_F) = \int_{E_{min}}^{E_{max}} DOS(E) f_e(E, T, E_F) dE$$N_h(T, E_F) = \int_{E_{min}}^{E_{max}} DOS(E) f_h(E, T, E_F) dE$$ 積分は、プログラム内で定義された integrate 関数（シンプソン法に類似した数値積分）または scipy.integrate.quad (コメントアウトされているが同等の機能) を用いて行われます。
イオン化不純物濃度: ドナー準位 $E_D$ とアクセプター準位 $E_A$ が存在する場合、それらのイオン化密度 $N_D^+$ と $N_A^-$ は、それぞれ以下のように計算されます。 $$N_D^+(T, E_F) = N_D \cdot (1 - f_e(E_D, T, E_F))$$N_A^-(T, E_F) = N_A \cdot f_e(E_A, T, E_F)$$ ここで $N_D$ と $N_A$ は、それぞれの不純物の総濃度です。
電荷中性条件: 半導体中の電荷の総和は中性であるという条件は、以下の式で表されます。 $$N_e(T, E_F) + N_A^-(T, E_F) - N_h(T, E_F) - N_D^+(T, E_F) - N_0 = 0$$ ここで $N_0$ は、プログラム開始時に `UseEF0` が設定されている場合に計算される、中性状態でのキャリア過剰濃度です。この方程式を解いて、$E_F$ を求めます。このプログラムでは、scipy.optimize.newton 関数または二分法（自作の tkequation.Equation クラスが提供）が使用されます。
ホール係数 ($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$ を算出します。 $$E_a = -\frac{k_B}{e} \frac{d(\ln N)}{d(1/T)}$$

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

本プログラムの実行には、以下の非標準Pythonライブラリが必要です。

numpy: 数値計算全般に利用されます。
scipy: 数値積分 (integrate)、最適化 (optimize.newton)、補間 (interpolate.interp1d) に利用されます。
matplotlib: 計算結果をグラフとしてプロットするために利用されます。
tklib: このプログラム固有のユーティリティ関数やVASP関連の読み込み関数が含まれるカスタムライブラリです。pip install tklib でインストールできない場合、このプログラムと同時に提供されるか、別途入手が必要です。

これらのライブラリは、以下の pip コマンドでインストールできます（tklib を除く）。

pip install numpy scipy matplotlib

tklib については、以下のモジュールがインポートされています。

tklib.tkfile
tklib.tkre
tklib.tkutils
tklib.tksci.tksci
tklib.tksci.tkmatrix
tklib.tkcrystal.tkcif
tklib.tkcrystal.tkcrystal
tklib.tkcrystal.tkvasp
tklib.tksci.tkequation
tklib.tkcsv

これらはユーザーの環境に合わせて適切に設定されている必要があります。

必要な入力ファイル

本プログラムは以下のファイルを読み込みます。

DOSCARファイル:
- ファイル名: コマンドライン引数 DOSCAR_path で指定。デフォルトは "TotalDOS-SnSe.dat"。
- 形式: VASPのDOSCARファイル、またはそれに準ずるテキストファイル。1行目にヘッダがあり、2行目以降がエネルギーと状態密度のデータである必要があります。プログラムは、np.genfromtxt を用いて、1行目をスキップし、1列目をエネルギー、2列目を状態密度として読み込みます。
- 内容: 各エネルギー点における状態密度データ。
POSCARファイル:
- ファイル名: コマンドライン引数 CAR_path で指定されたディレクトリ内にある、VASPのPOSCARファイル（例: POSCAR または CONTCAR）。tkVASP クラスが自動的に探索します。
- 形式: VASPのPOSCAR形式。
- 内容: 結晶構造情報。特に、単位セルの体積を計算するために使用されます。単位セルの体積は、状態密度を単位体積あたりの状態密度にスケーリングするために不可欠です。

補足:

CAR_path はVASP関連ファイルが存在するディレクトリを指定します。デフォルトは . (現在のディレクトリ) です。
DOSCAR_path は CAR_path の相対パスまたは絶対パスで指定可能です。

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

本プログラムは、直接的にはファイルへのデータ書き出しを行いません。 CSV出力用の savecsv 関数が定義されていますが、メインの実行パス (exec_T や exec_EF) からは呼び出されていません。

主な出力は以下の形式で提供されます。

標準出力 (コンソール):
- プログラムの進行状況、読み込んだファイル情報、計算設定。
- バンドエッジ ($E_V$, $E_C$)、バンドギャップ ($E_g$) の値。
- 追加のドナー/アクセプター準位の情報。
- モード 'T' の場合: 各温度におけるフェルミ準位 ($E_F$)、キャリア密度 ($N_e$, $N_h$, $N_D^+$, $N_A^-$)、ホール係数 ($R_H$)、および活性化エネルギー ($E_a$) の計算結果が表形式で出力されます。
- モード 'EF' の場合: 各フェルミ準位におけるキャリア密度 ($N_e$, $N_h$, $N_D^+$, $N_A^-$)、ホール係数 ($R_H$)、および $T_0$ ($N_e$), $T_0$ ($N_h$)、有効状態密度 ($N_C$, $N_V$) の計算結果が表形式で出力されます。
matplotlibによるグラフ表示: プログラムの実行後、計算結果を視覚的に表示するmatplotlibのウィンドウが立ち上がります。
- モード 'T' の場合:
  - $T$ vs $E_F$
  - $E$ vs $DOS$ (EF0位置を破線で表示)
  - $T$ vs $R_H$
  - $1000/T$ vs $\log(N)$ (電子、正孔、イオン化ドナー、イオン化アクセプター密度)
  - $1000/T$ vs $E_a$ (活性化エネルギー)
- モード 'EF' の場合:
  - $E$ vs $DOS$ (EF0位置を破線で表示)
  - $E_F$ vs $\log(N)$ (電子、正孔、イオン化ドナー、イオン化アクセプター密度、ボルツマン近似)
  - $E_F$ vs $R_H$
  - $E_F$ vs $T_0$ ($N_e$), $T_0$ ($N_h$)

グラフウィンドウは、ユーザーがEnterキーを押すまで表示され続けます。

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

本プログラムは、以下の形式で実行されます。

python vasp_ef.py mode (file args)

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

(i) 温度依存性モード (`mode = 'T'`)

指定された温度範囲でフェルミ準位とキャリア密度の温度依存性を計算します。

python vasp_ef.py T CAR_path DOSCAR_path Tmin Tmax nT (EF0 T0)

mode: T (温度依存性モード)
CAR_path: VASP関連ファイルがあるディレクトリ (例: . または /path/to/vasp_calc)
DOSCAR_path: 状態密度ファイルへのパス (例: DOS.dat または TotalDOS-SnSe.dat)
Tmin: 計算開始温度 (K)
Tmax: 計算終了温度 (K)
nT: 温度ステップ数 (TminからTmaxまでnT個の点で計算)
(EF0): (オプション) 電荷中性での初期フェルミ準位 (eV)。UseEF0 が設定されている場合に使用されます。デフォルトは 0.2。
(T0): (オプション) EF0 が測定された温度 (K)。UseEF0 が設定されている場合に使用されます。デフォルトは 300。

(ii) フェルミ準位依存性モード (`mode = 'EF'`)

指定されたフェルミ準位範囲でキャリア密度のフェルミ準位依存性を計算します。

python vasp_ef.py EF CAR_path DOSCAR_path nEF EF0 T0

mode: EF (フェルミ準位依存性モード)
CAR_path: VASP関連ファイルがあるディレクトリ (例: . または /path/to/vasp_calc)
DOSCAR_path: 状態密度ファイルへのパス (例: DOS.dat または TotalDOS-SnSe.dat)
nEF: フェルミ準位のステップ数
EF0: 計算を行う固定温度におけるフェルミ準位の中心値 (eV)。これがバンドギャップ中央に近くない場合、dEFmin と dEFmax (プログラム内のグローバル変数) によって調整されます。
T0: 計算を行う固定温度 (K)

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

以下の例では、現在のディレクトリに TotalDOS-SnSe.dat (DOSデータ) と POSCAR (結晶構造データ) が存在することを想定しています。

例1: 温度依存性モードの実行

300 Kから600 Kまで、11ステップでキャリア密度とフェルミ準位の温度依存性を計算します。

python vasp_ef.py T . TotalDOS-SnSe.dat 300 600 11

実行結果の説明: 標準出力には、以下のような情報が順次表示されます。

VASPファイルの読み込みパスとセル体積。
DOSデータのエネルギー範囲。
検出されたバンドエッジ $E_V$, $E_C$ とバンドギャップ $E_g$。
設定されたドナー/アクセプター準位。
各温度ステップ (iT) におけるフェルミ準位の計算過程 (収束情報)。
各温度におけるフェルミ準位 ($E_F$)、電子濃度 ($N_e$)、正孔濃度 ($N_h$)、イオン化ドナー濃度 ($N_D^+$)、イオン化アクセプター濃度 ($N_A^-$)、ホール係数 ($R_H$) の詳細。
活性化エネルギー ($E_a$) を含む要約表。

さらに、matplotlib によって6つのサブプロットを含むグラフウィンドウが表示されます。このウィンドウは、ユーザーがEnterキーを押すまで表示され続けます。

例2: フェルミ準位依存性モードの実行

300 Kで、フェルミ準位を一定範囲で変化させたときのキャリア密度の依存性を50ステップで計算します。

python vasp_ef.py EF . TotalDOS-SnSe.dat 50 0.2 300