`vasp_defect.py` プログラム技術ドキュメント

プログラムの動作

vasp_defect.py は、VASP (Vienna Ab initio Simulation Package) の計算結果とユーザーが用意した欠陥形成エンタルピーのExcelファイルに基づいて、半導体のキャリア密度とフェルミ準位を解析するためのPythonスクリプトです。

主な機能:

VASP計算結果の読み込み: 完璧な結晶および欠陥結晶のPOSCAR、OUTCAR、EIGENVAL、DOSCARファイルから、結晶構造情報（体積）、バンド端（価電子帯上端 $E_V$、伝導帯下端 $E_C$）、VASPで計算されたフェルミ準位 $E_{F0}$、全状態密度 (DOS) を読み込みます。
欠陥形成エンタルピーの読み込み: Excel形式の入力ファイルから、様々な電荷状態を持つ欠陥の形成エンタルピー、サイト密度、エントロピー情報を読み込みます。
キャリア密度と欠陥密度の計算: フェルミ・ディラック統計を用いて、特定の温度 $T$ およびフェルミ準位 $E_F$ における電子、正孔、および各欠陥のキャリア密度を計算します。欠陥密度は、欠陥が特定の温度で凍結されているか、それとも電子温度に応じて再計算されるかを考慮して計算されます。
電荷中性条件の解決: 全電荷密度 $\Delta Q = N_h - N_e + \sum_i q_i N_{D,i}$ がゼロとなるフェルミ準位 $E_F$ を、二分法（bisection method）またはニュートン法（Newton method）を用いて数値的に探索します。
遷移準位の特定: 欠陥の電荷状態が変化するフェルミ準位（遷移準位）を計算し、出力します。
結果の出力とプロット: 計算されたキャリア密度、欠陥密度、フェルミ準位などをExcelファイルに出力し、matplotlib を用いてDOS、$\Delta Q$ の $E_F$ 依存性、欠陥形成エンタルピーの $E_F$ 依存性、キャリア密度やフェルミ準位の $T$ 依存性などのグラフをプロットして表示します。

解決する課題:

VASP計算だけでは得られない、温度やフェルミ準位に応じた半導体のキャリア輸送特性や電気的特性、欠陥の安定な電荷状態の解析を可能にします。特に、欠陥の凍結効果を考慮した解析ができるため、実験データとの比較や、より現実的な材料特性の予測に役立ちます。

原理

このプログラムは、半導体中の電子、正孔、および欠陥の平衡キャリア密度とフェルミ準位を計算するために、以下の物理原理と数値アルゴリズムを使用します。

状態密度 (DOS): VASPのDOSCARファイルから、エネルギー $E$ に対する状態密度 $DOS(E)$ を読み込みます。この $DOS(E)$ は、scipy.interpolate.interp1d を用いて補間関数として扱われます。プログラム内部では、価電子帯上端 $E_V$ を基準 ($E_V = 0$) としてエネルギーが正規化されます。
フェルミ・ディラック分布関数: 電子と正孔のエネルギー準位占有確率は、フェルミ・ディラック分布関数 $f_e(E, T, E_F)$ および $f_h(E, T, E_F)$ によって記述されます。 $$f_e(E, T, E_F) = \frac{1}{\exp\left(\frac{(E - E_F)e}{k_B T}\right) + 1}$$f_h(E, T, E_F) = 1 - f_e(E, T, E_F)$$ ここで、$E$ はエネルギー、$T$ は温度、$E_F$ はフェルミ準位、$e$ は素電荷、$k_B$ はボルツマン定数です。
キャリア密度: 電子密度 $N_e$ および正孔密度 $N_h$ は、DOSとフェルミ・ディラック分布関数を特定のエネルギー範囲で積分することによって計算されます。 $$N_e(T, E_F, E_0, E_1) = \int_{E_0}^{E_1} DOS(E) f_e(E, T, E_F) dE$$N_h(T, E_F, E_0, E_1) = \int_{E_0}^{E_1} DOS(E) f_h(E, T, E_F) dE$$ 積分は、プログラム内で実装された台形公式 (`integrate_trapezoid`) を用いて数値的に実行されます。エネルギー範囲 $E_0, E_1$ は、計算効率のために通常はバンド端の近傍に設定されます。
欠陥密度: 欠陥の平衡密度 $N_{D,i}$ は、その形成エネルギー $\Delta E_i$ と、異なる電荷状態を持つ同じ欠陥サイトの分配関数 $Z_S$ に基づいて計算されます。 $$N_{D,i} = \frac{N_{0,i}}{V_{cell} \cdot 10^{-24}} \frac{\exp\left(-\frac{\Delta E_i}{k_B T_e}\right)}{Z_S}$$ ここで、$N_{0,i}$ は欠陥サイトの数（単位胞あたり）、$V_{cell}$ は単位胞体積 (A$^3$)、$\Delta E_i = \Delta H_{0,i} + q_i E_F$ は欠陥の形成エンタルピー $\Delta H_{0,i}$ と電荷 $q_i$ およびフェルミ準位 $E_F$ を含む形成エネルギー、$T_e$ は電子温度です。 $Z_S = \sum_{j} \exp\left(-\frac{\Delta E_j}{k_B T_e}\right)$ は、同じサイトにある全ての電荷状態 $j$ にわたる分配関数です。

