vasp_correction_VBM.py 技術ドキュメント
プログラムの動作
vasp_correction_VBM.py は、VASP計算を用いて欠陥を含む結晶のVBM (Valence Band Maximum) 補正量を推定するためのPythonスクリプトです。このプログラムは、欠陥スーパーセル計算において生じる有限サイズ効果によるVBMのエネルギーシフトを補正することを目的としています。
主な機能は以下の通りです。
VASP計算結果の読み込み: 理想結晶モデルと欠陥結晶モデルそれぞれのVASP計算出力ファイル (
POSCAR,OUTCAR,DOSCAR,EIGENVAL) を読み込みます。平均コアポテンシャルの計算: 理想結晶と欠陥結晶における各原子サイトの平均コアポテンシャル (\(V_{core}\)) をOUTCARファイルから抽出します。
欠陥サイトの特定と除外: 結晶構造の比較に基づいて空孔、格子間原子、置換欠陥を特定し、指定された距離 (
Rmin_potential) 以内の欠陥近傍サイトをポテンシャル補正計算から除外します。VBM補正量の決定: 欠陥から十分に離れたサイトの\(V_{core}\)の差(\(\Delta V_{core}\))の分布を解析し、その最頻値 (
'mode'アルゴリズム) または平均値 ('average'アルゴリズム) をVBM補正量 (\(dE_{VBM}\)) として決定します。DOSのエネルギーシフト: 欠陥結晶のDOSを\(dE_{VBM}\)でエネルギーシフトし、理想結晶のDOSと比較してプロットします。
サマリーとログの出力: 計算結果と設定をまとめたサマリーファイルと、実行ログファイルを出力します。
グラフ表示: \(V_{core}\)と\(\Delta V_{core}\)の欠陥からの距離依存性、およびDOSを可視化します。
原理
本プログラムは、欠陥計算におけるスーパーセルの有限サイズ効果によって生じるVBMのエネルギーシフトを補正する「静電ポテンシャルアライメント法」に基づいています。
静電ポテンシャルの差の利用: 理想結晶と欠陥結晶の2つのVASP計算結果から、それぞれの系の各原子サイトにおける平均コアポテンシャル \(V_{core}\) をOUTCARファイルから読み出します。これらのポテンシャルの差 \(\Delta V_{core}\) を、欠陥からの距離 \(r\) の関数として調べます。 $\(\Delta V_{core}(i) = V_{core, defect}(i) - V_{core, ideal}(i)\)\( ここで、\)V_{core, defect}(i)\( は欠陥スーパーセルの \)i\( 番目の原子サイトの平均コアポテンシャル、\)V_{core, ideal}(i)\( は理想スーパーセルの対応する \)i$ 番目の原子サイトの平均コアポテンシャルです。
欠陥領域の除外: 欠陥そのものやその周辺の原子サイトは、ポテンシャルが大きく乱されるため、\(\Delta V_{core}\) の計算には不適切です。したがって、以下の基準でサイトを除外します。
原子構成の不一致: 理想モデルと欠陥モデルで原子種が異なるサイト（置換欠陥など）は除外されます。
欠陥からの距離: ユーザーが指定する
Rmin_potential(アンダーシュテーム) 以下の距離にあるサイトは、欠陥の影響が強いとみなし除外されます。また、Rmin_vacancyは欠陥の特定に使用されます。
VBM補正量 \(dE_{VBM}\) の決定: 欠陥の影響が小さい遠方のサイト群について、\(\Delta V_{core}\) の値を集計します。
'mode'アルゴリズム: 遠方サイトの \(\Delta V_{core}\) 値の分布を、ガウス関数で畳み込み (convolution) して平滑化します。この分布の最頻値を \(dE_{VBM}\) とします。畳み込みには、Wgで指定されたFWHM (半値全幅) を持つガウス関数を使用します。 $\(G(x; \mu, \sigma) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)\)\( ここで \)\sigma = W_g / (2\sqrt{2\ln 2})$ です。'average'アルゴリズム: 遠方サイトの \(\Delta V_{core}\) の平均値を繰り返し計算し、平均値からdVthで指定された閾値を超える偏差を持つサイトを順次除外していきます。このプロセスを収束するまで繰り返し、最終的に得られた平均値を \(dE_{VBM}\) とします。
DOSの補正: 決定された \(dE_{VBM}\) を用いて、欠陥モデルのDOSのエネルギー軸をシフトします。 $\(E_{corrected} = E_{original} - dE_{VBM}\)$ これにより、理想結晶のVBMと欠陥結晶のVBMが一致するように調整されます。DOSのプロット時にも、
WG_DOSで指定されたFWHMを持つガウス関数で平滑化が行われます。
必要な非標準ライブラリとインストール方法
このプログラムは以下の非標準ライブラリを使用しています。
numpy: 数値計算のための基本的なライブラリです。
scipy: 科学技術計算のためのライブラリで、主に補間 (
interp1d) や畳み込み (convolution) に使用されます。matplotlib: グラフ描画のためのライブラリです。
tklib: VASPの出力ファイルを解析したり、結晶構造を操作したりするためのカスタムライブラリです。
tkApplication,tkVASP,tkCrystalなどがこれに含まれます。
インストール方法:
numpy, scipy, matplotlib はPythonの標準的なパッケージマネージャ pip を使用してインストールできます。
pip install numpy scipy matplotlib
tklib は一般的なPythonパッケージリポジトリには登録されていないカスタムライブラリです。別途入手し、Pythonのモジュール検索パス上に配置するか、sys.pathに追加する必要があります。通常、配布元から提供されるか、ローカルにクローンされたディレクトリをPython環境に認識させる必要があります。
必要な入力ファイル
プログラムを実行するには、以下のVASP計算結果ファイルが格納されたディレクトリを、それぞれ理想結晶モデル用と欠陥結晶モデル用に準備する必要があります。
理想結晶モデルのディレクトリ (例:
ideal_bulk_dir):POSCAR: 理想結晶の構造情報ファイル。OUTCAR: VASPの標準出力ファイル。原子ごとの平均コアポテンシャル情報（average of the local potential）が含まれている必要があります。DOSCAR: 状態密度 (DOS) 情報ファイル。EIGENVAL: バンドエネルギー情報ファイル。
欠陥結晶モデルのディレクトリ (例:
defect_charged_dir):POSCAR: 欠陥を含むスーパーセルの構造情報ファイル。OUTCAR: VASPの標準出力ファイル。理想結晶と同様に原子ごとの平均コアポテンシャル情報が含まれている必要があります。DOSCAR: 状態密度 (DOS) 情報ファイル。
生成される出力ファイル
プログラムの実行により、以下のファイルと出力が生成されます。
VBM_correction-out.txt:実行中の標準出力（コンソールに表示される内容）がリダイレクトされ、ログファイルとしてこのファイルに保存されます。
VBM補正の過程、読み込まれた情報、計算されたポテンシャル値、決定された\(dE_{VBM}\)などが含まれます。
VBM_correction-summary.prm:INI形式のサマリーファイルで、主要な結果が保存されます。
例:
[results] Car_dir1 = /path/to/ideal_bulk_dir Car_dir2 = /path/to/defect_charged_dir EV = 2.500000 EC = 5.000000 Eg = 2.500000 dEVBM = -0.234567 defect_atom0 = V defect_site0 = Ga
dEVBMが計算されたVBM補正量です。
グラフウィンドウ:
プログラムが終了するまで、複数のグラフが表示されます。
Vcoreの距離依存性: 欠陥中心からの距離に対して、理想結晶と欠陥結晶それぞれの\(V_{core}\)プロット。
\(\Delta V_{core}\)の距離依存性: 欠陥中心からの距離に対して、\(\Delta V_{core}\)をプロット。計算に含める/除外されるサイトが色分けされます。
\(\Delta V_{core}\)の分布: 欠陥近傍を除外したサイトの\(\Delta V_{core}\)の分布（ヒストグラムをガウス関数で平滑化したもの）が表示されます。
DOSプロット: 理想結晶と\(dE_{VBM}\)でシフトされた欠陥結晶のTotal DOSが重ねてプロットされます。バンドエッジとフェルミエネルギーも示されます。
コマンドラインでの使用例 (Usage)
基本的な使用方法は以下の通りです。
python vasp_correction_VBM.py mode CAR_dir(ideal) CAR_dir(defect) dVth WGauss Rmin_vacancy Rmin_potential Emin Emax [defect_position] [WG_DOS]
引数:
mode: \(dE_{VBM}\)を決定するためのアルゴリズム。'mode': 最頻値アルゴリズム (デフォルト)。'average': 平均値アルゴリズム。
CAR_dir(ideal): 理想結晶モデルのVASP計算結果ディレクトリのパス。CAR_dir(defect): 欠陥結晶モデルのVASP計算結果ディレクトリのパス。dVth: ポテンシャル差の閾値 (eV)。'average'モードで使用され、平均からの偏差がこの値を超えるサイトは計算から除外されます。'mode'モードでは最終的な平均計算に使用されます。WGauss: \(\Delta V_{core}\)の分布を平滑化するガウス関数のFWHM (半値全幅) (eV)。Rmin_vacancy: 欠陥を特定するための最小距離 (アンダーシュテーム)。この距離よりも原子が離れていれば空孔と見なされます。Rmin_potential: サイトポテンシャルを計算に含める最小距離 (アンダーシュテーム)。欠陥からこの距離以内にあるサイトは\(\Delta V_{core}\)の計算から除外されます。Emin: DOSプロットのエネルギー範囲の最小値 (eV)。Emax: DOSプロットのエネルギー範囲の最大値 (eV)。[defect_position]: (オプション) 欠陥中心の分数座標 (例:0.25,0.25,0.25または0.25 0.25 0.25)。指定しない場合、プログラムが自動で欠陥サイトの最初の位置を原点として設定します。[WG_DOS]: (オプション) DOSプロットの平滑化に使用するガウス関数のFWHM (半値全幅) (eV)。
コマンドラインでの具体的な使用例
ここでは、'mode'アルゴリズムを使用して、カレントディレクトリのサブディレクトリにある理想結晶 (./ideal_bulk) と欠陥結晶 (./defect_charged) の計算結果を解析する例を示します。
python vasp_correction_VBM.py mode ./ideal_bulk ./defect_charged 0.1 0.1 0.5 5.0 0.0 8.0 0.25,0.25,0.25 0.05
引数の内訳:
mode:'mode'アルゴリズムを使用。./ideal_bulk: 理想結晶モデルのディレクトリ。./defect_charged: 欠陥結晶モデルのディレクトリ。0.1:dVth(ポテンシャル差の閾値) を \(0.1\) eV に設定。0.1:WGauss(\(\Delta V_{core}\)分布の平滑化幅) を \(0.1\) eV に設定。0.5:Rmin_vacancy(空孔検出距離) を \(0.5\) Å に設定。5.0:Rmin_potential(ポテンシャル計算対象外距離) を \(5.0\) Å に設定。0.0:Emin(DOSプロット下限) を \(0.0\) eV に設定。8.0:Emax(DOSプロット上限) を \(8.0\) eV に設定。0.25,0.25,0.25: 欠陥中心位置を分数座標 \((0.25, 0.25, 0.25)\) に指定。0.05:WG_DOS(DOS平滑化幅) を \(0.05\) eV に設定。
実行結果の説明:
このコマンドを実行すると、まず標準出力に解析の進行状況、読み込まれたファイル情報、ポテンシャルの差、欠陥情報などが詳細に表示されます。同時に、これらの出力は ./defect_charged/VBM_correction-out.txt ファイルにも記録されます。
最終的に、欠陥のVBM補正量 \(dE_{VBM}\) が決定され、その値やバンドエッジ情報などが ./defect_charged/VBM_correction-summary.prm に保存されます。
そして、3種類のグラフが新しいウィンドウで表示されます。
欠陥からの距離に対する各原子サイトの\(V_{core}\)および\(\Delta V_{core}\)のプロット。
\(\Delta V_{core}\)の分布を示すプロット。
補正前の欠陥モデルDOSと補正後の理想モデルDOSが重ねて表示されるプロット。 これらのグラフウィンドウを閉じると、プログラムの実行が終了します。