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

プログラムの動作

vasp_correction_bandfilling.py は、第一原理計算ソフトウェアVASPで計算された欠陥モデルの電子構造において発生するバンドフィリング補正エネルギーを推定するPythonスクリプトです。

目的: VASPを用いた欠陥計算では、有限サイズの超格子モデルを採用するため、理想的なバンド構造からの逸脱（バンドフィリング）が生じ、これが欠陥形成エネルギーの誤差の原因となることがあります。本プログラムは、このバンドフィリングに起因するエネルギー補正項を定量的に評価し、欠陥形成エネルギーの正確な決定を支援します。

主な機能:

VASP計算結果の解析: 理想結晶モデルと欠陥結晶モデルそれぞれのVASP計算ディレクトリから、POSCAR (結晶構造)、OUTCAR (総エネルギー、フェルミ準位、スピン状態)、DOSCAR (状態密度)、EIGENVAL (バンドエネルギー、占有数) を読み込みます。
VBM補正の適用: 欠陥結晶モデルのバンドエネルギーおよびフェルミ準位に対し、外部から与えられたVBM補正値 (dEVBM) を適用し、エネルギー基準を調整します。
バンドエッジの特定: 理想結晶モデルの価電子帯上端 (VBM) と伝導帯下端 (CBM) を特定します。
バンドフィリング補正エネルギーの計算: 欠陥結晶の各k点におけるバンドエネルギーと占有数を理想結晶のVBMおよびCBMと比較し、価電子帯の不完全な占有（正孔）および伝導帯への過剰な占有（電子）によるエネルギー補正項 ($\Delta E_h$ および $\Delta E_e$) を算出します。
結果の可視化: 状態密度 (DOS) およびバンド構造のプロットを生成し、理想結晶と補正された欠陥結晶の電子構造を比較表示します。
サマリー出力: 計算された補正エネルギーなどの主要な情報をINI形式のサマリーファイルに保存します。

解決する課題: 欠陥計算における有限サイズの超格子モデルによって、本来バンドギャップ内に出現するはずの欠陥準位が伝導帯や価電子帯に混じってしまったり、バンドエッジが歪んだりすることがあります。これにより、欠陥形成エネルギーの計算に誤差が生じます。このプログラムは、特にバンドエッジ周辺における電子の過不足が引き起こすエネルギー的な寄与を定量化し、欠陥形成エネルギーの補正項を導出することで、より信頼性の高い欠陥形成エネルギー評価を可能にします。

原理

本プログラムは、VASP計算結果からバンドフィリング補正エネルギーを算出します。この補正は、欠陥モデルのバンド構造が理想結晶のバンド構造から逸脱している場合に生じるエネルギー変化を評価するものです。

バンドフィリング補正の概念: バンドフィリング補正は、主に2つの項に分けられます。一つは価電子帯内の正孔による補正 ($\Delta E_h$)、もう一つは伝導帯内の電子による補正 ($\Delta E_e$) です。これらの補正は、欠陥モデルにおいて理想結晶の価電子帯上端 (VBM) や伝導帯下端 (CBM) から逸脱した準位の占有状態に基づいて計算されます。

