vasp_defect.py ドキュメント

概要

vasp_defect.py は、VASPによる第一原理計算の結果（状態密度やバンド構造）および欠陥形成エネルギーのデータに基づき、半導体中の電子・正孔キャリア密度、欠陥密度、フェルミ準位を計算・解析するスクリプトです。

本スクリプトは、温度変化に伴うフェルミ準位の自己無撞着計算、欠陥の電荷状態による形成エネルギーの変化、電荷遷移準位の特定、およびこれらの物理量のグラフ表示機能を提供します。

依存ライブラリ

本プログラムは以下のライブラリに依存しています。

  • 標準ライブラリ: os, sys, glob, re, math

  • 非標準ライブラリ:

    • tklib （カスタムライブラリ群）

    • numpy

    • scipy

    • matplotlib

    • tkinter

入出力仕様

入力ファイル

完全結晶のVASP計算結果ファイル（CAR_path ディレクトリ内）

  • INCAR

  • POSCAR

  • CONTCAR

  • OUTCAR

  • EIGENVAL

  • DOSCAR

欠陥結晶モデルのVASP計算結果ファイル（CAR_path_defect ディレクトリ内）

  • POSCAR

欠陥形成エネルギーデータ

  • dH_path (Excelファイル。欠陥の原子、サイト、電荷、エントロピー、基準サイト数、ドープ濃度、および各ポイントでの形成エンタルピーを記録)

出力ファイル

実行モード（EF または T）に応じて、以下のファイルが出力されます（ファイル名は dH_path のベース名やパラメータに基づいて自動生成されます）。

  • テキスト形式のログ・計算結果ファイル: -out-{Point}.txt, -defect-out.txt

  • キャリア密度およびフェルミ準位の計算結果（Excel）: -N-EF-{Point}.xlsx （EF モード）または -T-{Point}.xlsx （T モード）

  • 遷移準位データ（Excel）: -dH-EF-all-{Point}.xlsx, -dH-EF-min-{Point}.xlsx

コマンドラインオプション

本スクリプトはコマンドラインから実行され、第一引数でモード（EF または T）を指定します。

EFモード（フェルミ準位依存性の解析） 引数順序:

  1. mode (EF)

  2. plot_mode (min または all)

  3. dH_path

  4. CAR_path

  5. CAR_path_defect

  6. iPoint

  7. T(electron)

  8. T(defect)

  9. EF(defect)

  10. EFmin

  11. EFmax

  12. nEF

  13. dHmin (オプション)

  14. dHmax (オプション)

  15. ignore_warning (オプション)

Tモード（温度依存性の解析） 引数順序:

  1. mode (T)

  2. dH_path

  3. CAR_path

  4. CAR_path_defect

  5. iPoint

  6. T(defect)

  7. EF(defect)

  8. Tmin

  9. Tmax

  10. nT

  11. ignore_warning (オプション)

物理的な理論とアルゴリズム（推測）

フェルミ・ディラック分布とキャリア密度

フェルミ・ディラック分布関数は、エネルギー \(E\) における電子の占有確率 \(f_e\) と正孔の占有確率 \(f_h\) を以下のように定めます。

\[f_e(E, T, E_F) = \frac{1}{\exp \left( \frac{E - E_F}{k_B T} \right) + 1}\]
\[f_h(E, T, E_F) = 1 - f_e(E, T, E_F)\]

ここで、\(k_B\) はボルツマン定数、\(T\) は電子温度 \([K]\)、\(E_F\) はフェルミ準位 \([eV]\) です。 伝導帯の電子密度 \(N_e\) \([cm^{-3}]\) および価電子帯の正孔密度 \(N_h\) \([cm^{-3}]\) は、状態密度 \(D(E)\) と分布関数の積を積分することで計算されます。

\[N_e = \int_{E_C}^{E_C + \Delta E} D(E) f_e(E, T, E_F) dE\]
\[N_h = \int_{E_V - \Delta E}^{E_V} D(E) f_h(E, T, E_F) dE\]

欠陥形成エネルギーと欠陥密度

電荷 \(q\) を持つ欠陥の形成エネルギー \(\Delta E\) は、フェルミ準位 \(E_F\) に依存して以下のように計算されます。

\[\Delta E = \Delta H_0 + q E_F\]

ここで \(\Delta H_0\) は欠陥の形成エンタルピーです。 熱平衡状態における各欠陥状態の密度 \(n\) \([cm^{-3}]\) は、分配関数 \(Z_S\) を用いてボルツマン分布に従い計算されます。

\[Z_S = \sum_{q} \exp \left( - \frac{\Delta E_q}{k_B T} \right)\]
\[n = \frac{N_0}{V} \frac{\exp \left( - \frac{\Delta E}{k_B T} \right)}{Z_S}\]

ここで \(N_0\) は対象サイトの基準数、\(V\) はセルの体積を \([cm^3]\) 単位に換算した値です。

欠陥凍結のモデル

