`vasp_correction_charge.py` 技術ドキュメント

プログラムの動作

vasp_correction_charge.py は、VASP計算で得られた結果から、周期境界条件下の帯電欠陥計算に不可欠な有限サイズ補正（帯電欠陥のエネルギー補正および双極子相互作用の補正）を推定するためのPythonスクリプトです。

主な機能:

双極子補正情報の読み取り: IDIPOL=4 を設定したVASP計算のOUTCARファイルから、双極子および四重極子によるエネルギー補正に関する情報を読み取ります。特に、電荷系に対するエネルギー補正 Ecorr(charge) を抽出します。
誘電率の読み取り:
- 電子誘電率: LOPTICS=.TRUE. を設定したVASP計算のOUTCARファイルから、電子的な誘電率テンソルの各成分を抽出し、その平均値を計算します。数値として直接入力することも可能です。
- イオン誘電率: LEPSILON=.TRUE. を設定したVASP計算のOUTCARファイルから、格子振動によるイオン的な誘電率テンソルの各成分を抽出し、その平均値を計算します。数値として直接入力することも可能です。
欠陥電荷の計算: 欠陥を含む系のPOTCAR、POSCAR、OUTCARファイルから、価電子の合計数とVASPが計算した全電子数 NELECT を基に、欠陥の電荷 $q$ を計算します。
有限サイズ補正の計算: 読み取った Ecorr(charge)、計算された電荷 $q$、および総誘電率 $\epsilon_{total}$ を用いて、帯電欠陥に対する最終的なエネルギー補正 dEcorr を計算します。
結果の出力: 計算された誘電率、電荷、補正エネルギーなどのサマリー情報をINI形式のファイルとして保存し、同時に標準出力にも表示します。

解決する課題: VASPのような周期境界条件を用いた第一原理計算では、帯電した欠陥を扱う場合、周期的なセル間に存在する電荷間のクーロン相互作用によって、計算される形成エネルギーがセルサイズに依存してしまいます。本プログラムは、Makov-Payneの補正モデルなどを参考に、この有限サイズ効果によるエネルギー誤差を見積もり、補正することを目的としています。

原理

本プログラムは、VASPのOUTCARファイルから直接得られる情報と、結晶の誘電率を用いて、帯電欠陥の形成エネルギーに対する有限サイズ補正を計算します。

電荷 $q$ の決定: 系の全価電子数 $n_{Ve}$ は、POSCARに記述された原子の種類と数、およびPOTCARに記述された各原子の価電子数から決定されます。VASPのOUTCARから読み取られる全電子数 NELECT との差分から、系の全電荷 $q$ が以下の式で計算されます。 $$q = n_{Ve} - \text{NELECT}$$
誘電率の決定:
- 電子誘電率 $\epsilon_{el}$: LOPTICS=.TRUE. を設定したVASP計算のOUTCARファイルから、REAL DIELECTRIC FUNCTION のエネルギー $E=0$ における実部テンソルの対角成分の平均値として取得されます。
- イオン誘電率 $\epsilon_{ion}$: LEPSILON=.TRUE. を設定したVASP計算のOUTCARファイルから、MACROSCOPIC STATIC DIELECTRIC TENSOR IONIC の対角成分の平均値として取得されます。
- 総誘電率 $\epsilon_{total}$: プログラムでは、電子誘電率の平均値 eps_av_el とイオン誘電率の平均値 epsionav の和として計算されます。 $$\epsilon_{total} = \epsilon_{el,av} + \epsilon_{ion,av}$$ ここで、$av$ はテンソルの対角成分の平均を表します。
VASPの Ecorr(charge) の利用と補正エネルギーの計算: VASPの IDIPOL=4 を設定した計算では、OUTCARに energy correction for charged system として Ecorr(charge) が出力されます。この値は、G. Makov and M. C. Payne のモデルに基づいた有限サイズ補正の一部であり、欠陥電荷 $q$ とスーパーセルの線形サイズ $L$（体積の立方根）を用いて、以下のような形式で表されます。 $$E_{corr(charge)} = q^2 \cdot \frac{aM}{L}$$ ここで、$aM$ はマデルング定数に相当する係数で、プログラムはこの式を基に $aM$ を逆算します。 $$aM = \frac{E_{corr(charge)} \cdot L}{q^2}$$ 最終的な補正エネルギー `dEcorr` は、この $aM$ と計算された総誘電率 $\epsilon_{total}$ を用いて、以下の式で計算されます。 $$dE_{corr} = \frac{1}{3} \cdot \frac{q^2 \cdot aM}{L \cdot \epsilon_{total}}$$ この補正式は、一般的なMakove-Payne補正 ($E_{corr} = -\frac{q^2 \alpha}{2\epsilon L}$) に類似していますが、コードに記述された具体的な係数 ($1/3$) と誘電率 $\epsilon_{total}$ の定義に基づいて計算されます。

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