VBM補正: まず、欠陥結晶モデルの計算結果 (E, EF, Eups, Edns など) に対して、理想結晶のVBMのシフト量 dEVBM を適用します。これにより、欠陥結晶のエネルギー基準が理想結晶のそれに合わせられます。 $$ E_{defect}' = E_{defect} - \Delta E_{VBM} $$ ここで、$E_{defect}'$ はVBM補正後の欠陥結晶のエネルギーです。
理想結晶のバンドエッジ: 理想結晶のVBM ($E_{V,ideal}$) とCBM ($E_{C,ideal}$) は、EIGENVAL または OUTCAR から読み取られ、バンドフィリング補正の基準点として使用されます。
正孔によるエネルギー補正 ($\Delta E_h$): 欠陥結晶のバンドにおいて、理想結晶のVBM ($E_{V,ideal}$) 以下に存在する準位 $E_i$ が完全に占有されていない場合（すなわち、$N_i < N_{emax}$）、その正孔によるエネルギー的な寄与を計算します。 $$ \Delta E_h = \sum_{k, i \text{ s.t. } E_{k,i} \le E_{V,ideal}} w_k (N_{emax} - N_{k,i}) (E_{k,i} - E_{V,ideal}) $$ ここで、$w_k$ はk点の重み、$N_{k,i}$ はk点 $k$ の準位 $i$ の占有数、$N_{emax}$ は一つの準位が占有できる最大電子数 (通常1.0または2.0 for spin-degenerate) です。この項は、VBMより低いエネルギー準位に正孔が存在することで系が安定化する（エネルギーが低下する）効果を負の値として考慮します。
電子によるエネルギー補正 ($\Delta E_e$): 欠陥結晶のバンドにおいて、理想結晶のCBM ($E_{C,ideal}$) 以上に存在する準位 $E_j$ が占有されている場合（すなわち、$N_j > 0$）、その電子によるエネルギー的な寄与を計算します。 $$ \Delta E_e = \sum_{k, j \text{ s.t. } E_{k,j} \ge E_{C,ideal}} w_k N_{k,j} (E_{k,j} - E_{C,ideal}) $$ この項は、CBMより高いエネルギー準位に電子が存在することで系が不安定化する（エネルギーが上昇する）効果を正の値として考慮します。
合計バンドフィリング補正エネルギー: 最終的なバンドフィリング補正エネルギー $\Delta E_{tot}$ は、これら2つの項の合計として計算されます。 $$ \Delta E_{tot} = \Delta E_h + \Delta E_e $$ ISPIN=1 (非スピン分極計算) の場合、正孔および電子の寄与はスピン縮退を考慮して2倍されます。

アルゴリズムの概要:

VASP出力ファイルから、理想結晶と欠陥結晶それぞれのフェルミ準位 ($E_F$)、総エネルギー ($ETOT$)、状態密度 ($DOS(E)$)、バンドエネルギー ($E_{k,i}$) と占有数 ($N_{k,i}$) を読み込む。
理想結晶のVBM ($E_{V,ideal}$) とCBM ($E_{C,ideal}$) を find_band_edges_from_eigenval または gbandedges メソッドで特定する。
欠陥結晶のバンドエネルギーを、ユーザーが指定したVBM補正値 dEVBM を用いてシフトする。
シフトされた欠陥結晶の各k点、各バンドについて以下の条件判定を行う：
- $E_{k,i} \le E_{V,ideal}$ かつ $N_{k,i} < N_{emax}$ の場合、正孔による寄与として $\Delta E_h$ に $w_k (N_{emax} - N_{k,i}) (E_{k,i} - E_{V,ideal})$ を加算する。
- $E_{k,i} \ge E_{C,ideal}$ かつ $N_{k,i} > 0$ の場合、電子による寄与として $\Delta E_e$ に $w_k N_{k,i} (E_{k,i} - E_{C,ideal})$ を加算する。
スピン分極計算でない場合は $\Delta E_h$, $\Delta E_e$, $\Delta E_{tot}$ を2倍する。
matplotlib を用いて、DOSとバンド構造をプロットし、VBM、CBM、フェルミ準位、欠陥準位の位置を視覚化する。
結果をサマリーファイル (.prm) に出力する。

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

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

numpy: 数値計算全般に利用されます。
scipy: 科学計算ライブラリで、特に畳み込み処理 (scipy.interpolate.interp1d, tkconvolution 内で利用) に使用されます。
matplotlib: グラフのプロットに利用されます。
tklib: 作者が開発したと思われるカスタムライブラリで、VASPファイルの読み書き、結晶構造の処理、ユーティリティ機能など、プログラムの中核をなす機能を提供しています。

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

pip install numpy scipy matplotlib

tklib のインストール方法: tklib は標準的なpipでインストールできる公開ライブラリではない可能性があります。このプログラムを実行するためには、tklibライブラリのソースコードを別途入手し、Pythonのモジュール検索パス（例: スクリプトと同じディレクトリ、またはPYTHONPATHに追加されたディレクトリ）に配置する必要があります。具体的な入手方法は、tklibの提供元にご確認ください。

必要な入力ファイル

プログラムの実行には、以下のVASP計算結果ファイルが、それぞれ指定されたディレクトリ内に存在する必要があります。

