vasp_correction_bandfilling.py ドキュメント

概要

vasp_correction_bandfilling.py は、VASP計算結果から欠陥モデルのバンドフィリング補正（Band Filling Correction）を推定するスクリプトです。理想結晶と欠陥結晶の計算結果（DOSや固有値エネルギー）を比較し、バンドアラインメントの補正値である dEVBM を適用した上で、占有状態の違いによる電子および正孔のエネルギー補正量を評価します。

依存ライブラリ

本スクリプトは標準ライブラリに加えて、以下の非標準ライブラリを使用します。

  • numpy

  • scipy

  • matplotlib

  • tklib （ユーザー定義のカスタムライブラリ群）

理論背景とアルゴリズム

欠陥を含むスーパーセルの計算において、浅い欠陥や高濃度の欠陥が存在すると、伝導帯や価電子帯にキャリア（電子または正孔）が入り込み、系全体のエネルギーが人為的に上昇する場合があります。これを補正するのがバンドフィリング補正です（推測）。

1. エネルギー準位のアラインメント補正

欠陥結晶の固有値エネルギー \(E_{i,k}^{(defect)}\) は、理想結晶とのポテンシャルのズレを補正するために、与えられた VBM（価電子帯上端）のシフト量 \(dE_{VBM}\) を用いて以下の通り調整されます（推測）。

\[E_{i,k}^{(defect)} \leftarrow E_{i,k}^{(defect)} - dE_{VBM}\]

2. 正孔によるエネルギー補正 \(dE_h\)

理想結晶の価電子帯上端 \(E_{V1}\) よりも低いエネルギー準位において、状態が完全に占有されていない（空きがある）場合、そこに正孔が存在するとみなし、そのエネルギー上昇分を計算します（推測）。

\[dE_h = \sum_{k} w_k \sum_{i} (N_{emax} - n_{i,k}) (E_{V1} - E_{i,k}^{(defect)})\]

（※コード内では \(dE_h\) に負の符号を含めるような形で加算されています。具体的には \(E_{i,k}^{(defect)} - E_{V1}\) に \(N_{emax} - n_{i,k}\) を乗じて加算しています）

3. 電子によるエネルギー補正 \(dE_e\)

理想結晶の伝導帯下端 \(E_{C1}\) よりも高いエネルギー準位において、電子が占有している場合、そのエネルギー上昇分を計算します（推測）。

\[dE_e = \sum_{k} w_k \sum_{i} n_{i,k} (E_{i,k}^{(defect)} - E_{C1})\]

ここで、\(w_k\) はk点の重み、\(n_{i,k}\) は状態の占有数、\(N_{emax}\) は状態の最大占有数（通常は 1.0）です。

4. 総補正エネルギー \(dE_{tot}\)

スピン分極の有無（ISPIN）に応じて縮退度を考慮し、トータルの補正エネルギーを算出します（推測）。

\[dE_{tot} = c_{spin} (dE_h + dE_e)\]

ここで、\(c_{spin}\) はスピン非分極（ISPIN=1）の場合は \(2\)、スピン分極（ISPIN=2）の場合は \(1\) となります。

コマンドライン引数

スクリプトの実行には、以下のコマンドライン引数を指定します。

python vasp_correction_bandfilling.py mode CAR_dir(ideal) CAR_dir(defect) dEVBM WG Emin Emax EF0

  • mode: 実行モード。現在は BF のみ有効です。

  • CAR_dir(ideal): 理想結晶のVASP出力ファイル群が格納されたディレクトリパス。

  • CAR_dir(defect): 欠陥結晶のVASP出力ファイル群が格納されたディレクトリパス。

  • dEVBM: VBM補正量 \( [eV] \) 。

  • WG: DOSをプロットする際のガウス関数の半値全幅 \( [eV] \) （WG_DOS に対応）。

  • Emin: プロット時のエネルギー下限値 \( [eV] \) 。

  • Emax: プロット時のエネルギー上限値 \( [eV] \) 。

  • EF0: バンドエッジを探索する際の基準エネルギー値 \( [eV] \) 。

入出力仕様

入力ファイル

指定された各ディレクトリ（CAR_dir1 および CAR_dir2）から、以下のVASP出力ファイルが読み込まれます。

  • INCAR

  • POSCAR

  • OUTCAR

  • DOSCAR

  • EIGENVAL

出力ファイル

スクリプトを実行すると、欠陥結晶のディレクトリ内に以下のファイルおよびグラフが出力されます。

  • BF_correction-out.txt 標準出力の内容がリダイレクトされるログファイルです。

  • BF_correction-summary.prm 計算結果のサマリーが保存されるINI形式のファイルです。以下の情報が記録されます。

    • 参照ディレクトリのパス（Car_dir1, Car_dir2）

    • トータルエネルギー（Etot1, Etot2）

    • 全電子数（totalNe）

    • 正孔および電子のエネルギー補正量（dEh, dEe）

    • 総エネルギー補正量（dEtot）

  • グラフプロット Matplotlibを用いたウィンドウが立ち上がり、DOSの比較、状態の占有数、およびバンド構造（エネルギー固有値）の可視化が行われます。

関数詳細

usage()

スクリプトの正しい使用方法やコマンドライン引数のフォーマットを標準出力に表示します。

updatevars()

コマンドライン引数（sys.argv）を読み込み、以下のグローバル変数に値を設定します。

  • mode

  • CAR_dir1

  • CAR_dir2

  • dEVBM

  • WG_DOS

  • Emin

  • Emax

  • EF0

不正なモードが指定された場合は、エラーメッセージを表示してスクリプトを終了します。

BF_correction(mode, CAR_path1, CAR_path2)

バンドフィリング補正のメイン処理を行う関数です。

  1. 理想結晶と欠陥結晶のそれぞれのディレクトリから、結晶構造（POSCAR）、計算情報（OUTCAR）、状態密度（DOSCAR）、および固有値（EIGENVAL）を読み込みます。

  2. 欠陥セルの大きさが理想セルの何倍（nx, ny, nz）かを推定します。

  3. EF0 などのパラメータを元に、バンドエッジ（\(E_V\), \(E_C\)）やHOMO/LUMOのエネルギーを特定します。

  4. 欠陥結晶の固有値やフェルミ準位に対し、dEVBM を用いたエネルギーシフトを適用します。

  5. 各k点および各エネルギーバンドを走査し、キャリアの占有状態から補正量（dEh, dEe）を計算します。

  6. 結果を BF_correction-summary.prm に保存し、MatplotlibによるDOSとバンド構造のグラフ描画を実行します。

main()

スクリプトのエントリポイントです。 updatevars() を呼び出して引数を処理したのち、標準出力を BF_correction-out.txt にリダイレクトする設定を行います。その後、実行モードに応じて BF_correction() を呼び出します。