このプログラムは、標準ライブラリの他に以下の非標準ライブラリに依存しています。

numpy: 数値計算のためのライブラリ。
scipy: 科学技術計算のためのライブラリ。scipy.interpolate.interp1d が使用されています。
matplotlib: グラフ描画のためのライブラリ。matplotlib.pyplot が使用されていますが、本プログラムの主要な機能では直接グラフ描画は行われません。
tklib: このプログラムは、tklib というカスタムライブラリに強く依存しています。このライブラリには、ファイル操作 (tkfile, tkinifile), 正規表現 (tkre), VASP関連ファイルパーシング (tkvasp), およびその他のユーティリティ関数が含まれています。

インストール方法:

numpy, scipy, matplotlib は pip コマンドでインストールできます。

pip install numpy scipy matplotlib

tklib は、このプログラムの実行環境に別途提供されるカスタムライブラリです。通常は、tklib ディレクトリをPythonのサイトパッケージパス、またはスクリプトが実行されるディレクトリからアクセス可能な場所に配置する必要があります。具体的なインストール方法は、tklib の提供元に確認してください。

必要な入力ファイル

本プログラムは、以下のVASP計算結果ディレクトリ、またはその中の特定のファイル、あるいは直接数値を入力として受け取ります。

CAR_IDIPOL_dir (IDIPOL=4 計算ディレクトリ):
- 内容: IDIPOL=4 を設定して実行されたVASP計算のディレクトリ、またはその中の OUTCAR ファイル。
- 期待されるデータ: OUTCAR ファイルには、DIPCOR: dipole corrections for dipol セクションが含まれており、特に energy correction for charged system の値が読み取られます。
- 例: ./defect_IDIPOL4/ または ./defect_IDIPOL4/OUTCAR
CAR_DIEL_el_dir (電子誘電率計算ディレクトリまたは値):
- 内容: LOPTICS=.TRUE. を設定して実行されたVASP計算のディレクトリ、またはその中の OUTCAR ファイル。あるいは、直接電子誘電率の平均値（浮動小数点数）。
- 期待されるデータ: OUTCAR ファイルには、REAL DIELECTRIC FUNCTION セクションが含まれており、エネルギーが 0.000000 eV の行から誘電率の実部テンソルの対角成分が読み取られます。
- 例: ./diel_el_LOPTICS/ または ./diel_el_LOPTICS/OUTCAR または 5.0
CAR_DIEL_ion_dir (イオン誘電率計算ディレクトリまたは値):
- 内容: LEPSILON=.TRUE. を設定して実行されたVASP計算のディレクトリ、またはその中の OUTCAR ファイル。あるいは、直接イオン誘電率の平均値（浮動小数点数）。
- 期待されるデータ: OUTCAR ファイルには、MACROSCOPIC STATIC DIELECTRIC TENSOR (incl. local field effects) と MACROSCOPIC STATIC DIELECTRIC TENSOR IONIC のセクションが含まれており、それぞれから誘電率テンソルの対角成分が読み取られます。
- 例: ./diel_ion_LEPSILON/ または ./diel_ion_LEPSILON/OUTCAR または 2.5
CAR_charge_dir (欠陥電荷計算ディレクトリ):
- 内容: 欠陥を含む系のVASP計算のディレクトリ、またはその中の OUTCAR, POSCAR, POTCAR ファイル。
- 期待されるデータ:
  - OUTCAR: NELECT (全電子数) の値。
  - POSCAR: 結晶構造、セル体積を決定するための情報（格子ベクトル、原子位置）。
  - POTCAR: 各原子種の価電子数。
- 例: ./defect_charged/ または ./defect_charged/OUTCAR (POSCAR, POTCARは同じディレクトリ内に存在することが期待される)

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

プログラムは実行時に以下のファイルを生成します。