欠陥の移動度にはエネルギー障壁が存在するため、ある温度 \(T_{def}\) 以下の温度 \(T_e\) においては、欠陥全体の濃度は \(T_{def}\) における値 \(N_{ds,Tdef}\) で凍結されると仮定されています。ただし、各欠陥間の電子的な電荷状態の再分配は電子温度 \(T_e\) で行われます。この場合、電子温度 \(T_e\) での分配関数 \(Z_{SA}\) を再計算し、以下のように電荷を帯びた欠陥密度を求めます。

\[n = N_{ds,Tdef} \frac{\exp \left( - \frac{\Delta E(T_e)}{k_B T_e} \right)}{Z_{SA}}\]

全電荷中性条件とフェルミ準位の探索

半導体全体は電気的に中性である必要があり、総電荷の偏差 \(\Delta Q\) がゼロになるようにフェルミ準位 \(E_F\) が決定されます。

\[\Delta Q(E_F) = N_h - N_e + \sum_{i} q_i n_i = 0\]

本スクリプトでは、関数 \(\Delta Q(E_F) = 0\) の解を見つけるため、最初に広い範囲で二分法（bisection method）を適用し、続いてニュートン法（Newton's method）を用いて精密な \(E_F\) を探索します。ニュートン法における微分は、微小幅 \(h\) を用いた中心差分によって計算されます。

\[\frac{\partial \Delta Q}{\partial E_F} \approx \frac{\Delta Q(E_F + h) - \Delta Q(E_F - h)}{2h}\]

クラスと関数の詳細

クラス

tkDefect

個々の欠陥サイトの特性（原子、サイト、電荷、エントロピー、基準サイト数、形成エンタルピーなど）を保持するデータクラスです。

tkDefects

複数の tkDefect オブジェクトを管理し、欠陥密度や遷移準位の計算を統括するクラスです。

  • get_names(): 登録されているユニークな欠陥名のリストを返します。

  • cal_dE(): 欠陥の形成エネルギー \(\Delta E\) を計算します。

  • calculate_total_defect_densities(): 欠陥凍結温度 \(T_{def}\) における全欠陥密度を計算します。

  • calculate_charged_defect_densities(): 凍結された全欠陥密度をもとに、電子温度 \(T_e\) における各電荷状態の欠陥密度を計算します。

  • calculate_densities(): 指定温度での電子・正孔密度および欠陥密度、総電荷を計算します。

  • dQ(): フェルミ準位 \(E_F\) に依存する全電荷の中性条件からの偏差を計算します。

  • diff(): dQ の \(E_F\) に対する数値微分を計算します。

  • find_EF(): 二分法およびニュートン法を用いて、\(\Delta Q = 0\) となるフェルミ準位 \(E_F\) を探索します。

  • find_all_transition_EF(), find_next_transition_EF(), find_order(), find_all_transitions(): 各欠陥の形成エネルギーの交点を計算し、最小の形成エンタルピーを持つ電荷状態が変化する遷移準位（Transition Level）を特定します。

  • read_excel(): 欠陥形成エネルギーのExcelファイルを読み込みます。

主要関数

usage()

コマンドライン引数の使用方法と例を標準出力に表示します。

updatevars()

sys.argv を解析し、グローバル変数（モード、ファイルパス、温度範囲、フェルミ準位範囲など）を更新します。

savecsv(), read_csv()

CSVファイルの読み書きを補助する関数群です。

DOS()

与えられたエネルギー \(E\) における状態密度を補間関数 fdos を用いて計算します。

fe(), fh()

フェルミ・ディラック分布に基づく電子および正孔の占有確率を計算します。

DOSfe(), DOSfh()

状態密度とフェルミ分布関数の積を返します。キャリア密度の積分における被積分関数です。

integrate_trapezoid(), integrate_trapezoid_by_list()

指定された関数またはデータリストを台形則（Trapezoidal rule）で数値積分します。

Ne(), Nh()

指定されたエネルギー範囲における電子密度および正孔密度を数値積分によって計算します。

FindBandEdges()

状態密度（DOS）データと閾値に基づいて、価電子帯上端 \(E_V\) と伝導帯下端 \(E_C\) を探索します。

read_files()

VASPの計算結果（POSCAR, OUTCAR, EIGENVAL, DOSCAR）を読み込み、セル体積、フェルミ準位、バンド端情報、状態密度などの物理パラメーターをグローバル変数に設定します。

check_atom_sites()

POSCAR ファイルと入力Excelファイルに記載された原子サイト数が一致するかを検証します。

exec_T()

温度依存性（T モード）の計算を実行します。各温度におけるキャリア密度、欠陥密度、およびフェルミ準位を自己無撞着に求め、Excelへの保存とグラフによる可視化を行います。

exec_EF()

フェルミ準位依存性（EF モード）の計算を実行します。指定されたフェルミ準位の範囲に対してキャリア密度や欠陥密度、総電荷の偏差を計算し、Excelへの保存とグラフによる可視化を行います。

main()

スクリプトのメインルーチン。設定の読み込み、ログの初期化を行い、指定されたモードに応じて exec_EF() または exec_T() を呼び出します。