理想結晶モデルのVASP計算ディレクトリ (CAR_dir1)
- INCAR: VASPの入力パラメータファイル。
- POSCAR: 理想結晶の構造ファイル。セル情報や原子の種類・数が読み込まれます。
- OUTCAR: VASPの計算結果が出力されるファイル。フェルミ準位 (EF)、全エネルギー (TOTEN または TOTEN_zero_sigma)、スピン分極 (ISPIN) などが読み込まれます。
- DOSCAR: 状態密度 (DOS) のデータファイル。全状態密度が読み込まれます。
- EIGENVAL: バンドエネルギーと占有数のデータファイル。バンドエッジの特定やバンドフィリング補正の計算に利用されます。
欠陥結晶モデルのVASP計算ディレクトリ (CAR_dir2)
- INCAR: VASPの入力パラメータファイル。
- POSCAR: 欠陥結晶の構造ファイル。セル情報や原子の種類・数が読み込まれます。
- OUTCAR: VASPの計算結果が出力されるファイル。フェルミ準位 (EF)、全エネルギー (TOTEN または TOTEN_zero_sigma)、スピン分極 (ISPIN) などが読み込まれます。
- DOSCAR: 状態密度 (DOS) のデータファイル。全状態密度が読み込まれます。
- EIGENVAL: バンドエネルギーと占有数のデータファイル。バンドフィリング補正の計算に利用されます。

期待されるファイル形式とデータ構造: これらのファイルはVASPが標準で出力するテキストファイル形式である必要があります。プログラムは tklib.tkcrystal.tkvasp モジュールを使用してこれらのファイルを解析します。

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

プログラムの実行により、以下のファイルと出力が生成されます。

ログファイル (BF_correction-out.txt): プログラムの実行中に標準出力 (コンソール) に表示されるすべての情報が、欠陥結晶モデルの計算ディレクトリ内に BF_correction-out.txt という名前で保存されます。これにより、プログラムの実行ログや計算の途中経過、最終結果を後から確認できます。
サマリーファイル (BF_correction-summary.prm): 計算された主要な結果が、欠陥結晶モデルの計算ディレクトリ内に BF_correction-summary.prm という名前のINI形式ファイルとして保存されます。内容は以下の通りです。
- Car_dir1: 理想結晶モデルのVASP計算ディレクトリパス。
- Car_dir2: 欠陥結晶モデルのVASP計算ディレクトリパス。
- Etot1: 理想結晶モデルの総エネルギー。欠陥モデルの超格子サイズに合わせてスケーリングされる場合があります。
- Etot2: 欠陥結晶モデルの総エネルギー。
- totalNe: 欠陥結晶モデルにおける総電子数。
- dEh: 正孔によるバンドフィリング補正エネルギー ($\Delta E_h$)。
- dEe: 電子によるバンドフィリング補正エネルギー ($\Delta E_e$)。
- dEtot: 合計バンドフィリング補正エネルギー ($\Delta E_{tot} = \Delta E_h + \Delta E_e$)。
例 (BF_correction-summary.prm の内容):
```
[results]
Car_dir1 = /path/to/ideal_vasp_calc
Car_dir2 = /path/to/defect_vasp_calc
Etot1 = -1234.56789
Etot2 = -1230.12345
totalNe = 128.0
dEh = -0.050000
dEe = 0.100000
dEtot = 0.050000
```
グラフウィンドウ: matplotlib を使用して、以下の情報を含むグラフウィンドウが対話的に表示されます。
- 左上 (ax1): 理想結晶と欠陥結晶 (VBM補正適用後) の総状態密度 (TDOS) のプロット。ガウス関数で平滑化されています。理想結晶のVBM, CBM, EF、欠陥結晶のEFが表示されます。
- 左下 (ax2): 欠陥結晶の電子数累積プロット。理想結晶のVBM, CBM、欠陥結晶のEFが表示されます。
- 右 (axband): 欠陥結晶 (VBM補正適用後) のバンド構造プロット。理想結晶のVBM, CBM、欠陥結晶のEFが表示され、各バンドの占有状態がマーカーで示されます。計算された $\Delta E_h$, $\Delta E_e$, $\Delta E_{tot}$ もタイトルに表示されます。
このグラフウィンドウは、ユーザーが手動で閉じない限り表示され続けます。plt.pause(0.1) が呼び出されているため、ウィンドウは即座に閉じず、ユーザーが内容を確認できます。ユーザーはグラフウィンドウの機能を利用して画像を保存できます。

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