プログラムでは、欠陥が特定の凍結温度 $T_{def}$ で凍結される場合（$T \le T_{def}$）、その時点での総欠陥密度が保持され、それ以上の温度では電荷状態の分布のみが電子温度 $T_e$ に応じて再調整されます。$T > T_{def}$ の場合は、欠陥密度も電子温度 $T_e$ で再計算されます。
電荷中性条件: 半導体中の平衡フェルミ準位 $E_F$ は、常に電荷中性条件を満たします。これは、電子、正孔、および全ての欠陥からの全電荷がゼロであることを意味します。 $$\Delta Q(T, E_F) = N_h(T, E_F) - N_e(T, E_F) + \sum_i q_i N_{D,i}(T, E_F) = 0$$ この非線形方程式を解くために、tkequation モジュールの二分法（bisection method）またはニュートン法（Newton method）が用いられます。
遷移準位: ある欠陥が複数の電荷状態を取りうるとき、フェルミ準位の変化に伴って最も安定な電荷状態が変化します。この電荷状態が変化するフェルミ準位を遷移準位と呼びます。これは、二つの異なる電荷状態 $q_1, q_2$ の欠陥形成エネルギーが等しくなる点 $E_F$ として定義されます。 $$\Delta H_{0,1} + q_1 E_F = \Delta H_{0,2} + q_2 E_F$$E_F = \frac{\Delta H_{0,2} - \Delta H_{0,1}}{q_1 - q_2}$$ プログラムは全ての可能な電荷状態間の遷移準位を計算し、出力します。

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

vasp_defect.py の実行には、以下のPython非標準ライブラリが必要です。

numpy: 数値計算を効率的に行うためのライブラリ。
scipy: 科学技術計算ライブラリ。特に数値積分 (scipy.integrate) や補間 (scipy.interpolate.interp1d)、数値最適化 (scipy.optimize) が使用されます。
matplotlib: データのプロットと可視化に使用されます。
tklib: カスタムライブラリであり、ファイル操作、正規表現、ユーティリティ関数、VASPファイルのパーシング、Excel操作、GUI機能など、このプログラムの多くの基盤機能を提供します。

インストール方法

numpy, scipy, matplotlib は pip を使ってインストールできます。

pip install numpy scipy matplotlib

tklib は標準の pip リポジトリにはありません。プログラム内のエラーメッセージによると、tklib ディレクトリがPythonのモジュールパス (PYTHONPATH) に含まれている必要があります。例えば、tklib が /path/to/tkProg/tklib/python にある場合、PYTHONPATH を設定する必要があります。

# bash または zsh の場合
export PYTHONPATH=$PYTHONPATH:/path/to/tkProg/tklib/python

# Windows のコマンドプロンプトの場合
set PYTHONPATH=%PYTHONPATH%;C:\path\to\tkProg\tklib\python

# Windows のPowerShellの場合
$env:PYTHONPATH = $env:PYTHONPATH + ";C:\path\to\tkProg\tklib\python"

上記のコマンドを実行後、プログラムを実行してください。この設定はシェルセッションが終了するとリセットされるため、永続的に設定したい場合は、環境変数設定ファイル（例: .bashrc, .zshrc, システムの環境変数設定）に追加してください。

