vasp_plot_epsilon.py ドキュメント

概要

vasp_plot_epsilon.py は、VASPによる計算出力ファイル（主に OUTCAR や DOSCAR）から、光学誘電関数、X線吸収スペクトル（XAS）、および状態密度（DOS）のスペクトルデータを読み取り、データの保存（Excel形式）およびグラフの描画を行うPythonスクリプトです。

依存ライブラリ

標準ライブラリに加え、以下の非標準ライブラリを使用しています。

numpy
scipy
matplotlib
tklib （ユーザー定義ライブラリ：ファイルI/O、データ処理、VASP形式のパーサなどの機能を提供）

コマンドライン引数

スクリプトの実行時には以下の引数を指定できます。

mode (位置引数、省略可、デフォルト: density) 実行モードを指定します。以下のいずれかを指定します。
- density: 密度-密度応答による誘電関数
- current: 電流-電流応答による誘電関数
- xas: コアホール効果を考慮したX線吸収スペクトル（ CH_LSPEC=.TRUE. での出力）
- dos: 全状態密度（ DOSCAR から読み取り）
CAR_dir (位置引数、省略可、デフォルト: .) 対象となるVASPの出力ファイルが存在するディレクトリのパスを指定します。
Emin (位置引数、省略可、デフォルト: None) プロットおよびデータ抽出を行うエネルギーの最小値 \(eV\) を指定します。
Emax (位置引数、省略可、デフォルト: None) プロットおよびデータ抽出を行うエネルギーの最大値 \(eV\) を指定します。
--sigma (オプション、デフォルト: 0.0) XASおよびDOSのポストプロセスとして適用する、ガウシアンブロードニングの標準偏差（広がり幅） \(eV\) を指定します。

入出力仕様

入力ファイル

OUTCAR mode が density, current, xas の場合に読み込まれます。誘電関数の実部・虚部のテンソル成分（XX, YY, ZZ, XY, YZ, ZX）が抽出されます。
DOSCAR mode が dos の場合に読み込まれます。ヘッダからフェルミエネルギーを読み取り、全状態密度のデータを取得します。
INCAR, POSCAR, CONTCAR tklib を通じてパスの存在確認や付随する情報の取得が行われます（直接データ抽出に使用されるかはコードからは一部確認できません）。

出力ファイル

実行された mode に応じて、以下のExcelファイルおよび画像ファイルが作業ディレクトリに出力されます。

density または current モード
- epsilon.xlsx: 各エネルギー点での誘電関数の実部 (e1) と虚部 (e2) の6つのテンソル成分を保存します。
xas モード
- xas.xlsx: エネルギー、生のXASテンソル成分、ブロードニング適用後のテンソル成分、等方的吸収強度 (xas_isotropic)、および対角成分の和 (xas_diagonal_sum) を保存します。
- xas.png: XASスペクトルのプロット画像。
dos モード
- dos.xlsx: フェルミエネルギーからの相対エネルギー、生の状態密度（スピン分極がある場合はアップとダウン）、およびブロードニング適用後の状態密度を保存します。
- dos.png: DOSスペクトルのプロット画像。

理論的背景と計算アルゴリズム（推測）

本スクリプトは、固体物理学において物質の光学特性や電子構造を評価するためのスペクトル解析を行います。

誘電関数（密度-密度および電流-電流）

VASPの線形応答計算（独立粒子近似）により得られるマクロな複素誘電関数テンソル \(\epsilon(\omega)\) を処理します。複素誘電関数は実部 \(\epsilon_1(\omega)\) と虚部 \(\epsilon_2(\omega)\) から成ります。

\[\epsilon(\omega) = \epsilon_1(\omega) + i\epsilon_2(\omega)\]

スクリプトは OUTCAR ファイルから以下のキーワードを検索してデータを抽出します。

IMAGINARY DIELECTRIC FUNCTION.+density または current
REAL DIELECTRIC FUNCTION.+density または current

これらのデータはテンソルの独立な6成分（XX, YY, ZZ, XY, YZ, ZX）として保持されます。

X線吸収スペクトル (XAS)

xas モードは、コアホールを導入した計算によって得られるX線吸収スペクトルを扱います。このとき、XASの強度は誘電関数テンソルの虚部 \(\epsilon_2(\omega)\) に比例するため、実質的に虚部を読み取っています。異方性を持つ物質の等方的なXAS強度 \(\sigma_{iso}\) は、テンソルの対角成分の平均として計算されます（推測）。

\[\sigma_{iso} = \frac{1}{3} (\sigma_{xx} + \sigma_{yy} + \sigma_{zz})\]

コード内では変数 isotropic がこの数式通りに計算されており、出力として保存されます。

状態密度 (DOS)

DOSCAR から全状態密度 \(D(E)\) を読み取ります。エネルギー軸は、ファイルから読み取ったフェルミエネルギー \(E_F\) を基準として再計算されます。

\[E_{plot} = E_{raw} - E_F\]

スピン分極計算の場合、アップスピンとダウンスピンのそれぞれについて状態密度を取得し、ダウンスピンの状態密度を負の値としてプロットします（推測）。

ガウシアンブロードニング (Gaussian Broadening)

XASおよびDOSデータに対して、スペクトルの微細な振動を平滑化し、実験的な分解能や寿命による広がりをシミュレートするためにガウシアンフィルタを適用します（推測）。関数 gaussian_broaden にて実装されています。ブロードニング後のスペクトル \(I_{broad}(E)\) は、元のスペクトル \(I(E)\) と標準偏差 \(\sigma\) のガウス関数との畳み込み積分によって得られます。

\[I_{broad}(E) = \int I(E') \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(E-E')^2}{2\sigma^2}\right) dE'\]

数値計算上では、エネルギーグリッドのステップ幅 \(\Delta E\) を考慮して、分散が \(\sigma / \Delta E\) となるように scipy.ndimage.gaussian_filter1d を適用しています。また、エネルギー間隔が不均一な場合は、一時的に均一なグリッドに線形補間を行ってから畳み込みを実施し、元のグリッドに逆補間するアルゴリズムが採用されています。