VASP DOSおよびキャリア濃度、Seebeck係数プロットツール

プログラムの動作

plot_VASP_DOS_Ne.py は、VASP (Vienna Ab-initio Simulation Package) の計算結果から得られる状態密度 (DOS) データを用いて、様々な電子輸送特性を解析し可視化するためのPythonスクリプトです。

このプログラムの主な目的は以下の通りです。

  1. 状態密度 (DOS) のプロット: VASPの DOSCAR ファイルから全状態密度 (Total DOS) を読み込み、エネルギーに対するDOSと積算電子数をプロットします。価電子帯上端 (EVBM)、伝導帯下端 (ECBM)、バンドギャップ、フェルミ準位、および全電子数 (NELECT) がグラフ上に示されます。

  2. キャリア濃度の計算とプロット: 特定の温度において、フェルミ準位の位置を変化させた場合の電子濃度 (\(n\)) と正孔濃度 (\(p\)) を計算し、フェルミエネルギーの関数としてプロットします。これにより、材料がN型またはP型半導体として振る舞う際のキャリア濃度変化を評価できます。

  3. Seebeck係数の計算とプロット: 定緩和時間近似 (Constant Relaxation Time Approximation, CRTA) を用いて、特定の温度におけるSeebeck係数をフェルミ準位の関数として計算し、プロットします。これにより、熱電材料の性能評価に必要なSeebeck係数とドーピングレベルの関係を調べることができます。

このプログラムは、VASP計算の出力から材料の電子構造と熱電特性を迅速に評価したい研究者やエンジニアにとって有用です。

原理

このプログラムは、主に半導体のフェルミ統計と熱電輸送理論に基づいています。

1. フェルミ・ディラック分布関数

電子が特定のエネルギー準位 \(E\) に存在する確率 \(f(E)\) は、フェルミ・ディラック分布関数によって記述されます。

\[f(E) = \frac{1}{e^{(E - E_f)/(k_B T)} + 1}\]

ここで、\(E_f\) はフェルミエネルギー、\(k_B\) はボルツマン定数、\(T\) は絶対温度です。

2. キャリア濃度

電子濃度 \(n\) と正孔濃度 \(p\) は、状態密度 \(D(E)\) とフェルミ・ディラック分布関数を用いて、次のように計算されます。

  • 電子濃度: 伝導帯に存在する電子の数を示します。 $\(n = \int D(E) f(E) dE\)$

  • 正孔濃度: 価電子帯に存在する正孔の数を示します。 $\(p = \int D(E) (1 - f(E)) dE\)$

積分は、DOSデータが与えられたエネルギー範囲全体にわたって実行されます。

3. Seebeck係数 (定緩和時間近似, CRTA)

Seebeck係数 \(S\) は、材料の熱電性能を示す重要な指標の一つであり、温度勾配によって生じる電圧を記述します。プログラムでは定緩和時間近似 (Constant Relaxation Time Approximation, CRTA) を用いて計算されます。CRTAでは、キャリアの緩和時間 \(\tau\) がエネルギーに依存しない定数であると仮定します。

Seebeck係数は以下の輸送積分 \(L_0\) と \(L_1\) を用いて計算されます。

  • 輸送積分 \(L_0\): $\(L_0 = \int (-\frac{\partial f}{\partial E}) D(E) dE\)$

  • 輸送積分 \(L_1\): $\(L_1 = \int (-\frac{\partial f}{\partial E}) (E - E_f) D(E) dE\)$

ここで、\(D(E)\) は状態密度、\(f\) はフェルミ・ディラック分布関数です。関数 \(-\frac{\partial f}{\partial E}\) は、フェルミ準位近傍でピークを持つデルタ関数に近い形状をしており、輸送現象に寄与するエネルギー範囲を特定します。

Seebeck係数 \(S\) は、これらの輸送積分を用いて次のように定義されます。

\[S = \frac{1}{k_B T} \frac{L_1}{L_0}\]

このプログラムでは、計算されたSeebeck係数は \(μV/K\) 単位に変換されてプロットされます。

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

このプログラムの実行には、以下の非標準Pythonライブラリが必要です。

  • NumPy: 数値計算（配列操作、数学関数）

  • Matplotlib: グラフ描画

  • SciPy: 科学技術計算（数値積分 simpson 関数）

  • tkvasp (tkVASP): VASPファイルの読み込みと解析（このプログラムのローカル環境に特化したライブラリ）

NumPy, Matplotlib, SciPy は pip を使用してインストールできます。

pip install numpy matplotlib scipy

tkvasp は、プログラムコード内で特定のローカルパス (d:/git/tkProg/tklib/python) を参照するように設定されています。このライブラリは標準の pip ではインストールできず、ユーザーが d:/git/tkProg のパスに tkProg リポジトリをクローンし、Python環境からアクセスできるようにする必要があります。

プログラムは以下の環境変数を設定して tkvasp ライブラリを見つけます。

os.environ['tkProg_Root'] = 'd:/git/tkProg'
os.environ['tklib'] = 'd:/git/tkProg/tklib/python'
sys.path.append('d:/git/tkProg/tklib/python')