必要な入力ファイル

プログラムの実行には、VASPの出力ファイルと欠陥形成エンタルピーを記述したExcelファイルが必要です。

欠陥形成エンタルピーファイル (dH_path)
- ファイル名: デフォルトは input.xlsx。コマンドライン引数で指定可能。
- 形式: Excel (.xlsx) ファイル。
- 内容: 欠陥の種類、サイト、電荷、エントロピー、サイト密度、ドープ密度、および異なる点（例えば、異なるひずみ状態や計算方法）における欠陥形成エンタルピーが含まれます。
  - ヘッダー行:
    - atom: 欠陥を構成する原子のラベル (例: Ga, As)
    - site: 欠陥が位置するサイトのラベル (例: Ga1, As1, Vac)
    - charge: 欠陥の電荷 (例: -1, 0, 1)
    - S/kB: 欠陥のエントロピー項（ボルツマン定数 $k_B$ 単位、オプション）
    - N0(sites): 単位胞あたりの欠陥サイトの数
    - Ndoped(cm^-3): (オプション) 特定の欠陥のドープ密度。これは使用されていないか、現在の実装では別の目的で使用されている可能性があります。
    - dH0@point1, dH0@point2, ...: 各計算点における欠陥形成エンタルピー ($\Delta H_0$、単位: eV)。iPoint 引数でどの点を参照するか指定します。
  - データ行: ヘッダーに対応する値。
VASP計算結果ディレクトリ (CAR_path)
- パス: 完璧な結晶モデルのVASP計算結果ファイルが格納されているディレクトリのパス。デフォルトはカレントディレクトリ (.)。
- 内容:
  - POSCAR: 完璧な結晶の構造ファイル。単位胞体積を取得するために使用されます。
  - OUTCAR: VASPのメイン出力ファイル。ISPIN情報 (スピン分極計算かどうか)、VASP計算におけるフェルミ準位 $E_{F0}$ を取得するために使用されます。
  - EIGENVAL: バンド構造計算の結果。価電子帯上端 $E_V$ と伝導帯下端 $E_C$ を決定するために使用されます。
  - DOSCAR: 全状態密度 (DOS) データ。電子・正孔キャリア密度の計算に使用されます。
欠陥結晶のPOSCARディレクトリ (CAR_path_defect)
- パス: 欠陥のある結晶モデルのVASP計算結果ファイルが格納されているディレクトリのパス。デフォルトはカレントディレクトリ (.)。
- 内容:
  - POSCAR: 欠陥のある結晶の構造ファイル。欠陥密度の計算において、完璧な結晶の体積に対する欠陥結晶の体積比を考慮するために使用されます。これにより、欠陥サイトの数を単位体積あたりにスケーリングします。

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

プログラムは、実行モード（EF または T）と入力ファイルのパスに基づいて、いくつかのログファイル、データファイル、およびグラフを生成します。

共通出力ファイル:

ログファイル:
- ファイル名: [dH_pathのファイル名]-defect-out.txt (例: input-defect-out.txt)
- 内容: プログラムの実行中に標準出力に表示される全ての情報（入力パラメータ、VASPファイルから読み取られた情報、バンド端、欠陥情報、計算の途中経過、最終結果など）が詳細に記録されます。

EF モードで生成されるファイル:

N-EFデータファイル:
- ファイル名: [dH_pathのファイル名]-N-EF-[iPointのラベル].xlsx (例: input-N-EF-Point_0.xlsx)
- 内容: 指定された温度 $T_0$ におけるフェルミ準位 $E_F$ の関数として、電子密度 ($N_e$)、正孔密度 ($N_h$)、各欠陥の密度 ($N_D$)、および全電荷 ($\Delta Q$) の数値データがExcel形式で保存されます。
dH-EF (全電荷状態) データファイル:
- ファイル名: [dH_pathのファイル名]-dH-EF-all-[iPointのラベル].xlsx (例: input-dH-EF-all-Point_0.xlsx)
- 内容: フェルミ準位 $E_F$ の関数として、全ての電荷状態における各欠陥の形成エンタルピー ($\Delta H$) がExcel形式で保存されます。
dH-EF (安定電荷状態) データファイル:
- ファイル名: [dH_pathのファイル名]-dH-EF-min-[iPointのラベル].xlsx (例: input-dH-EF-min-Point_0.xlsx)
- 内容: フェルミ準位 $E_F$ の関数として、各欠陥サイトで最も安定な電荷状態の形成エンタルピー ($\Delta H$) がExcel形式で保存されます。遷移準位の情報も含まれます。

