"""
VASPのDOSCARファイルを解析し、DOS、キャリア濃度、またはSeebeck係数をプロットするスクリプト。
詳細説明:
本スクリプトは、VASP計算ディレクトリからPOSCAR、DOSCAR、OUTCARファイルを読み込み、
バンド端(EVBM, ECBM)を特定します。
コマンドライン引数で指定されたモードに応じて、以下のいずれかのプロットを実行します。
- 全状態密度(DOS)と積分電子数 ('dos'モード)
- フェルミ準位に対する電子および正孔濃度 ('n'モード)
- フェルミ準位に対するSeebeck係数 ('seebeck'モード)
キャリア濃度とSeebeck係数の計算には、フェルミ・ディラック統計と定緩和時間近似(CRTA)が用いられます。
計算の数値的安定性を向上させるための工夫も含まれています。
関連リンク:
:doc:`plot_VASP_DOS_Ne_usage`
"""
import os
import sys
import argparse
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simpson
# tkvaspライブラリのインポート設定
os.environ['tkProg_Root'] = 'd:/git/tkProg'
os.environ['tklib'] = 'd:/git/tkProg/tklib/python'
os.environ['tklib'] = 'd:/git/tkProg/tklib/python' # 重複しているが既存コード変更禁止のためそのまま
sys.path.append('d:/git/tkProg/tklib/python')
from tklib.tkcrystal.tkvasp import tkVASP
[ドキュメント]
def fermi_dirac(E, Ef, T):
"""
フェルミ・ディラック分布関数を計算する。
詳細説明:
与えられたエネルギーE、フェルミエネルギーEf、温度Tに基づいてフェルミ・ディラック分布関数の値を計算する。
数値的な安定性を向上させるため、指数関数 `exp(x)` の引数 `x` が非常に大きくなることによる
オーバーフロー（`inf`）を避けるために、`arg_threshold` を用いて引数をクリッピングする。
これにより、`exp(x)` が `inf` になる領域では `1 / (inf + 1)` が `0` に収束するように振る舞う。
:param E: エネルギー (eV)。floatまたはnumpy.ndarray。
:param Ef: フェルミエネルギー (eV)。float。
:param T: 温度 (K)。float。
:returns: フェルミ・ディラック分布関数の値。floatまたはnumpy.ndarray。
"""
k_B = 8.617333262e-5 # Boltzmann定数 (eV/K)
arg = (E - Ef) / (k_B * T)
# 指数が大きすぎる場合のオーバーフローを避ける
# exp(x) -> inf の場合、1 / (inf + 1) -> 0 に収束
arg_threshold = 700.0
# NumPyのクリッピング機能を使用して引数の範囲を制限
clipped_arg = np.clip(arg, a_min=None, a_max=arg_threshold)
return 1.0 / (np.exp(clipped_arg) + 1.0)
[ドキュメント]
def calculate_concentrations(dos_info, fermi_energy_list, T):
"""
指定されたフェルミエネルギーのリストに対するキャリア濃度を計算する。
詳細説明:
DOS情報 (エネルギーと全DOS) とフェルミ・ディラック分布関数を用いて、
各フェルミエネルギーにおける電子濃度と正孔濃度を数値積分（シンプソン法）で算出する。
DOS値が負の場合は0として処理される。
:param dos_info: DOS情報を含む辞書。必須キーは'E' (エネルギー配列) と 'TotalDOS' (全DOS配列)。
例: `{'E': np.array([...]), 'TotalDOS': np.array([...])}`。
:param fermi_energy_list: 計算対象のフェルミエネルギーの配列 (eV)。numpy.ndarray。
:param T: 温度 (K)。float。
:returns:
- electron_concentrations (numpy.ndarray): 各フェルミエネルギーに対応する電子濃度の配列 (cm^-3)。
- hole_concentrations (numpy.ndarray): 各フェルミエネルギーに対応する正孔濃度の配列 (cm^-3)。
"""
energy = np.array(dos_info['E'])
total_dos = np.array(dos_info['TotalDOS'])
total_dos[total_dos < 0] = 0.0
electron_concentrations = []
hole_concentrations = []
for Ef in fermi_energy_list:
# VB内の正孔濃度 (1 - f(E))
f_holes = 1.0 - fermi_dirac(energy, Ef, T)
hole_conc = simpson(total_dos * f_holes, x=energy)
hole_concentrations.append(hole_conc)
# CB内の電子濃度 (f(E))
f_electrons = fermi_dirac(energy, Ef, T)
electron_conc = simpson(total_dos * f_electrons, x=energy)
electron_concentrations.append(electron_conc)
return np.array(electron_concentrations), np.array(hole_concentrations)
[ドキュメント]
def calculate_seebeck_coefficient(dos_info, fermi_energy_list, T):
"""
Constant Relaxation Time Approximation (CRTA) を用いてSeebeck係数を計算する。
詳細説明:
DOS情報とフェルミ・ディラック分布関数のエネルギー微分を用いて、輸送積分L0とL1を計算し、
それらの比からSeebeck係数を導出する。
数値的な安定性を確保するため、エネルギー微分の計算時に`exp(arg)`が大きくなりすぎる領域は0と見なす。
Seebeck係数に加え、各フェルミエネルギーにおける電子および正孔濃度も同時に計算し、
結果をコンソールに出力する。
:param dos_info: DOS情報を含む辞書。必須キーは'E' (エネルギー配列) と 'TotalDOS' (全DOS配列)。
例: `{'E': np.array([...]), 'TotalDOS': np.array([...])}`。
:param fermi_energy_list: 計算対象のフェルミエネルギーの配列 (eV)。numpy.ndarray。
:param T: 温度 (K)。float。
:returns: seebeck_coefficients (numpy.ndarray): 各フェルミエネルギーに対応するSeebeck係数の配列 (μV/K)。
"""
energy = np.array(dos_info['E'])
total_dos = np.array(dos_info['TotalDOS'])
total_dos[total_dos < 0] = 0.0
# フェルミ・ディラック関数のエネルギー微分
k_B = 8.617333262e-5 # Boltzmann定数 (eV/K)
seebeck_coefficients = []
print(f"{'E_F (eV)':>12} | {'Ne (cm^-3)':>12} | {'Nh (cm^-3)':>12} | {'S (uV/K)':>12}")
print("-" * 55)
for Ef in fermi_energy_list:
# オーバーフロー対策
arg = (energy - Ef) / (k_B * T)
# 指数が大きくなりすぎる領域は0と見なす
# exp(-x)はxが大きいと0に収束するため、exp(x)の引数xが小さい場合にのみ計算する
f_prime = np.zeros_like(arg)
valid_indices = arg < 500 # 経験的に安全な閾値
f_prime[valid_indices] = -1.0 / (k_B * T) * np.exp(arg[valid_indices]) / (np.exp(arg[valid_indices]) + 1)**2
# 輸送積分L_0とL_1を計算
L0 = simpson(f_prime * total_dos, x=energy)
L1 = simpson(f_prime * (energy - Ef) * total_dos, x=energy)
# Seebeck係数を計算
if L0 != 0:
seebeck = 1.0 / (k_B * T) * L1 / L0
else:
seebeck = 0.0
# 単位をV/KからμV/Kに変換
seebeck_coefficients.append(seebeck * 1.0e6)
# キャリア濃度も計算し、コンソールに出力
Ne = simpson(fermi_dirac(energy, Ef, T) * total_dos, x=energy)
Nh = simpson((1.0 - fermi_dirac(energy, Ef, T)) * total_dos, x=energy)
print(f"{Ef:12.6f} | {Ne:12.6g} | {Nh:12.6g} | {seebeck * 1.0e6:12.6f}")
return np.array(seebeck_coefficients)
[ドキュメント]
def main():
"""
スクリプトのエントリポイント。VASP DOSCARを読み込み、DOS、キャリア濃度、またはSeebeck係数をプロットする。
詳細説明:
コマンドライン引数を解析し、VASP計算ディレクトリからPOSCAR、DOSCAR、OUTCARファイルを読み込む。
読み込んだDOSデータからバンド端 (EVBM, ECBM) を特定し、VASPのフェルミ準位とバンドギャップを出力する。
選択されたモード ('dos', 'n', 'seebeck') に応じて、以下のいずれかの計算とプロットを実行する。
- 'dos'モード: 全状態密度と積分電子数を、左右のy軸を用いてプロットする。EVBM, ECBM, VASPフェルミ準位、
およびOUTCARから読み取れる場合は全電子数NELECTを垂直/水平線で表示する。
- 'n'モード: 指定された温度におけるフェルミエネルギーに対する電子および正孔濃度を対数スケールでプロットする。
- 'seebeck'モード: 指定された温度におけるフェルミエネルギーに対するSeebeck係数をプロットする。
:param args: コマンドライン引数。ArgumentParserによりパースされる。
'path': VASP計算ディレクトリへのパス。
'--mode': プロットモード ('dos', 'n', 'seebeck')。
'--T': 計算に使用する温度 (K)。
:returns: なし。グラフを表示するか、エラー時にスクリプトを終了する。
"""
parser = argparse.ArgumentParser(description='Plot DOS or carrier concentrations from VASP DOSCAR.')
parser.add_argument('path', type=str, help='Path to VASP calculation directory (e.g., ./vasp_run_dir/)')
parser.add_argument('--mode', type=str, default='dos', choices=['dos', 'n', 'seebeck'], help='Plot mode: "dos" for total DOS, "n" for carrier concentrations vs. Fermi level, "seebeck" for Seebeck coefficient.')
parser.add_argument('--T', type=float, default=300.0, help='Temperature in Kelvin for calculation (default: 300K).')
args = parser.parse_args()
# ファイルパスの定義
base_path = args.path
poscar_path = os.path.join(base_path, 'POSCAR')
doscar_path = os.path.join(base_path, 'DOSCAR')
outcar_path = os.path.join(base_path, 'OUTCAR')
if not all(os.path.exists(p) for p in [poscar_path, doscar_path]):
print(f"Error: POSCAR or DOSCAR not found in directory '{base_path}'")
sys.exit(1)
# VASPファイル読み込み
vasp_reader = tkVASP()
# POSCAR情報を読み込む (このスクリプトでは直接使用しないが、DOSCAR読み込みに必要となる場合がある)
poscar_inf = vasp_reader.read_poscar_inf(poscar_path)
# DOSCARを読み込み、エネルギーをVBMに正規化し、DOS単位を/cm3にする
doscar_inf = vasp_reader.read_doscar(doscar_path, normalize_E=True, unit='/cm3')
outcar_inf = {}
if os.path.exists(outcar_path):
# OUTCAR情報を読み込む
outcar_inf = vasp_reader.read_outcar_inf(outcar_path)
else:
print(f"Warning: OUTCAR not found in directory '{base_path}'. NELECT will be unavailable.")
if not doscar_inf or 'E' not in doscar_inf:
print("Error: Failed to read DOSCAR data.")
sys.exit(1)
# OUTCARから全電子数を取得
NELECT = outcar_inf.get('NELECT', None)
# EVBMとECBMをtkVASPのfind_band_edges_from_dos()で取得
dos_energy = doscar_inf['E']
dos_data = doscar_inf['TotalDOS']
EF_vasp = doscar_inf['Efermi'] # VASPのフェルミ準位
dos_threshold = 1.0e18 # バンド端を決定するためのDOS閾値 (states/cm^3/eV)
EVBM, ECBM = vasp_reader.find_band_edges_from_dos(dos_energy, dos_data, EF0=EF_vasp, DOSth=dos_threshold)
print(f"EVBM: {EVBM:.3f} eV")
print(f"ECBM: {ECBM:.3f} eV")
print(f"Band Gap: {ECBM - EVBM:.3f} eV")
if NELECT is not None:
print(f"NELECT: {NELECT}")
# プロットの実行
if args.mode == 'dos':
fig, ax1 = plt.subplots(figsize=(10, 6))
# DOSを左y軸にプロット
color = 'tab:blue'
ax1.set_xlabel('Energy (eV)')
ax1.set_ylabel('DOS (states/cm$^3$/eV)', color=color)
ax1.plot(dos_energy, dos_data, color=color, label='Total DOS')
ax1.fill_between(dos_energy, 0, dos_data, color='gray', alpha=0.3)
ax1.tick_params(axis='y', labelcolor=color)
# 積分電子数Nintを計算 (台形公式による累積積分)
integrated_dos = np.cumsum((dos_data[:-1] + dos_data[1:]) / 2.0 * np.diff(dos_energy))
integrated_dos = np.insert(integrated_dos, 0, 0.0) # 最初の要素は0で初期化
# 積分電子数を右y軸にプロット
ax2 = ax1.twinx()
color = 'tab:red'
ax2.set_ylabel('Integrated Electron Count (states/cm$^3$)', color=color)
ax2.plot(dos_energy, integrated_dos, color=color, linestyle='--', label='Integrated DOS')
ax2.tick_params(axis='y', labelcolor=color)
# 垂直線の追加
ax1.axvline(x=EVBM, color='blue', linestyle='--', label=f'EVBM ({EVBM:.3f} eV)')
ax1.axvline(x=ECBM, color='red', linestyle='--', label=f'ECBM ({ECBM:.3f} eV)')
ax1.axvline(x=doscar_inf['EF'], color='green', linestyle=':', label=f'Fermi Level ({doscar_inf["EF"]:.3f} eV)')
if NELECT is not None:
# 積分電子数がNELECTとなるエネルギー点を特定し、VASPフェルミ準位と比較すると良い場合もあるが、
# ここでは単にNELECTの値を水平線で示す。
ax2.axhline(y=NELECT, color='purple', linestyle='-.', label=f'NELECT ({NELECT})')
# グラフの装飾
plt.title('Total Density of States and Integrated Electron Count')
ax1.grid(True)
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')
plt.show()
elif args.mode == 'n':
# フェルミエネルギーの計算範囲を設定 (VBMから2eV下からCBMから2eV上まで)
fermi_level_range = np.linspace(EVBM - 2.0, ECBM + 2.0, 200)
e_conc, h_conc = calculate_concentrations(doscar_inf, fermi_level_range, args.T)
plt.figure(figsize=(10, 6))
plt.plot(fermi_level_range, e_conc, label='Electron Concentration ($n$)', color='red')
plt.plot(fermi_level_range, h_conc, label='Hole Concentration ($p$)', color='blue')
plt.axvline(x=(EVBM + ECBM) / 2.0, color='gray', linestyle='--', label='Midgap') # バンドギャップ中央に線を追加
# グラフの装飾
plt.title(f'Carrier Concentration vs. Fermi Level (T = {args.T} K)')
plt.xlabel('Fermi Energy ($E_F$) relative to $E_{VBM}$')
plt.ylabel('Concentration (cm$^{-3}$)')
plt.yscale('log') # キャリア濃度は対数スケールで表示
plt.grid(True, which="both", linestyle=':')
plt.legend()
plt.show()
elif args.mode == 'seebeck':
# フェルミエネルギーの計算範囲を設定 (VBMから2eV下からCBMから2eV上まで)
fermi_level_range = np.linspace(EVBM - 2.0, ECBM + 2.0, 200)
seebeck_coefficients = calculate_seebeck_coefficient(doscar_inf, fermi_level_range, args.T)
plt.figure(figsize=(10, 6))
plt.plot(fermi_level_range, seebeck_coefficients, color='green')
plt.axvline(x=EVBM, color='blue', linestyle='--', label=f'EVBM ({EVBM:.3f} eV)')
plt.axvline(x=ECBM, color='red', linestyle='--', label=f'ECBM ({ECBM:.3f} eV)')
plt.axhline(y=0.0, color='black', linestyle='-') # Seebeck係数ゼロの線
# グラフの装飾
plt.title(f'Seebeck Coefficient vs. Fermi Level (T = {args.T} K)')
plt.xlabel('Fermi Energy ($E_F$) relative to $E_{VBM}$')
plt.ylabel('Seebeck Coefficient ($μV/K$)')
plt.grid(True)
plt.legend()
plt.show()
if __name__ == "__main__":
main()