Electron_Helmholtz_energy.py 技術ドキュメント
プログラムの動作
Electron_Helmholtz_energy.py は、VASP計算の出力ファイルである DOSCAR から電子の状態密度 (DOS: Density of States) データを読み込み、有限温度における電子の内部エネルギー、エントロピー、およびヘルムホルツ自由エネルギーの温度依存性を計算するPythonスクリプトです。また、フェルミ準位の温度変化も計算し、結果を標準出力に表形式で表示するとともに、matplotlib を用いてグラフとして可視化します。
このプログラムの主な機能は以下の通りです。
DOSCARファイルからエネルギーとDOSデータを抽出。指定された温度範囲で、フェルミ・ディラック統計に基づいて電子の物理量を計算。
温度変化に伴うフェルミ準位の変動を電荷中性条件から決定。
電子の内部エネルギー変化 (\( \Delta U_e \))、エントロピー寄与 (\( -T \Delta S_e \))、およびヘルムホルツ自由エネルギー変化 (\( \Delta F_e = \Delta U_e - T \Delta S_e \)) を計算。
計算結果を温度ごとに表形式で出力。
計算された物理量の温度依存性およびDOSプロファイルをグラフ表示。
このプログラムは、金属や半導体などの電子系が有限温度で示す熱力学的特性を理解し、VASP計算の結果をさらに深く解析するために役立ちます。
原理
本プログラムは、電子の有限温度における熱力学的特性を、状態密度 \(g(E)\) とフェルミ・ディラック分布関数 \(f(E, T, E_F)\) に基づいて計算します。
フェルミ・ディラック分布関数: 電子が特定のエネルギー \(E\) を占める確率 \(f(E, T, E_F)\) は、フェルミ準位 \(E_F\) と温度 \(T\) に依存します。 $\(f(E, T, E_F) = \frac{1}{\exp\left(\frac{(E - E_F)e}{k_B T}\right) + 1}\)\( ここで、\)e\( は電気素量、\)k_B\( はボルツマン定数です。\)T=0\( K の場合、\)E < E_F\( で \)f(E)=1\(、\)E > E_F\( で \)f(E)=0$ となります。
電子数の計算: 特定のエネルギー範囲 \(E_0\) から \(E_1\) までの電子の総数 \(N_e\) は、状態密度 \(g(E)\) とフェルミ・ディラック分布関数を用いて次のように計算されます。 $\(N_e(T, E_F) = \int_{E_0}^{E_1} g(E) f(E, T, E_F) dE\)\( プログラムでは、0 Kにおける電子数 \)N_{e0}(E_{F0})\( を基準として、有限温度 \)T\( における電子数が \)N_{e0}\( と等しくなるようにフェルミ準位 \)E_F(T)\( を決定します。これは、ニュートン法 (`scipy.optimize.newton`) を用いて方程式 \)N_e(T, E_F) - N_{e0} = 0$ を解くことで実現されます。
電子の内部エネルギー: 電子の内部エネルギー \(U_e\) は、状態密度 \(g(E)\) とフェルミ・ディラック分布関数を用いて次のように計算されます。 $\(U_e(T, E_F) = \int_{E_0}^{E_1} E \cdot g(E) f(E, T, E_F) dE\)\( プログラムの出力では、0 Kにおける内部エネルギー \)U_e(0, E_{F0})\( を基準とした変化量 \) \Delta U_e(T) = U_e(T, E_F) - U_e(0, E_{F0}) $ が表示されます。
電子のエントロピー: 電子のエントロピー \(S_e\) は、以下の式で計算されます。 $\(S_e(T, E_F) = -k_B \int_{E_0}^{E_1} g(E) [f(E, T, E_F) \ln f(E, T, E_F) + (1-f(E, T, E_F)) \ln (1-f(E, T, E_F))] dE\)\( プログラムの出力では、\)-T S_e\( の形で、エネルギー単位 (eV) に変換された値が表示されます。\)T=0$ K の場合、エントロピーはゼロとなります。
電子のヘルムホルツ自由エネルギー: 電子のヘルムホルツ自由エネルギー \(F_e\) は、内部エネルギー \(U_e\) とエントロピー \(S_e\) を用いて定義されます。 $\(F_e(T, E_F) = U_e(T, E_F) - T S_e(T, E_F)\)\( プログラムの出力では、0 Kにおけるヘルムホルツ自由エネルギーを基準とした変化量 \) \Delta F_e(T) = F_e(T, E_F) - F_e(0, E_{F0}) \( が表示されます。\) \Delta F_e(T) = \Delta U_e(T) - T \Delta S_e(T) $ とも書けます。
数値計算手法:
状態密度補間:
DOSCARから読み込んだ離散的なDOSデータは、scipy.interpolate.interp1dを用いて3次スプライン補間 (kind='cubic') され、連続的な関数f_dos(E)として扱われます。数値積分: 電子数、内部エネルギー、エントロピーの計算には、
scipy.integrate.quad関数による適応型求積法が用いられます。積分範囲は、計算負荷と精度のバランスを考慮し、\( E_F \pm 6 k_B T \) 程度の範囲に設定されます。フェルミ準位の決定:
scipy.optimize.newtonを用いたニュートン法により、電荷中性条件を満たすフェルミ準位が反復的に計算されます。
必要な非標準ライブラリとインストール方法
本プログラムの実行には、以下の非標準Pythonライブラリが必要です。
NumPy: 数値計算
SciPy: 数値積分、最適化、補間
Matplotlib: グラフ描画
これらのライブラリは、Pythonのパッケージ管理システム pip を用いてインストールできます。
pip install numpy scipy matplotlib
必要な入力ファイル
プログラムは、カレントディレクトリに存在する DOSCAR という名前のファイルを読み込みます。このファイルはVASPの出力形式に従っている必要があります。
ファイル名:
DOSCARファイル形式:
最初の6行はヘッダとしてスキップされます。
7行目以降にデータが記述されます。
各行は、エネルギー (eV)、DOS (states/unit cell/eV)、および積分DOSの順にスペースまたはタブ区切りで並んでいます。
プログラムは1列目（エネルギー）と2列目（DOS）のデータを使用します。
DOSCAR ファイルの例:
# VASP DOSCAR header (6 lines)
# Lines 1-5 contain general information.
# Line 6 specifies number of ions, # of spin channels, EMIN, EMAX, NEDOS, EFERMI, etc.
#
# EFERMI value from OUTCAR of VASP (e.g. 7.5883 eV) should be used for EF0 in the script.
#
# Data starts from line 7:
# Energy DOS IDOS
-5.000 0.100 0.100
-4.990 0.105 0.205
...
7.580 1.020 85.340 # EF0に近い値
7.590 1.050 86.390
...
15.000 0.050 120.000
プログラムは、上記の例で # VASP DOSCAR header から # Data starts from line 7: までの6行をスキップし、-5.000, 0.100 や 7.580, 1.020 のような数値を読み込みます。
生成される出力ファイル
プログラムは以下の情報を生成します。
標準出力 (コンソール): 温度、フェルミ準位、電子の内部エネルギー変化、電子のエントロピー寄与 (\( -T \Delta S_e \))、電子のヘルムホルツ自由エネルギー変化、および電子数の整合性チェック値を含む表形式のデータが出力されます。
出力例:
quad() parameters releps= 1e-05 nmaxdiv= 150 Temperature range: 0 - 1500, 50 K step DOS configuration Read from [DOSCAR] E range: -5 - 15, 0.01 eV step EF at 0K: 7.5883 eV total Ne at 0 K: 85.3458 total E at 0 K: 500.234 eV T(K) EF(eV) dE_ele(eV/u.c.) dS_ele*T(eV/u.c.) (approx) dF_ele(eV/u.c.) Ne(0 K) (check) 0 7.588300 0.000000 0.000000( 0.000000) 0.000000 85.3458(85.3458) 50 7.588299 0.000005 -0.000005( -0.000005) 0.000009 85.3458(85.3458) 100 7.588290 0.000021 -0.000020( -0.000020) 0.000041 85.3458(85.3458) ...
1500 7.587841 0.004739 -0.004739( -0.004739) 0.009478 85.3458(85.3458)
```
各列の意味は以下の通りです。
* T(K): 温度 (K)
* EF(eV): その温度におけるフェルミ準位 (eV)
* dE_ele(eV/u.c.): 0 Kからの電子内部エネルギーの変化量 (eV/unit cell)
* dS_ele*T(eV/u.c.) (approx): 電子エントロピーの寄与 (\( -T \Delta S_e \)) (eV/unit cell)。括弧内はエントロピー項の近似値。
* dF_ele(eV/u.c.): 電子ヘルムホルツ自由エネルギーの変化量 (\( \Delta F_e \)) (eV/unit cell)
* Ne(0 K) (check): 0 Kにおける電子総数（積分範囲内）。括弧内は、計算された EF での電子総数（整合性チェック用）。これらは一致するはずです。
グラフウィンドウ: プログラムの実行後、
matplotlibによってグラフウィンドウがポップアップ表示されます。このウィンドウには、以下の3つのサブプロットが含まれます。左上: 温度 (\(T\)) に対する電子の内部エネルギー変化 (\( \Delta U_e \))、エントロピー寄与 (\( -T \Delta S_e \))、およびヘルムホルツ自由エネルギー変化 (\( \Delta F_e \)) のプロット。
右上: 温度 (\(T\)) に対するフェルミ準位 (\(E_F\)) のプロット。
左下: エネルギー (\(E\)) に対するDOSプロット。0 Kのフェルミ準位 (\(E_{F0}\)) が赤色の破線で示されます。
グラフウィンドウは、ユーザーがEnterキーを押すまで表示され続けます。ファイルとして自動的に保存はされませんが、ウィンドウの機能を利用して手動で画像を保存できます。
コマンドラインでの使用例 (Usage)
本プログラムはコマンドライン引数を取らず、スクリプトを実行するだけで動作します。
python Electron_Helmholtz_energy.py
プログラムを実行する前に、必要な入力ファイル DOSCAR がスクリプトと同じディレクトリに存在することを確認してください。
コマンドラインでの具体的な使用例
ここでは、仮想的な DOSCAR ファイルが存在する状況を想定し、プログラムの実行例を示します。
まず、以下の内容を持つ DOSCAR ファイルを、Electron_Helmholtz_energy.py と同じディレクトリに作成します。
（実際にはVASPの計算結果を使用してください。ここでは例示のため簡略化しています。）
1
5
1
0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00
0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00
-5.00 15.00 2001 7.5883
-5.0000 0.0500 0.0500
-4.9900 0.0510 0.1010
-4.9800 0.0520 0.1530
...
7.5800 1.0000 85.3400
7.5900 1.0200 86.3600
...
14.9900 0.0100 120.0000
15.0000 0.0100 120.0100
（... の部分は省略されており、実際には -5.0000 から 15.0000 まで 0.01 eV刻みで2001行のデータがあるものとします。EF0 の値は 7.5883 eVに設定されています。）
次に、ターミナルまたはコマンドプロンプトを開き、Electron_Helmholtz_energy.py と DOSCAR があるディレクトリに移動します。
cd /path/to/your/calculation/directory
以下のコマンドを実行します。
python Electron_Helmholtz_energy.py
実行すると、以下のような情報がコンソールに出力され始めます。
quad() parameters
releps= 1e-05
nmaxdiv= 150
Temperature range: 0 - 1500, 50 K step
DOS configuration
Read from [DOSCAR]
E range: -5 - 15, 0.01 eV step
EF at 0K: 7.5883 eV
total Ne at 0 K: 85.3458
total E at 0 K: 500.234 eV
T(K) EF(eV) dE_ele(eV/u.c.) dS_ele*T(eV/u.c.) (approx) dF_ele(eV/u.c.) Ne(0 K) (check)
0 7.588300 0.000000 0.000000( 0.000000) 0.000000 85.3458(85.3458)
50 7.588299 0.000005 -0.000005( -0.000005) 0.000009 85.3458(85.3458)
100 7.588290 0.000021 -0.000020( -0.000020) 0.000041 85.3458(85.3458)
150 7.588277 0.000047 -0.000045( -0.000045) 0.000092 85.3458(85.3458)
200 7.588257 0.000084 -0.000080( -0.000080) 0.000164 85.3458(85.3458)
... (途中省略) ...
1450 7.587848 0.004739 -0.004739( -0.004739) 0.009478 85.3458(85.3458)
1500 7.587841 0.004739 -0.004739( -0.004739) 0.009478 85.3458(85.3458)
Press ENTER to exit>>
同時に、計算結果を示すグラフウィンドウが表示されます。ウィンドウには、左上に \( \Delta U_e, -T \Delta S_e, \Delta F_e \) の温度依存性、右上に \(E_F\) の温度依存性、左下にDOSプロットがそれぞれ表示されます。
コンソールに「Press ENTER to exit>>」と表示されたら、Enterキーを押すことでプログラムが終了し、グラフウィンドウも閉じます。