vasp_correction_bandfilling.py をコマンドラインから実行する際の書式は以下の通りです。

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

引数の説明:

mode: 実行モードを指定します。現在のバージョンでは 'BF' (Band Filling) のみが有効です。
CAR_dir(ideal): 理想結晶のVASP計算が行われたディレクトリへのパス。
CAR_dir(defect): 欠陥結晶のVASP計算が行われたディレクトリへのパス。
dEVBM: VBM補正値 (単位: eV)。欠陥結晶のバンドエネルギーをこの値だけシフトします。
WG: DOS (状態密度) のガウス関数による平滑化に使用するFWHM (半値全幅) の値 (単位: eV)。
Emin: DOSおよびバンド構造プロットのエネルギー範囲の下限値 (単位: eV)。
Emax: DOSおよびバンド構造プロットのエネルギー範囲の上限値 (単位: eV)。
EF0: バンドエッジ (VBM/CBM) を探索する際の初期フェルミ準位の目安値 (単位: eV)。

例 (usage() 関数の出力):

Usage:
  (a)  python vasp_correction_bandfilling.py mode CAR_dir(ideal) CAR_dir(defect) dEVBM WG Emin Emax EF0
         mode: ['BF']
     ex: python vasp_correction_bandfilling.py BF . . 0.1 0.1 -10.0 10.0 0.1

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

ここでは、仮想的なVASP計算ディレクトリを用いて、vasp_correction_bandfilling.py の具体的な実行例とその結果を説明します。

シナリオ:

理想結晶のVASP計算結果は ./ideal_vasp_calc ディレクトリにあるとします。
欠陥結晶のVASP計算結果は ./defect_vasp_calc ディレクトリにあるとします。
VBM補正値は 0.2 eV とします。
DOSのガウス平滑化FWHMは 0.05 eV とします。
プロットのエネルギー範囲は -5.0 eV から 5.0 eV とします。
バンドエッジ探索のための初期EF0は 0.0 eV とします。

実行コマンド:

python vasp_correction_bandfilling.py BF ./ideal_vasp_calc ./defect_vasp_calc 0.2 0.05 -5.0 5.0 0.0

実行結果の説明:

コンソール出力: コマンドを実行すると、ターミナルまたはコマンドプロンプトに計算の進行状況や解析結果の詳細が表示されます。これには、読み込まれたファイルのパス、理想結晶と欠陥結晶のVASP情報 (フェルミ準位、総エネルギー、ISPINなど)、検出されたバンドエッジのエネルギー、および計算されたバンドフィリング補正エネルギー ($\Delta E_h$, $\Delta E_e$, $\Delta E_{tot}$) が含まれます。
ログファイル (BF_correction-out.txt): 上記コンソールに出力された内容と同一のテキストが、./defect_vasp_calc/BF_correction-out.txt というファイル名で保存されます。
サマリーファイル (BF_correction-summary.prm): ./defect_vasp_calc/BF_correction-summary.prm ファイルが生成され、以下のような内容で保存されます (数値は例です)。
```
[results]
Car_dir1 = ./ideal_vasp_calc
Car_dir2 = ./defect_vasp_calc
Etot1 = -2500.123456
Etot2 = -2495.678901
totalNe = 256.0
dEh = -0.025000
dEe = 0.075000
dEtot = 0.050000
```
この例では、正孔による寄与 $\Delta E_h = -0.025$ eV、電子による寄与 $\Delta E_e = 0.075$ eV、合計のバンドフィリング補正エネルギー $\Delta E_{tot} = 0.050$ eV が算出されています。
グラフウィンドウの表示: matplotlib によって、理想結晶と補正された欠陥結晶のDOSおよびバンド構造を示すインタラクティブなグラフウィンドウが表示されます。このウィンドウには、以下の情報が視覚的に表示されます。
- 理想結晶 (黒線) と欠陥結晶 (赤線) の平滑化されたDOS。
- 理想結晶のVBM (赤破線)、CBM (赤破線)。
- 理想結晶のEF (青破線)、補正後の欠陥結晶のEF (太い青破線)。
- 欠陥結晶の各k点におけるバンドエネルギーと占有状態。
- グラフのタイトル部分には、計算された $\Delta E_h$, $\Delta E_e$, $\Delta E_{tot}$ の値が明記されます。
このグラフウィンドウを閉じることでプログラムは終了します。