T モードで生成されるファイル:

Tデータファイル:
- ファイル名: [dH_pathのファイル名]-T-[iPointのラベル].xlsx (例: input-T-Point_0.xlsx)
- 内容: 温度 $T$ の関数として、平衡フェルミ準位 ($E_F$)、電子密度 ($N_e$)、正孔密度 ($N_h$)、各欠陥の密度 ($N_D$)、および全電荷 ($\Delta Q$) の数値データがExcel形式で保存されます。

グラフウィンドウ:

matplotlib を用いて、以下のグラフが表示されます。
- DOSプロット: エネルギーに対するDOS曲線と、バンド端 ($E_V, E_C$)、および平衡フェルミ準位 ($E_{F,eq}$) を示す垂直線。
- $\Delta Q$ プロット (EFモード): フェルミ準位 $E_F$ に対する全電荷 $\Delta Q$ の絶対値の対数（$\log_{10}|\Delta Q|$）。
- $E_F$ プロット (Tモード): 温度 $T$ または $1000/T$ に対する平衡フェルミ準位 $E_F$。
- $\Delta H$ プロット (EFモード): フェルミ準位 $E_F$ に対する各欠陥の形成エンタルピー ($\Delta H$)。安定な電荷状態のエンタルピーと、オプションで全ての電荷状態のエンタルピーが表示されます。
- キャリア密度プロット: フェルミ準位 $E_F$ （EFモード）または温度 $T$ ($1000/T$) （Tモード）に対する、電子密度 ($N_e$)、正孔密度 ($N_h$)、および各欠陥の密度 ($N_D$) の対数プロット。

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

vasp_defect.py は、実行モードとパラメータをコマンドライン引数として受け取ります。

基本形式:

python vasp_defect.py mode (other args)

EF モード (フェルミ準位依存性解析):

指定された電子温度 $T$ で、フェルミ準位 $E_F$ をスキャンしてキャリア密度と欠陥密度を計算します。

python vasp_defect.py EF [min|all] dH_path CAR_path CAR_path_defect iPoint T(electron) T(defect) EF(defect) EFmin EFmax nEF (dHmin) (dHmax) (ignore_warning)

mode: EF (必須)
[min|all]: プロットモード。min は安定な電荷状態の欠陥形成エンタルピーのみをプロット、all は全ての電荷状態をプロットします。
dH_path: 欠陥形成エンタルピーExcelファイルへのパス (例: input.xlsx)
CAR_path: 完璧な結晶のVASP計算結果ディレクトリへのパス (例: .)
CAR_path_defect: 欠陥結晶のVASP計算結果ディレクトリへのパス (例: .)
iPoint: 欠陥形成エンタルピーExcelファイル内の、どの計算点（列）を使用するかを示すインデックス (0から始まる整数)。
T(electron): 電子の温度 (K)。
T(defect): 欠陥が凍結される温度 (K)。$T(electron) \le T(defect)$ の場合、欠陥密度はこの温度で凍結されます。
EF(defect): 欠陥凍結温度 $T(defect)$ におけるフェルミ準位 (eV)。eq を指定すると、この温度での平衡フェルミ準位が計算されます。数値 (例: 0.5) を直接指定することも可能です。
EFmin: フェルミ準位のプロット範囲の最小値 (eV)。
EFmax: フェルミ準位のプロット範囲の最大値 (eV)。
nEF: $EFmin$ から $EFmax$ までのプロット点数。
(dHmin): (オプション) 欠陥形成エンタルピープロットのY軸最小値 (eV)。
(dHmax): (オプション) 欠陥形成エンタルピープロットのY軸最大値 (eV)。
(ignore_warning): (オプション) サイト数ミスマッチの警告を無視するかどうか (0: 無視しない, 1: 無視する)。

T モード (温度依存性解析):