したがって、このパスに tkProg リポジトリが存在し、tklib/python ディレクトリ内に tklib.tkcrystal.tkvasp モジュールが配置されている必要があります。

必要な入力ファイル

プログラムは、VASP計算が実行されたディレクトリを指定し、その中にある以下のファイルを読み込みます。

  1. POSCAR: 結晶構造情報が含まれるファイル。tkVASP ライブラリがセルの体積などの情報を読み込むために使用します。

  2. DOSCAR: VASPのDOS計算結果が含まれるファイル。エネルギー、Total DOS、フェルミエネルギー (EFERMI) などの情報が含まれます。このファイルはプログラムの主要なデータソースです。

  3. OUTCAR (オプション): VASP計算の出力ファイル。全電子数 (NELECT) などの追加情報を読み込むために使用されます。このファイルが存在しない場合でもプログラムは実行されますが、NELECT情報は利用できません。

これらのファイルは、指定されたVASP計算ディレクトリ内に配置されている必要があります。

生成される出力ファイル

plot_VASP_DOS_Ne.py は、ファイルを直接保存する機能は持っていません。 代わりに、計算結果をグラフィカルなプロットとして画面上に表示します。

  • --mode dos: エネルギー対Total DOSと積算電子数のグラフが表示されます。

  • --mode n: フェルミエネルギー対電子濃度および正孔濃度のグラフが表示されます。濃度は対数スケールで表示されます。

  • --mode seebeck: フェルミエネルギー対Seebeck係数のグラフが表示されます。

また、--mode seebeck を実行した場合、各フェルミエネルギーにおける電子濃度、正孔濃度、Seebeck係数の値がコンソールに表形式で出力されます。

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

plot_VASP_DOS_Ne.py は、コマンドライン引数を受け取ります。

python plot_VASP_DOS_Ne.py <path_to_vasp_directory> [--mode <mode>] [--T <temperature>]

  • <path_to_vasp_directory>: 必須引数。VASP計算結果 (POSCAR, DOSCAR, OUTCAR) が格納されているディレクトリのパスを指定します。

  • --mode <mode>: オプション引数。プロットの種類を指定します。

    • dos: Total DOSと積算電子数をプロットします。（デフォルト）

    • n: キャリア濃度 (電子・正孔) をフェルミレベルの関数としてプロットします。

    • seebeck: Seebeck係数をフェルミレベルの関数としてプロットします。

  • --T <temperature>: オプション引数。キャリア濃度やSeebeck係数計算に使用する温度をケルビン単位で指定します。デフォルトは 300.0 Kです。

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

ここでは、./vasp_run_dir というディレクトリにVASP計算結果ファイルが存在することを仮定します。

1. Total DOSと積算電子数のプロット

VASP計算ディレクトリ内の DOSCAR から状態密度と積算電子数をプロットします。

python plot_VASP_DOS_Ne.py ./vasp_run_dir --mode dos

実行結果の説明: 新しいMatplotlibウィンドウが開き、エネルギー (eV) を横軸に、Total DOS (states/cm\(^3\)/eV) を左y軸に、積算電子数 (states/cm\(^3\)) を右y軸にしたグラフが表示されます。価電子帯上端 (EVBM)、伝導帯下端 (ECBM)、VASPのフェルミ準位、および OUTCAR から読み取られた場合はNELECTが垂直線や水平線として表示され、バンド構造に関する情報が視覚的に提供されます。

2. キャリア濃度のプロット

フェルミエネルギーに対する電子濃度と正孔濃度を300Kでプロットします。

python plot_VASP_DOS_Ne.py ./vasp_run_dir --mode n --T 300

実行結果の説明: 新しいMatplotlibウィンドウが開き、フェルミエネルギー (eV) を横軸に、濃度 (cm\(^{-3}\)) を対数スケールにした縦軸に、電子濃度 (\(n\)) と正孔濃度 (\(p\)) の曲線が表示されます。中性点（Midgap）が垂直線で示され、フェルミ準位がどこにあるときにN型またはP型半導体として振る舞うかの関係が示されます。

3. Seebeck係数のプロットとコンソール出力

フェルミエネルギーに対するSeebeck係数を500Kでプロットし、詳細な計算結果をコンソールにも出力します。

python plot_VASP_DOS_Ne.py ./vasp_run_dir --mode seebeck --T 500

実行結果の説明: まず、以下のような表形式のデータがコンソールに出力されます。

 E_F (eV) | Ne (cm^-3) | Nh (cm^-3) |    S (uV/K)
-------------------------------------------------------
  -2.000000 | 1.00000e-06 | 1.00000e+20 |  123.456789
  -1.980000 | 1.00000e-06 | 9.00000e+19 |  120.123456
  ...
   0.000000 | 1.00000e+18 | 1.00000e+18 |    0.000000
  ...
   2.000000 | 1.00000e+20 | 1.00000e-06 | -123.456789

その後、新しいMatplotlibウィンドウが開き、フェルミエネルギー (eV) を横軸に、Seebeck係数 (\(μV/K\)) を縦軸にしたグラフが表示されます。EVBMとECBMが垂直線で示され、Seebeck係数のエネルギー依存性が視覚的に評価できます。正の値は正孔による支配、負の値は電子による支配を示します。