charge_correction-out.txt:
- ファイル名: CAR_IDIPOL_dir として指定されたディレクトリ内に生成されます。例: ./defect_IDIPOL4/charge_correction-out.txt
- 内容: プログラムの実行中に標準出力に表示されるすべての情報（進捗メッセージ、読み取られたVASP情報、計算結果など）がリダイレクトされて保存されるログファイルです。

charge_correction-summary.prm:

ファイル名: CAR_IDIPOL_dir として指定されたディレクトリ内に生成されます。例: ./defect_IDIPOL4/charge_correction-summary.prm
内容: 計算された誘電率、電荷、補正エネルギーなどの主要な結果をINI形式でまとめた要約ファイルです。

形式:

[results]
Car_dir1 = <CAR_IDIPOL_dirの絶対パス>
Car_dir2 = <CAR_DIEL_el_dirの絶対パス>
Car_dir3 = <CAR_DIEL_ion_dirの絶対パス>
Car_dir4 = <CAR_charge_dirの絶対パス>
eps_av_el = <電子誘電率の平均値 (float)>
epsionav = <イオン誘電率の平均値 (float)>
eps_tot = <総誘電率の平均値 (float)>
q = <欠陥電荷 (float)>
L = <セルサイズL (float, Å)>
Error_charge = <OUTCARから読み取られたEcorr(charge) (float, eV)>
dEcorr = <計算された最終補正エネルギー (float, eV)>

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

プログラムの実行形式は以下の2種類があります。

誘電率をOUTCARファイルから読み取る場合:

python vasp_correction_charge.py CAR_IDIPOL_dir CAR_DIEL_el_dir CAR_DIEL_ion_dir CAR_charge_dir

電子誘電率を数値で直接指定する場合:
```
python vasp_correction_charge.py CAR_IDIPOL_dir CAR_DIEL_el_dir_value CAR_DIEL_ion_dir CAR_charge_dir
```
ここで、CAR_DIEL_el_dir_value は電子誘電率の数値（例: 5.0）です。

引数の説明:

CAR_IDIPOL_dir: IDIPOL=4 で計算されたVASPの出力ディレクトリまたはOUTCARファイル。
CAR_DIEL_el_dir / CAR_DIEL_el_dir_value: LOPTICS=.TRUE. で計算されたVASPの出力ディレクトリまたはOUTCARファイル、または電子誘電率の数値。
CAR_DIEL_ion_dir: LEPSILON=.TRUE. で計算されたVASPの出力ディレクトリまたはOUTCARファイル。
CAR_charge_dir: 欠陥電荷を計算するVASPの出力ディレクトリ（OUTCAR, POSCAR, POTCAR を含む）。

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

以下に、具体的なディレクトリ名と数値を指定した実行例を示します。

シナリオ:

IDIPOL=4 計算結果: ./data/defect_IDIPOL4
電子誘電率 (LOPTICS=.TRUE.) 計算結果: ./data/diel_el_optics
イオン誘電率 (LEPSILON=.TRUE.) 計算結果: ./data/diel_ion_epsilon
欠陥電荷計算結果: ./data/defect_charged

例1: 全ての誘電率をVASPのOUTCARから読み取る場合

python vasp_correction_charge.py ./data/defect_IDIPOL4 ./data/diel_el_optics ./data/diel_ion_epsilon ./data/defect_charged

実行結果の説明:

./data/defect_IDIPOL4/charge_correction-out.txt に、実行ログが保存されます。

./data/defect_IDIPOL4/charge_correction-summary.prm に、以下のような内容の要約ファイルが生成されます（数値は仮のものです）。

[results]
Car_dir1 = /path/to/data/defect_IDIPOL4
Car_dir2 = /path/to/data/diel_el_optics
Car_dir3 = /path/to/data/diel_ion_epsilon
Car_dir4 = /path/to/data/defect_charged
eps_av_el = 4.500000
epsionav = 2.000000
eps_tot = 6.500000
q = -1.000000
L = 10.000000
Error_charge = -0.500000
dEcorr = -0.025641

標準出力にも同様の計算過程と結果が出力されます。

例2: 電子誘電率を数値で直接指定する場合 （電子誘電率が実験値などで $5.0$ と既知の場合）

python vasp_correction_charge.py ./data/defect_IDIPOL4 5.0 ./data/diel_ion_epsilon ./data/defect_charged