指定された欠陥凍結温度 $T(defect)$ で、温度 $T$ をスキャンして平衡フェルミ準位とキャリア密度を計算します。

python vasp_defect.py T dH_path CAR_path CAR_path_defect iPoint T(defect) EF(defect) Tmin Tmax nT (ignore_warning)

mode: T (必須)
dH_path: 欠陥形成エンタルピーExcelファイルへのパス (例: input.xlsx)
CAR_path: 完璧な結晶のVASP計算結果ディレクトリへのパス (例: .)
CAR_path_defect: 欠陥結晶のVASP計算結果ディレクトリへのパス (例: .)
iPoint: 欠陥形成エンタルピーExcelファイル内の、どの計算点（列）を使用するかを示すインデックス (0から始まる整数)。
T(defect): 欠陥が凍結される温度 (K)。$T(electron) \le T(defect)$ の場合、欠陥密度はこの温度で凍結されます。
EF(defect): 欠陥凍結温度 $T(defect)$ におけるフェルミ準位 (eV)。eq を指定すると、この温度での平衡フェルミ準位が計算されます。数値 (例: 0.5) を直接指定することも可能です。
Tmin: 温度のプロット範囲の最小値 (K)。
Tmax: 温度のプロット範囲の最大値 (K)。
nT: $Tmin$ から $Tmax$ までのプロット点数。
(ignore_warning): (オプション) サイト数ミスマッチの警告を無視するかどうか (0: 無視しない, 1: 無視する)。

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

`EF` モードの例

コマンド:

python vasp_defect.py EF min input.xlsx . . 0 300 300 eq -1.0 4.0 101

引数の説明:

EF: フェルミ準位の依存性を解析するモード。
min: プロットモード。安定な電荷状態の欠陥形成エンタルピーのみをプロット。
input.xlsx: 欠陥形成エンタルピーデータが記述されたExcelファイル。
.: 完璧な結晶のVASP計算結果（DOSCAR, OUTCARなど）があるカレントディレクトリ。
.: 欠陥結晶のPOSCARがあるカレントディレクトリ。
0: input.xlsx の dH0@point0 列 (0から始まるインデックス) のデータを使用。
300: 電子の温度を 300 K に設定。
300: 欠陥が凍結される温度を 300 K に設定。この場合、電子温度と欠陥凍結温度が等しいため、欠陥は凍結されず、電子温度に応じて再計算されます。
eq: 欠陥凍結温度における平衡フェルミ準位を計算して使用。
-1.0: フェルミ準位の計算・プロット範囲の最小値 -1.0 eV。
4.0: フェルミ準位の計算・プロット範囲の最大値 4.0 eV。
101: フェルミ準位の計算点数を 101 点に設定。

実行結果の概説:

このコマンドを実行すると、プログラムはまずinput.xlsx、カレントディレクトリ内のVASPファイル（POSCAR, OUTCAR, EIGENVAL, DOSCAR）を読み込みます。次に、電子温度300 Kのもと、フェルミ準位を -1.0 eV から 4.0 eV まで101ステップで変化させながら、各フェルミ準位における電子、正孔、およびinput.xlsxで定義された各欠陥のキャリア密度を計算します。サイト数ミスマッチの警告がある場合、ユーザーは継続または停止を選択できます。欠陥凍結温度も300 Kであるため、欠陥密度は電子温度に応じた平衡状態で計算されます。

以下のファイルがカレントディレクトリに生成されます:

input-defect-out.txt: 実行ログ。読み込んだデータ、バンド端、平衡フェルミ準位、計算過程、最終的なキャリア密度と欠陥密度などが詳細に記録されます。
input-N-EF-Point_0.xlsx: フェルミ準位に対する電子・正孔・欠陥密度、全電荷のデータがExcel形式で出力されます。
input-dH-EF-all-Point_0.xlsx: 全ての電荷状態の欠陥形成エンタルピー対EFのデータ。
input-dH-EF-min-Point_0.xlsx: 安定な電荷状態の欠陥形成エンタルピー対EFのデータ。
matplotlibによるグラフウィンドウが表示され、DOS、$\log_{10}|\Delta Q|$ 対 $E_F$、$\Delta H$ 対 $E_F$、キャリア密度対 $E_F$ のプロットが示されます。

`T` モードの例