vasp_defect.py ダウンロード/コピー

vasp_defect.py をダウンロード

vasp_defect.py

vasp_defect.py

"""
概要:
    VASPのDOSCARファイルに基づきキャリア密度とフェルミ準位を計算するスクリプト。
詳細説明:
    このスクリプトは、VASP計算から得られたDOSCAR、OUTCAR、EIGENVAL、POSCARファイル、
    および欠陥形成エネルギーのExcelファイルを利用して、半導体の電子・正孔キャリア密度、
    欠陥密度、およびフェルミ準位の温度依存性またはフェルミ準位依存性を計算・解析します。
    特に、温度変化に伴うフェルミ準位の自己無撞着計算や、欠陥の電荷状態による形成エネルギーの変化、
    遷移準位の特定、およびこれらの特性のグラフ表示をサポートします。
関連リンク:
    vasp_defect_usage
"""
import os
import sys
import glob
import re
from math import exp, sqrt, log

try:
    import tklib.tkimport as imp
except Exception as e:
    print()
    print("######################################################################")
    print("###########  ERROR ERROR ERROR ERROR ERROR ERROR #####################")
    print("######################################################################")
    print(f"# Failed to import [tklib.tkimport] module ({e}).")
    print(f"#  Add [tkProg]{os.sep}tklib{os.sep}python to PYTHONPATH variable")
    print(f"#  Current PYTHONPATH:", sys.path)
    print("######################################################################")
    input("Press ENTER to terminate>>\n")
    exit()

np    = imp.import_lib("numpy", stop_by_error = False)
scipy = imp.import_lib("scipy", stop_by_error = False)
mpl   = imp.import_lib("matplotlib", stop_by_error = False)
tk    = imp.import_lib("tkinter", stop_by_error = False)

from numpy import arange
from scipy import integrate         # 数値積分関数 integrateを読み込む
from scipy import optimize          # newton関数はscipy.optimizeモジュールに入っている
from scipy.interpolate import interp1d
from matplotlib import pyplot as plt

from tklib.tkfile import tkFile
import tklib.tkre as tkre
from tklib.tkutils import mprint, IsDir, IsFile, SplitFilePath, modify_path, delete_file, minmax_xy
from tklib.tkutils import terminate, safe_getelement, pint, pfloat, getarg, getintarg, getfloatarg
from tklib.tkutils import colors
from tklib.tkapplication import tkApplication
from tklib.tkparams import tkParams
from tklib.tksci.tksci import Reduce01, Round
from tklib.tksci.tkmatrix import make_matrix1, make_matrix2, make_matrix3
from tklib.tkcrystal.tkcif import tkCIF, tkCIFData
from tklib.tkcrystal.tkcrystal import tkCrystal
from tklib.tkcrystal.tkvasp import tkVASP
from tklib.tksci import tkequation
import tklib.tkcsv as csv # csvモジュールを使用するためにtkcsvをcsvとしてimport
from tklib.tkexcel import tkExcel
from tklib.tkgraphic.tkplotevent import tkPlotEvent
from tklib.tkgui.tksimple_gui import get_window_from_plt, CustomDialog_with_config


# constants
pi   = 3.14159265358979323846
h    = 6.6260755e-34    # Js";
hbar = 1.05459e-34      # "Js";
c    = 2.99792458e8     # m/s";
e    = 1.60218e-19      # C";
kB   = 1.380658e-23     # JK<sup>-1</sup>";
me   = 9.1093897e-31    # kg";
log10 = log(10.0)

#==================================
# global variables
#==================================
Debug = 0

# option flag to stop ISPIN=2 visualization
stop_spin_polarized = False

#mode: 'EF', 'T'
mode = 'EF'
#mode = 'T'

# plot_mode: 'all', 'min'
plot_mode = 'all'
plot_mode = 'min'

# index of point to be plotted
iPoint = 0

# files
outputfile  = None
outfp       = None
dH_path     = 'input.xlsx'
CAR_path    = '.'
DOSCAR_path = None
CAR_path_defect = None

T0    = 300      # K, Electron tempeature
Tdef  = 300      # K, Defect frozen tempearture
EFdef = 'eq'     # eV, EF at Tdef to freeze defect densities. 'eq' will calculate EF,eq at Tdef

# EF range for mode = 'EF'
dEFmin = -0.2  # measured from EV, eV
dEFmax =  0.2  # measured from EC, eV
nEF    = 101

# Temperature range for mode = 'T'
Tmin  = 300  # K
Tmax  = 600
nT    = 11
Tstep = None

# Plot configuration
view_Emin  = -5.0    # eV
view_Emax  = 10.0
view_dHmax =  5.0
view_Nmin  =  1.0e5  # cm-3
colors = ['#000000', '#ff0000', '#00aa00', '#0000ff', '#aaaa00', '#ff00ff', '#00ffff', '#aa0000', '#00aa00', '#0000aa']
#colors = ['k', 'r', 'g', 'b', 'y', 'm', 'c']
fontsize = 16
legend_fontsize = 8

# integration parameters
# integrator: 'interpolate', 'raw'
integrator = 'interpolate'
#integrator = 'raw'

# Integration range w.r.t. kBT
Einteg0 = -3.0
Einteg1 =  3.0
nrange  =  6.0
iEinteg0 = None
iEinteg1 = None

# Bisection method parameters
# Initial search range
#eps_bisec          = 1.0e-3
#nmaxiter_bisection = 300
eps_bisec          = 1.0e-5
nmaxiter_bisection = 100

# Newton method parmeters
dump_newton     = 0.5
h_newton        = 1.0e-8
nmaxiter_newton = 0
eps_newton      = 1.0e-5

iprintinterval = 1


# Global variables
Vcell      = None  # A^3
E_raw      = None
dos_raw    = None
dos_norm   = None
nDOS       = None

E_norm     = None
Enorm_min  = None
Enorm_max  = None
Enorm_step = None

Emin    = None
Emax    = None
Estep   = None

EV      = None
EC      = None

EFmin = None
EFmax = None
EFstep = None

dHmin = None
dHmax = None

# exponent threshold
nexp = 100.0

ignore_warning = False


app = tkApplication()


#==================================
# Definition of tkDefect class
#==================================
class tkDefect():
    """
    概要:
        個々の欠陥サイトの特性を保持するクラス。
    詳細説明:
        欠陥の原子ラベル、サイトラベル、電荷、エントロピー、基準サイト数、
        ドープ濃度、および複数の形成エンタルピー値を管理します。
    """
    def __init__(self, **args):
        """
        概要:
            tkDefectクラスのコンストラクタ。
        詳細説明:
            欠陥の各種プロパティを初期化します。
        引数:
            :param args: 欠陥のプロパティを格納するキーワード引数辞書。
            :type args: dict
        """
        self.labels       = safe_getelement(args, "labels")
        self.plabels      = safe_getelement(args, "plabels")
        self.names        = safe_getelement(args, "names")

        self.atom_label    = safe_getelement(args, "atom_label")
        self.site_label    = safe_getelement(args, "site_label")
        self.charge_label  = safe_getelement(args, "charge_label")
        self.entropy_label = safe_getelement(args, "entropy_label")
        self.N0_label      = safe_getelement(args, "N0_label")
        self.Ndoped_label  = safe_getelement(args, "Ndoped_label")
        self.dH0_labels    = safe_getelement(args, "dH0_labels")

        self.atom    = safe_getelement(args, "atom")
        self._site    = safe_getelement(args, "site")
        self.name    = safe_getelement(args, "name")
        self.charge  = safe_getelement(args, "charge")
        self.entropy = safe_getelement(args, "entropy")
        self.N0      = safe_getelement(args, "N0")
        self.Ndoped  = safe_getelement(args, "Ndoped")
        self.dH0s    = safe_getelement(args, "dH0s")

    @property
    def site(self):
        """
        概要:
            欠陥サイトのプロパティを取得します。
        戻り値:
            :returns: 欠陥サイトの名前。
            :rtype: str
        """
        return self._site

class tkDefects():
    """
    概要:
        複数の欠陥情報を管理するクラス。
    詳細説明:
        Excelファイルから欠陥形成エネルギーを読み込み、セル体積を考慮して
        欠陥密度、キャリア密度、フェルミ準位の計算を行います。
        また、欠陥の遷移準位の探索と保存も担当します。
    """
    def __init__(self, **args):
        """
        概要:
            tkDefectsクラスのコンストラクタ。
        詳細説明:
            欠陥に関する各種リストやパスを初期化します。
        引数:
            :param args: 欠陥のプロパティを格納するキーワード引数辞書。
            :type args: dict
        """
        self.path         = None
        self.file_version = None
        self.V            = None    # in A^3
        self.labels       = safe_getelement(args, "labels")
        self.plabels      = safe_getelement(args, "plabels")
        self.atoms        = safe_getelement(args, "atoms")
        self.sites        = safe_getelement(args, "sites")
        self.names        = safe_getelement(args, "names")
        self.defects      = []

    def get_names(self) -> list[str]:
        """
        概要:
            登録されているユニークな欠陥名のリストを取得します。
        戻り値:
            :returns: ユニークな欠陥名のリスト。
            :rtype: list
        """
        names = {}
        for i in range(self.ndefects):
            d = self.defects[i]
            names[d.name] = 1
        return list(names.keys())

    def cal_dE(self, iPoint: int, idefect: int, T: float, EF: float) -> float:
        """
        概要:
            特定の欠陥の形成エネルギー(dE)を計算します。
        詳細説明:
            指定されたポイント、欠陥ID、温度、フェルミ準位に基づいて、
            欠陥の形成エンタルピーと電荷によるエネルギー項を合計します。
        引数:
            :param iPoint: 計算に使用する形成エンタルピーのインデックス。
            :type iPoint: int
            :param idefect: 欠陥リスト内の欠陥のインデックス。
            :type idefect: int
            :param T: 温度 (K)。
            :type T: float
            :param EF: フェルミ準位 (eV)。
            :type EF: float
        戻り値:
            :returns: 計算された欠陥形成エネルギー (eV)。
            :rtype: float
        """
        d = self.defects[idefect]
        dE = d.dH0s[iPoint] + d.charge * EF
#        print("dE=", d.entropy, d.entropy * kB / e * T, d.dH0s[iPoint], dE)
        return dE

    '''
    def cal_dG(self, iPoint, idefect, T, EF):
        d = self.defects[idefect]
        dG = d.dH0s[iPoint] + d.charge * EF - d.entropy * kB / e * T
#        print("dG=", d.entropy, d.entropy * kB / e * T, d.dH0s[iPoint], dG)
        return dG
    '''

    def calculate_total_defect_densities(self, Tdef: float, EF: float, ECmin: float, ECmax: float, EVmin: float, EVmax: float) -> tuple[dict, list, dict, list, float]:
        """
        概要:
            指定された欠陥凍結温度とフェルミ準位で、全欠陥の密度を計算します。
        詳細説明:
            各欠陥サイトの分配関数を考慮し、個々の欠陥の電荷状態に応じた密度と、
            全欠陥による総電荷を算出します。
            戻り値は各サイトの分配関数辞書、各欠陥状態の密度リスト、Tdefにおける総欠陥密度の辞書、
            各欠陥状態の形成エネルギーリスト、全欠陥による総電荷のタプルです。
        引数:
            :param Tdef: 欠陥が凍結している温度 (K)。
            :type Tdef: float
            :param EF: フェルミ準位 (eV)。
            :type EF: float
            :param ECmin: 伝導帯のエネルギー下限 (eV)。
            :type ECmin: float
            :param ECmax: 伝導帯のエネルギー上限 (eV)。
            :type ECmax: float
            :param EVmin: 価電子帯のエネルギー下限 (eV)。
            :type EVmin: float
            :param EVmax: 価電子帯のエネルギー上限 (eV)。
            :type EVmax: float
        戻り値:
            :returns: 分配関数辞書、密度リスト、総密度辞書、形成エネルギーリスト、総電荷のタプル。
            :rtype: tuple
        """
#        print("*calculate_total_defect_densities at Tdef = {} K:".format(Tdef))
        kBTedef = kB * Tdef / e

# Calculate distribution function
        ZS = {}
        for id in range(self.ndefects):
#            print("calculate_defect_densities: id=", id, self.defects[id].site)
            d = self.defects[id]
            ZS[d.site] = 1.0
#        exit()            

        for id in range(self.ndefects):
            d = self.defects[id]
#            dH = d.dH0s[iPoint] + d.charge * EF
            dE = self.cal_dE(iPoint, id, Tdef, EF)
            Kdef = dE / kBTedef
#            dG = self.cal_dG(iPoint, id, Tdef, EF)
#            Kdef = dG / kBTedef
            if Kdef > nexp:
                Kdef = nexp
            elif Kdef <= -nexp:
                Kdef = -nexp
            ZS[d.site] += exp(-Kdef)

# Calculate defect densities
        Nds = []
        dEs = []
#        dGs = []
        Qtot = 0.0
        for id in range(self.ndefects):
            d = self.defects[id]
#            dH = d.dH0s[iPoint] + d.charge * EF
            dE = self.cal_dE(iPoint, id, Tdef, EF)
            Ke = dE / kBTedef
#            dG = self.cal_dG(iPoint, id, Tdef, EF)
#            Ke = dG / kBTedef
            if Ke > nexp:
                Ke = nexp
            elif Ke < -nexp:
                Ke = -nexp
            n = d.N0 / (self.V * 1.0e-24) * exp(-Ke) / ZS[d.site]
            Nds.append(n)
            dEs.append(dE)
#            dGs.append(dG)
            Qtot += d.charge * n
#            print("T, K=", id, d.atom, d.site, d.charge, dH, Tdef, Kdef)
#            print("   n=", dH, n, d.charge, Qtot)

        NdsTdef = {}
        for id in range(self.ndefects):
            d = self.defects[id]
            label  = "{}_{}".format(d.atom, d.site)
            labelq = "{}^{}".format(label, d.charge)
            if labelq in NdsTdef.keys(): # Original: if label in NdsTdef.keys():
#                NdsTdef[label] += Nds[id]
                NdsTdef[labelq] += Nds[id]
            else:
#                NdsTdef[label] = Nds[id]
                NdsTdef[labelq] = Nds[id]

        return ZS, Nds, NdsTdef, dEs, Qtot
#        return ZS, Nds, NdsTdef, dGs, Qtot

    def calculate_charged_defect_densities(self, Te: float, NdsTdef: dict, EF: float, ECmin: float, ECmax: float, EVmin: float, EVmax: float) -> tuple[dict, list, dict, list, float]:
        """
        概要:
            電子温度と凍結された欠陥密度を用いて、電荷を帯びた欠陥密度を計算します。
        詳細説明:
            まず、各サイトの全欠陥密度を合計し、電子温度での各電荷状態の分配関数を計算します。
            その後、各電荷状態の欠陥密度と、全欠陥による総電荷を算出します。
            戻り値は各サイトの分配関数辞書、各欠陥状態の密度リスト、Tdefにおける総欠陥密度辞書、
            各欠陥状態の形成エネルギーリスト、全欠陥による総電荷のタプルです。
        引数:
            :param Te: 電子温度 (K)。
            :type Te: float
            :param NdsTdef: 欠陥凍結温度Tdefで計算された総欠陥密度の辞書 (cm^-3)。
            :type NdsTdef: dict
            :param EF: フェルミ準位 (eV)。
            :type EF: float
            :param ECmin: 伝導帯のエネルギー下限 (eV)。
            :type ECmin: float
            :param ECmax: 伝導帯のエネルギー上限 (eV)。
            :type ECmax: float
            :param EVmin: 価電子帯のエネルギー下限 (eV)。
            :type EVmin: float
            :param EVmax: 価電子帯のエネルギー上限 (eV)。
            :type EVmax: float
        戻り値:
            :returns: 分配関数辞書、密度リスト、総密度辞書、形成エネルギーリスト、総電荷のタプル。
            :rtype: tuple
        """
#        print("*calculate_charged_defect_densities at Te = {} K:".format(Te))
        kBTe = kB * Te / e

# calculate total defect densities including all charged states
        NdsTdef_tot = {}
        for id in range(self.ndefects):
            d = self.defects[id]
            label = f"{d.atom}_{d.site}"
            labelq = f"{d.atom}_{d.site}^{d.charge}"

            if label not in NdsTdef_tot.keys():
                NdsTdef_tot[label] = NdsTdef[labelq]
            else:
                NdsTdef_tot[label] += NdsTdef[labelq]

        ZSA = {}
        for label, N in NdsTdef_tot.items():
            ZSA[label] = 0.0
        '''
        for id in range(self.ndefects):
            d = self.defects[id]
            label = "{}_{}".format(d.atom, d.site)
            ZSA[label] = 0.0
        '''

        for id in range(self.ndefects):
            d = self.defects[id]
            label = "{}_{}".format(d.atom, d.site)

#            dH = d.dH0s[iPoint] + d.charge * EF
            dE = self.cal_dE(iPoint, id, Te, EF)
            Ke = dE / kBTe
#            dG = self.cal_dG(iPoint, id, Te, EF)
#            Ke = dG / kBTe
            if Ke > nexp:
                Ke = nexp
            elif Ke <= -nexp:
                Ke = -nexp

            ZSA[label] += exp(-Ke)

# Calculate defect densities
        Nds = make_matrix1(self.ndefects)
        dEs = make_matrix1(self.ndefects)
#        dGs = make_matrix1(self.ndefects)
        Qtot = 0.0
        for id in range(self.ndefects):
            d = self.defects[id]
            label = f"{d.atom}_{d.site}"
            labelq = f"{d.atom}_{d.site}^{d.charge}"

#            dH = d.dH0s[iPoint] + d.charge * EF
            dE = self.cal_dE(iPoint, id, Te, EF)
            Ke = dE / kBTe
#            dG = self.cal_dG(iPoint, id, Te, EF)
#            Ke = dG / kBTe
#           print("k=", kBTe, dH, d.N0, d.N0 * exp(-dH / kBTedef), self.V)
            if Ke > nexp:
               Ke = nexp
            elif Ke <= -nexp:
                Ke = -nexp

#            n = NdsTdef[labelq] * exp(-Ke) / ZSA[label]  # cm^-3
#            n = NdsTdef[label] * exp(-Ke) / ZSA[label]  # cm^-3
            n = NdsTdef_tot[label] * exp(-Ke) / ZSA[label]  # cm^-3
            Nds[id] = n
            dEs[id] = dE
#            dGs[id] = dG
            Qtot += d.charge * n

#        print("NdsTdef=", NdsTdef)
#        print("NdsTdef_tot=", NdsTdef_tot)
#        print("ZSA=", ZSA)
#        print("Nds=", Nds)
#        exit()
        
        return ZSA, Nds, NdsTdef, dEs, Qtot
#        return ZSA, Nds, NdsTdef, dGs, Qtot

    def calculate_defect_densities(self, T: float, Tdef: float, NdsTdef: dict | None, EF: float, ECmin: float, ECmax: float, EVmin: float, EVmax: float) -> tuple[dict, list, dict, list, float]:
        """
        概要:
            欠陥の密度を計算します。
        詳細説明:
            NdsTdefがNoneの場合は、凍結されていない欠陥の総密度を計算します。
            NdsTdefが与えられている場合は、凍結された総欠陥密度を基に、電荷を帯びた欠陥密度を計算します。
            戻り値は分配関数辞書、密度リスト、総密度辞書、形成エネルギーリスト、総電荷のタプルです。
        引数:
            :param T: 電子温度 (K)。
            :type T: float
            :param Tdef: 欠陥が凍結している温度 (K)。
            :type Tdef: float
            :param NdsTdef: 欠陥凍結温度Tdefで計算された総欠陥密度の辞書。Noneの場合、未凍結として計算。
            :type NdsTdef: dict | None
            :param EF: フェルミ準位 (eV)。
            :type EF: float
            :param ECmin: 伝導帯のエネルギー下限 (eV)。
            :type ECmin: float
            :param ECmax: 伝導帯のエネルギー上限 (eV)。
            :type ECmax: float
            :param EVmin: 価電子帯のエネルギー下限 (eV)。
            :type EVmin: float
            :param EVmax: 価電子帯のエネルギー上限 (eV)。
            :type EVmax: float
        戻り値:
            :returns: 分配関数辞書、密度リスト、総密度辞書、形成エネルギーリスト、総電荷のタプル。
            :rtype: tuple
        """
        if NdsTdef is None:
            return self.calculate_total_defect_densities(Tdef, EF, ECmin, ECmax, EVmin, EVmax)
        else:
            return self.calculate_charged_defect_densities(T, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax)

    def calculate_densities(self, T: float, Tdef: float, NdsTdef: dict | None, EF: float, ECmin: float, ECmax: float, EVmin: float, EVmax: float) -> tuple[dict, float, float, list, dict, list, float]:
        """
        概要:
            電子、正孔、および欠陥のキャリア密度と総電荷を計算します。
        詳細説明:
            電子温度が欠陥凍結温度以下の場合、欠陥密度はTdefで凍結されます。
            TがTdefより高い場合、欠陥密度はTで再計算されます。
            最終的に各サイトの分配関数、電子密度、正孔密度、各欠陥密度リスト、総密度辞書、
            各欠陥の形成エネルギーリスト、および総電荷のタプルを返します。
        引数:
            :param T: 電子温度 (K)。
            :type T: float
            :param Tdef: 欠陥が凍結している温度 (K)。
            :type Tdef: float
            :param NdsTdef: 欠陥凍結温度Tdefで計算された総欠陥密度の辞書。Noneの場合、未凍結として計算。
            :type NdsTdef: dict | None
            :param EF: フェルミ準位 (eV)。
            :type EF: float
            :param ECmin: 伝導帯のエネルギー下限 (eV)。
            :type ECmin: float
            :param ECmax: 伝導帯のエネルギー上限 (eV)。
            :type ECmax: float
            :param EVmin: 価電子帯のエネルギー下限 (eV)。
            :type EVmin: float
            :param EVmax: 価電子帯のエネルギー上限 (eV)。
            :type EVmax: float
        戻り値:
            :returns: 分配関数辞書、電子密度、正孔密度、欠陥密度リスト、総密度辞書、形成エネルギーリスト、総電荷のタプル。
            :rtype: tuple
        """
# If T <= Tdef, defect densities are frozen at Tdef
# If T > Tdef, defect densities are recalculated at T
        if T <= Tdef:
            Z, Nds, NdsTdef, dGs, Qtot = self.calculate_defect_densities(T, Tdef, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax)
        else:
            Z, Nds, NdsTdef, dGs, Qtot = self.calculate_defect_densities(T, T, None, EF, ECmin, ECmax, EVmin, EVmax)

        ne = Ne(T, EF, ECmin, ECmax)
        nh = Nh(T, EF, EVmin, EVmax)
        Qtot += nh - ne

        return Z, ne, nh, Nds, NdsTdef, dGs, Qtot

    def dQ(self, T: float, Tdef: float, NdsTdef: dict | None, EF: float, ECmin: float, ECmax: float, EVmin: float, EVmax: float) -> float:
        """
        概要:
            全電荷の中性条件からの偏差(dQ)を計算します。
        詳細説明:
            calculate_densitiesメソッドを呼び出して、電子、正孔、欠陥の密度から総電荷を計算し、
            その値を返します。この関数はフェルミ準位探索のための目的関数として使用されます。
        引数:
            :param T: 電子温度 (K)。
            :type T: float
            :param Tdef: 欠陥が凍結している温度 (K)。
            :type Tdef: float
            :param NdsTdef: 欠陥凍結温度Tdefで計算された総欠陥密度の辞書。
            :type NdsTdef: dict | None
            :param EF: フェルミ準位 (eV)。
            :type EF: float
            :param ECmin: 伝導帯のエネルギー下限 (eV)。
            :type ECmin: float
            :param ECmax: 伝導帯のエネルギー上限 (eV)。
            :type ECmax: float
            :param EVmin: 価電子帯のエネルギー下限 (eV)。
            :type EVmin: float
            :param EVmax: 価電子帯のエネルギー上限 (eV)。
            :type EVmax: float
        戻り値:
            :returns: 全電荷の偏差(単位電荷)。
            :rtype: float
        """
        Z, ne, nh, Nds, NdsTdef, dGs, Qtot = self.calculate_densities(T, Tdef, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax)
        return Qtot

    def diff(self, h: float, T: float, Tdef: float, NdsTdef: dict | None, EF: float, ECmin: float, ECmax: float, EVmin: float, EVmax: float) -> float:
        """
        概要:
            全電荷の中性条件からの偏差(dQ)のフェルミ準位による数値微分を計算します。
        詳細説明:
            中心差分法を用いてdQ/dEFを計算します。これはニュートン法のヤコビアンとして使用されます。
        引数:
            :param h: 微小なEFの摂動量 (eV)。
            :type h: float
            :param T: 電子温度 (K)。
            :type T: float
            :param Tdef: 欠陥が凍結している温度 (K)。
            :type Tdef: float
            :param NdsTdef: 欠陥凍結温度Tdefで計算された総欠陥密度の辞書。
            :type NdsTdef: dict | None
            :param EF: フェルミ準位 (eV)。
            :type EF: float
            :param ECmin: 伝導帯のエネルギー下限 (eV)。
            :type ECmin: float
            :param ECmax: 伝導帯のエネルギー上限 (eV)。
            :type ECmax: float
            :param EVmin: 価電子帯のエネルギー下限 (eV)。
            :type EVmin: float
            :param EVmax: 価電子帯のエネルギー上限 (eV)。
            :type EVmax: float
        戻り値:
            :returns: dQのEFによる数値微分。
            :rtype: float
        """
        dQp = self.dQ(T, Tdef, NdsTdef, EF + h, ECmin, ECmax, EVmin, EVmax)
        dQm = self.dQ(T, Tdef, NdsTdef, EF - h, ECmin, ECmax, EVmin, EVmax)
        return (dQp - dQm) / 2.0 / h

    def find_EF(self, T: float, Tdef: float, NdsTdef: dict | None, ECmin: float, ECmax: float, EVmin: float, EVmax: float,
            callbackfunc: callable = None, 
            eps_bisec: float = 1.0e-1, nmaxiter_bisection: int = 300, initial: list[float] = [-1.0, 4.0],
            eps_newton: float = 1.0e-6, nmaxiter_newton: int = 300, dump_newton: float = 0.3, h_newton: float = 1.0e-6,
            is_print: bool = False) -> tuple[float | None, float | None]:
        """
        概要:
            全電荷中性条件を満たすフェルミ準位を探索します。
        詳細説明:
            二分法とニュートン法を組み合わせて、全電荷がゼロとなるフェルミ準位を探索します。
            最初に広い範囲で二分法を実行し、その後、より精密なニュートン法を適用します。
            戻り値は収束したフェルミ準位と、そのときの全電荷の値のタプルです。
            収束しない場合はNoneのタプルを返します。
        引数:
            :param T: 電子温度 (K)。
            :type T: float
            :param Tdef: 欠陥が凍結している温度 (K)。
            :type Tdef: float
            :param NdsTdef: 欠陥凍結温度Tdefで計算された総欠陥密度の辞書。
            :type NdsTdef: dict | None
            :param ECmin: 伝導帯のエネルギー下限 (eV)。
            :type ECmin: float
            :param ECmax: 伝導帯のエネルギー上限 (eV)。
            :type ECmax: float
            :param EVmin: 価電子帯のエネルギー下限 (eV)。
            :type EVmin: float
            :param EVmax: 価電子帯のエネルギー上限 (eV)。
            :type EVmax: float
            :param callbackfunc: 探索の各ステップで呼び出されるコールバック関数。
            :type callbackfunc: callable | None
            :param eps_bisec: 二分法の収束判定閾値。
            :type eps_bisec: float
            :param nmaxiter_bisection: 二分法の最大反復回数。
            :type nmaxiter_bisection: int
            :param initial: 二分法の初期探索範囲のリスト。
            :type initial: list
            :param eps_newton: ニュートン法の収束判定閾値。
            :type eps_newton: float
            :param nmaxiter_newton: ニュートン法の最大反復回数。
            :type nmaxiter_newton: int
            :param dump_newton: ニュートン法のダンプ係数。
            :type dump_newton: float
            :param h_newton: ニュートン法の数値微分のための微小量。
            :type h_newton: float
            :param is_print: 途中経過をコンソールに出力するかどうか。
            :type is_print: bool
        戻り値:
            :returns: 収束したフェルミ準位と全電荷の値のタプル。
            :rtype: tuple
        """
        EF = (initial[0] + initial[1]) / 2.0

        if nmaxiter_bisection > 0:
            solver = tkequation.Equation(
                debug_explicit = 1,
                method    = 'bisection',
                func      = lambda EF: self.dQ(T, Tdef, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax),
                xa        = initial[0],
                xb        = initial[1],
                nmaxiter  = nmaxiter_bisection,
                eps       = eps_bisec,
                callback  = callbackfunc,
                isprint   = 0
                )
            x = solver.solve()
            if solver.iter >= 0:
                if is_print:
                    print("  Convergence reached in bisection proc at iter = {}, x = {:10.6g}, f = {:8.4g}"
                            .format(solver.iter, solver.x, solver.f))
            else:
                print("")
                print("*** Error: Convergence is not reached")
                return None, None

            EF  = solver.x
            dQh = solver.f

        if nmaxiter_newton > 0:
            solver = tkequation.Equation(
                debug_explicit = 1,
                method    = 'newton', 
                func      = lambda EF: self.dQ(T, Tdef, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax),
                diff1func = lambda EF: self.diff(h_newton, T0, Tdef, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax), 
                x0        = EF,
                dump      = dump_newton,
                nmaxiter  = nmaxiter_newton,
                eps       = eps_newton,
                callback  = callbackfunc,
                isprint   = 0
                )
            x = solver.solve()
            if solver.iter >= 0:
                if is_print:
                    print("  Convergence reached in bisection proc at iter = {}, x = {:10.6g}, f = {:8.4g}"
                            .format(solver.iter, solver.x, solver.f))
            else:
                print("")
                print("*** Error: Convergence is not reached")
                return None, None

            EF  = solver.x
            dQh = solver.f

        EF  = solver.x
        dQh = solver.f
        
        return EF, dQh

    def find_all_transition_EF(self, idx0: int, groupdata: list[dict]) -> list[float]:
        """
        概要:
            特定の欠陥状態から他のすべての欠陥状態への遷移準位を計算します。
        詳細説明:
            欠陥形成エネルギーの交点から、異なる電荷状態間の遷移準位を特定します。
            電荷が同じ場合は非常に大きな値を返します。
        引数:
            :param idx0: 基準となる欠陥状態のインデックス。
            :type idx0: int
            :param groupdata: 特定の欠陥名に属する全欠陥状態のデータリスト。
            :type groupdata: list
        戻り値:
            :returns: 各欠陥状態への遷移準位のリスト。
            :rtype: list
        """
        print(f"  find_all_transition_EF: idx0={idx0}")
        ngroup = len(groupdata)
        d0 = groupdata[idx0]
        q0 = d0["charge"]
        H0 = d0["dH0"]

        EF_list = []
        for idx1 in range(ngroup):
            d1 = groupdata[idx1]
            q1 = d1["charge"]

            if idx1 == idx0 or q1 == q0:
                EF_list.append(1.0e100)
                continue

            H1 = d1["dH0"]
            EF1 =  (H0 - H1) / (q1 - q0)
            EF_list.append(EF1)
        
        return EF_list
          
    def find_next_transition_EF(self, idx0: int, prevEF: float, EFmin: float, EFmax: float, groupdata: list[dict]) -> tuple[int | None, float]:
        """
        概要:
            現在のフェルミ準位より大きく、指定範囲内の次の遷移準位を探索します。
        詳細説明:
            find_all_transition_EFを使用してすべての遷移準位を計算し、
            現在のフェルミ準位よりも大きく、かつ探索範囲内で最小の遷移準位を特定します。
            戻り値は次の遷移準位のインデックスとフェルミ準位値のタプルです。
        引数:
            :param idx0: 基準となる欠陥状態のインデックス。
            :type idx0: int
            :param prevEF: 現在のフェルミ準位 (eV)。
            :type prevEF: float
            :param EFmin: フェルミ準位の最小探索値 (eV)。
            :type EFmin: float
            :param EFmax: フェルミ準位の最大探索値 (eV)。
            :type EFmax: float
            :param groupdata: 特定の欠陥名に属する全欠陥状態のデータリスト。
            :type groupdata: list
        戻り値:
            :returns: 次の遷移準位のインデックスとフェルミ準位値のタプル。
            :rtype: tuple
        """
        print(f"  find_next_transition_EF(): idx0={idx0}  prevEF={prevEF:.4f}  EFmin={EFmin:.4f}  EFmax={EFmax:.4f}")
        EF_list = self.find_all_transition_EF(idx0, groupdata)
        print(f"  prevEF={prevEF} EF_list=", EF_list)

        nextidx = None
        nextEF = EFmax
        for i, EF in enumerate(EF_list):
            if EF < EFmin or EFmax < EF or EF <= prevEF: continue

            if EF < nextEF:
                nextidx = i
                nextEF = EF

#        print(f"     next: idx={nextidx} at EF={nextEF}")
        return nextidx, nextEF

    def find_order(self, groupdata: list[dict], idx: int, EF: float, dH: float) -> tuple[int, int, float]:
        """
        概要:
            特定のフェルミ準位における、ある欠陥状態の形成エンタルピーの順位を計算します。
        詳細説明:
            同じ欠陥名に属するすべての電荷状態について形成エンタルピーを計算し、
            入力されたペアがその中で何番目に低い形成エンタルピーを持つかを決定します。
            戻り値は順位、最小の形成エンタルピーを持つインデックス、およびその値のタプルです。
        引数:
            :param groupdata: 特定の欠陥名に属する全欠陥状態のデータリスト。
            :type groupdata: list
            :param idx: 現在の欠陥状態のインデックス。
            :type idx: int
            :param EF: フェルミ準位 (eV)。
            :type EF: float
            :param dH: 現在の欠陥状態の形成エンタルピー (eV)。
            :type dH: float
        戻り値:
            :returns: 順位、最小エンタルピーのインデックス、最小エンタルピー値のタプル。
            :rtype: tuple
        """
        print(f"  find_oder() for idx={idx}, EF={EF}, dH={dH}:")
        ndata = len(groupdata)
        dH2 = []
        q_list = []
        for id2 in range(ndata):
            d2  = groupdata[id2]
            q2  = d2["charge"]
            H2  = d2["dH0"]
            dH2.append(H2 + q2 * EF)
            q_list.append(q2)
#            print(f"find_order: id2={id2}: {d2['name']} q={d2['charge']} dH={d2['dH0']}")
#            print("  id2={:2}: EF1={:8.3g} dH1={:8.3g} dH2={:8.3g}".format(id2, EF1, dH1, dH2[id2]))

        print("    dH2=", [float(f"{v:8.6g}") for v in dH2], " for q=", q_list)
        iorder = 0
        for id2 in range(ndata):
#            print("order:", iorder, dH2[id2], dH)
            if dH2[id2] < dH - 1.0e-8:
                print(f"  current: idx={idx} ({EF:.6f}, {dH:.6g}): for id2={id2}: dH2[id2]={dH2[id2]:.6g} < dH={dH:.6g}")
                iorder += 1

        idxmin = idx
        dHmin = dH
        for id2 in range(ndata):
            if dH2[id2] < dHmin:
                dHmin  = dH2[id2]
                idxmin = id2

        print(f"     *found dHmin={dHmin:8.6g} at EF={EF}: for idxmin={idxmin:2}, iorder={iorder:2}")

# iorder: (EF, dH)が同じグループで何番目か
# idxmin, dHmin: 同じグループでEFにおいて最小のdHをもつidx,dH
        return iorder, idxmin, dHmin


    def find_all_transitions(self, name: str, groupdata: list[dict], EFmin: float, EFmax: float) -> tuple[list[list], list[list]]:
        """
        概要:
            特定の欠陥名に属するすべての遷移点と最小形成エンタルピー曲線上の点を探索します。
        詳細説明:
            フェルミ準位をEFminからEFmaxまで変化させながら、各電荷状態の形成エンタルピーを計算し、
            形成エンタルピーが最小となる電荷状態が変化する遷移準位を特定します。
            戻り値はすべての遷移点のリストと、最小形成エンタルピー曲線上の主要な点のリストのタプルです。
        引数:
            :param name: 欠陥名。
            :type name: str
            :param groupdata: 特定の欠陥名に属する全欠陥状態のデータリスト。
            :type groupdata: list
            :param EFmin: フェルミ準位の最小探索値 (eV)。
            :type EFmin: float
            :param EFmax: フェルミ準位の最大探索値 (eV)。
            :type EFmax: float
        戻り値:
            :returns: すべての遷移点のリストと最小形成エンタルピー曲線上の点のリストのタプル。
            :rtype: tuple
        """
        print("  find_all_transitions():")
        eps = 1.0e-6
        ngroupdata = len(groupdata)
        allpts = []
        minpts = []

        idxmin = 0
        d0 = groupdata[idxmin]
        q0 = d0["charge"]
        H0 = d0["dH0"]

        print()
        print(f"defect name: {d0['name']}")

        iorder, idxmin, dHmin = self.find_order(groupdata, idxmin, EF = EFmin, dH = H0 + q0 * EFmin)
        allpts.append([EFmin, dHmin, q0, groupdata[idxmin]["charge"], 0, idxmin])
        minpts.append([EFmin, dHmin, q0, groupdata[idxmin]["charge"], 0, idxmin])
        print(f"  ** first point #1: (EF, dH)=({EFmin:.6f}, {dHmin:.6g}) q={q0} idx=0")
        print()

        prevEF = EFmin
        prevq = None
        count = 2
        while 1:
            nextidx, nextEF = self.find_next_transition_EF(idxmin, prevEF, EFmin, EFmax, groupdata)
            if nextidx is None: break
            print(f"  nextidx={nextidx}  EF={nextEF:.6g}")

            d0 = groupdata[idxmin]
            q0 = d0["charge"]
            H0 = d0["dH0"]
            dH1 = H0 + q0 * nextEF
            d1 = groupdata[nextidx]
            q1 = d1['charge']
            iorder, idxmin, dHmin = self.find_order(groupdata, idxmin, nextEF, dH1)
            print(f"  ** next point #{count}: iorder={iorder} (EF, dH)=({nextEF:.6f}, {dHmin:.6g}) idx={nextidx} next_q={q1}")
            print()

            allpts.append([nextEF, dH1, q0, q1, iorder, idxmin])
#            if prevq is None or (prevq > q0):
            if iorder == 0:
               minpts.append([nextEF, dH1, q0, q1, iorder, idxmin])

            prevEF = nextEF
            prevq = q0
            count += 1

        # 最後のEFmax点での最小dHを計算し追加
        # ここでは、最後に決定されたidxmin (EFmaxにおける最小dHの電荷状態) を使用
        if ngroupdata > 0:
            last_idx = idxmin # 最後の遷移点での最小dHを持つidx
            last_d = groupdata[last_idx]
            last_q = last_d["charge"]
            last_H = last_d["dH0"]
            dH_at_EFmax = last_H + last_q * EFmax
            iorder_at_EFmax, idxmin_at_EFmax, dHmin_at_EFmax = self.find_order(groupdata, last_idx, EFmax, dH_at_EFmax)
        else:
            dH_at_EFmax = 1.0e100 # 欠陥データがない場合は大きな値を設定
            iorder_at_EFmax, idxmin_at_EFmax, dHmin_at_EFmax = 0, 0, 1.0e100 # デフォルト値を設定
            
        allpts.append([EFmax, dHmin_at_EFmax, last_q, groupdata[idxmin_at_EFmax]["charge"], iorder_at_EFmax, idxmin_at_EFmax])
        minpts.append([EFmax, dHmin_at_EFmax, last_q, groupdata[idxmin_at_EFmax]["charge"], iorder_at_EFmax, idxmin_at_EFmax])

#        if name == 'O_Ge1':
#        if name == 'V_O1':
#            print("allpts=", allpts)
#            print("minpts=", minpts)
#            exit()

        return allpts, minpts

    def get_transitionlevel_lists(self) -> tuple[list[str], list[list[list]], list[list[list]]]:
        """
        概要:
            全てのユニークな欠陥名に対して、遷移準位のリストを生成します。
        詳細説明:
            各欠陥名について遷移準位を計算し、すべての形成エンタルピー曲線上の点と、
            最小曲線上の点を取得します。戻り値は欠陥名のリスト、全遷移点のリストのリスト、
            および最小曲線上の点のリストのリストのタプルです。
        戻り値:
            :returns: 欠陥名リスト、全遷移点リスト、最小曲線上の点リストのタプル。
            :rtype: tuple
        """
        names = self.get_names()
        allpts_list = []
        minpts_list = []
        for iname in range(len(names)):
            name = names[iname]
            groupdata = self.get_groupdata(name, iPoint)
            allpts, minpts = self.find_all_transitions(name, groupdata, EFmin, EFmax)
            allpts_list.append(allpts)
            minpts_list.append(minpts)
#        exit()

        return names, allpts_list, minpts_list
    
    def save_allpts(self, path: str) -> int | None:
        """
        概要:
            全ての欠陥状態の形成エンタルピーとフェルミ準位の関係をExcelファイルに保存します。
        詳細説明:
            各欠陥の電荷状態について、フェルミ準位の範囲における形成エンタルピーの線形関係を計算し、
            Excelシートに書き出します。
        引数:
            :param path: 出力Excelファイルのパス。
            :type path: str
        戻り値:
            :returns: ファイル書き込みが成功すれば1、失敗すればNone。
            :rtype: int | None
        """
        allwb = tkExcel(path, 'w')
        if not allwb:
            return None

        for id in range(self.ndefects):
            d = self.defects[id]
            name = d.name
            q    = d.charge
            allwb.Print([name, q])
            allwb.Print(["EF(eV)", "dH(eV)"])
            allwb.Print([EFmin, d.dH0s[iPoint] + q * EFmin])
            allwb.Print([EFmax, d.dH0s[iPoint] + q * EFmax])

        allwb.Close()
        return 1

    def save_minpts(self, path: str) -> int | None:
        """
        概要:
            各欠陥タイプにおける最小形成エンタルピー曲線上の点をExcelファイルに保存します。
        詳細説明:
            フェルミ準位に対して形成エンタルピーが最小となる曲線上の遷移準位点とその情報を
            Excelシートに書き出します。
        引数:
            :param path: 出力Excelファイルのパス。
            :type path: str
        戻り値:
            :returns: ファイル書き込みが成功すれば1、失敗すればNone。
            :rtype: int | None
        """
        minwb = tkExcel(path, 'w')
        if not minwb:
            return None

        names = self.get_names()
        for iname in range(len(names)):
            name = names[iname]
            groupdata = self.get_groupdata(name, iPoint)
            allpts, minpts = self.find_all_transitions(name, groupdata, EFmin, EFmax)

            minwb.Print([name])
            minwb.Print(["EF(eV)", "dH(eV)", "q(-)", "q(+)"])

            for ipt in range(len(minpts)):
                minwb.Print([minpts[ipt][0], minpts[ipt][1], minpts[ipt][2], minpts[ipt][3]])

        minwb.Close()
        return 1

    def get_groupdata(self, name: str, iPoint: int) -> list[dict]:
        """
        概要:
            指定された欠陥名に属する全電荷状態のデータを取得します。
        詳細説明:
            指定された名前と一致する欠陥を抽出し、それらの情報を辞書形式でリストにまとめます。
            リストは電荷の降順にソートされます。
        引数:
            :param name: 抽出する欠陥名。
            :type name: str
            :param iPoint: 使用する形成エンタルピーのインデックス。
            :type iPoint: int
        戻り値:
            :returns: 抽出された欠陥データのリスト。
            :rtype: list
        """
        groupdata = []
        for i in range(self.ndefects):
            d = self.defects[i]
            if name != d.name: continue
            
            p = {
                'defects': d,
                'atom'   : d.atom,
                'site'   : d.site,
                'name'   : d.name,
                'charge' : d.charge,
                'dH0'    : d.dH0s[iPoint]
                }
            groupdata.append(p)
        groupdata.sort(key = lambda x: -x["charge"])
        return groupdata

    def SetVolume(self, V: float) -> float:
        """
        概要:
            ユニットセルの体積を設定します。
        引数:
            :param V: ユニットセルの体積 (A^3)。
            :type V: float
        戻り値:
            :returns: 設定されたユニットセルの体積 (A^3)。
            :rtype: float
        """
        self.V = V
        return self.V
        
    def Volume(self) -> float | None:
        """
        概要:
            ユニットセルの体積を取得します。
        戻り値:
            :returns: 設定されたユニットセルの体積 (A^3)。
            :rtype: float | None
        """
        return self.V
    
    def add(self, defect: 'tkDefect'):
        """
        概要:
            欠陥オブジェクトをリストに追加します。
        引数:
            :param defect: 追加するtkDefectオブジェクト。
            :type defect: tkDefect
        """
        self.defects.append(defect)

    def get(self, idx: int):
        """
        概要:
            指定されたインデックスの欠陥オブジェクトを取得します。
        引数:
            :param idx: 取得する欠陥のインデックス。
            :type idx: int
        戻り値:
            :returns: 指定されたインデックスのtkDefectオブジェクト。
            :rtype: tkDefect | None
        """
        try:
            return self.defects[idx]
        except:
            print("")
            print("Warning in tkDefects.get: No defect for index=", idx)
            print("")
            return None

    def Print(self, PrintLabels: bool = False, outfp: object = None):
        """
        概要:
            登録されている欠陥情報を標準出力または指定されたファイルポインタに出力します。
        詳細説明:
            欠陥の数、ポイント数、セル体積、および各欠陥の詳細情報を整形して表示します。
        引数:
            :param PrintLabels: ラベル情報を出力するかどうか。
            :type PrintLabels: bool
            :param outfp: 出力先のファイルポインタ。
            :type outfp: object | None
        """
        mprint("# of defect points: ", self.npoints, fp = outfp)
        mprint("# of defects      : ", self.ndefects, fp = outfp)
        if self.V:
            mprint("Unit cell volume: {:10.6g} A^3".format(self.V))
        if PrintLabels:
            mprint("labels : ", self.labels, fp = outfp)
            mprint("plabels: ", self.plabels, fp = outfp)
            mprint("names  : ", self.names, fp = outfp)
        mprint(" {:>6} ({:>2} {:>4})  {:>5} {:^6} {:^8} {:^12}  ".format("Defect", "at", "site", "q", "dS/kB", "N0(sites)", "Ndoped(cm^-3)"), end = '', fp = outfp)
        for j in range(self.npoints):
             mprint("   {:^10}".format("dH0(eV) @ {}".format(self.plabels[j])), end = '', fp = outfp)
        mprint("", fp = outfp)

        for i in range(self.ndefects):
            d = self.get(i)
            mprint(f" {d.name:6} ({d.atom:2} {d.site:4}): {d.charge:5} {d.entropy:6.2g} {d.N0:8.4g} {d.Ndoped:12.4g}  ", end = '', fp = outfp)
            for j in range(self.npoints):
                mprint("   {:10.6g}".format(d.dH0s[j]), end = '', fp = outfp)
            mprint("", fp = outfp)

    def read_excel(self, infile: str, V: float | None) -> 'tkDefects':
        """
        概要:
            欠陥形成エンタルピーを含むExcelファイルを読み込みます。
        詳細説明:
            指定されたExcelファイルから欠陥に関する情報を抽出し、ファイル形式のバージョンを
            自動判別して適切にデータをパースします。セル体積Vが与えられた場合はそれを設定します。
        引数:
            :param infile: 欠陥形成エンタルピーを含むExcelファイルのパス。
            :type infile: str
            :param V: ユニットセルの体積 (A^3)。
            :type V: float | None
        戻り値:
            :returns: 読み込まれたデータで構築されたtkDefectsオブジェクト自身。
            :rtype: tkDefects
        """
        self.path = infile
        
        if V:
            self.V = V

        try:
            wb = tkExcel(infile, 'r')
        except:
            app.terminate("Error to read [{}]".format(infile), pause = True)

        version = wb.Get(0, 0)
#        print("version:", version)

        if 'Version' in version:
            line0    = 0
            col0     = 1
            nlabels0 = 6
        else:
            line0    = 0
            col0     = 0
            nlabels0 = 5

        labels   = []
        idx = col0
        while 1:
            v = wb.Get(line0, idx)
            if v is None:
                break

            labels.append(v)
            idx += 1

        nlabels = len(labels)

        values = []
        for i in range(nlabels):
            values.append([])

#       print("")
        iline = line0 + 1
        is_last = False
        while 1:
            for irow in range(nlabels):
                v = wb.Get(iline, irow + col0)
#                print("iline=", iline, v)
                if v is None or v == '':
                    is_last = True
                    break

                if irow <= 1 + col0:
                    values[irow].append(v)
                else:
                    values[irow].append(pfloat(v))

            if is_last:
                break
            iline += 1

        wb.Close()

#        print("Labels:", labels)
#        print("values=", values)

        ndefects = len(values[0])
        npoints  = nlabels - nlabels0
        plabels  = labels[nlabels0:]
        labels   = labels[0:nlabels0]

        atoms     = values[0]
        sites     = values[1]
        charges   = values[2]
        if 'Version' in version:
            entropies = values[3]
            N0s       = values[4]
            Nds       = values[5]
            pvalues   = values[6:]
        else:
            entropies = make_matrix1(ndefects, 0.0)
            N0s       = values[3]
            Nds       = values[4]
            pvalues   = values[5:]

        defects = self
        defects.file_version = version
        defects.ndefects = ndefects
        defects.npoints  = npoints
        defects.labels   = labels
        defects.plabels  = plabels
        defects.atoms    = atoms
        defects.sites    = sites
        defects.names    = []
        print("atoms=", defects.atoms)
        print("sites=", defects.sites)
        for i in range(defects.ndefects):
            defects.names.append(defects.atoms[i] + '_' + defects.sites[i])
        for i in range(defects.ndefects):
            d = tkDefect()
            d.atom_label    = labels[0]
            d.site_label    = labels[1]
            d.charge_label  = labels[2]
            if 'Version' in version:
                d.entropy_label = labels[3]
                d.N0_label      = labels[4]
                d.Ndoped_label  = labels[5]
            else:
                d.entropy_label = 'S/kB'
                d.N0_label      = labels[3]
                d.Ndoped_label  = labels[4]
            d.dH0_labels    = plabels
            d.atom          = defects.atoms[i]
            d._site         = defects.sites[i]
            d.name          = defects.names[i]
            d.charge        = charges[i]
            d.entropy       = entropies[i]
            d.N0            = N0s[i]
#            if self.V:
#                d.N0     = N0s[i] / (self.V * 1.0e-24)
#            else:
#                d.N0     = N0s[i]
            d.Ndoped = Nds[i]
            d.dH0s   = []
            for j in range(npoints):
                d.dH0s.append(pvalues[j][i])
            defects.add(d)

        return defects


#=============================
# Treat argments
#=============================
def usage(app: 'tkApplication' = None):
    """
    概要:
        スクリプトのコマンドライン引数と使用方法を表示します。
    詳細説明:
        プログラムの実行に必要な引数とその例を示します。
        特にモードに応じた引数の指定方法を説明します。
    引数:
        :param app: tkApplicationのインスタンス。
        :type app: tkApplication | None
    """
    global CAR_path
    
    if CAR_path == '':
        CAR_path = '.'
    
    print("")
    print("Usage: Variables in () are optional")
    print("  Input files: INCAR, POSCAR, POTCAR, OUTCAR, DOSCAR, EIGENVAL of perfect crystal model")
    print("               Excel (.xlsx) file of defect formation energies")
    print(" python {} mode (other args)".format(sys.argv[0], ))
    print("   NOTE: If T(electron) <= T(defect) defect densities are frozen at T(defect)")
    print("         If T(electron) >  T(defect) defect densities are calculated at T(electron)")
    print(" (i) python {} EF [min|all] dH_path CAR_path CAR_path(defect) iPoint T(electron) T(defect) EF(defect) EFmin EFmax nEF".format(sys.argv[0]))
    print("     ex: python {} EF min {} {} {} {} {} {} {} {} {} {}"
                .format(sys.argv[0], dH_path, CAR_path, CAR_path_defect, iPoint, T0, Tdef, EFdef, EFmin, EFmax, nEF))
    print("     ex: python {} EF all {} {} {} {} {} {} {} {} {} {}"
                .format(sys.argv[0], dH_path, CAR_path, CAR_path_defect, iPoint, T0, Tdef, EFdef, EFmin, EFmax, nEF))
    print(" (ii) python {} T dH_path, CAR_path iPoint T(defect) Tmin Tmax nT".format(sys.argv[0]))
    print("     ex: python {} T {} {} {} {} {} {} {} {} {}"
                .format(sys.argv[0], dH_path, CAR_path, CAR_path_defect, iPoint, Tdef, EFdef, Tmin, Tmax, nT))

def updatevars():
    """
    概要:
        コマンドライン引数からグローバル変数を更新します。
    詳細説明:
        sys.argvを解析し、実行モードに応じて欠陥データパス、VASP計算パス、
        温度範囲、フェルミ準位範囲などのグローバル設定値を更新します。
    """
    global mode, plot_mode
    global dH_path, CAR_path, Vcell, CAR_path_defect
    global nDOS, Emin, Emax, Estep
    global T0, Tdef, EFdef, EF0
    global EFmin, EFmax, nEF, EFstep, dHmin, dHmax
    global Tmin, Tmax, nT, Tstep
    global iPoint
    global ignore_warning

    argv = sys.argv
    if len(argv) == 1:
        usage()
        exit()

    mode        = getarg( 1, mode)
    if mode == 'EF':
        plot_mode   = getarg( 2, plot_mode)
        dH_path     = getarg( 3, dH_path)
        CAR_path    = getarg( 4, CAR_path)
        CAR_path_defect = getarg( 5, CAR_path_defect)
        iPoint = getintarg  ( 6, iPoint)
        T0     = getfloatarg( 7, T0)
        Tdef   = getfloatarg( 8, Tdef)
        EFdef  = getarg     ( 9, EFdef)
        EFmin  = getarg(10, EFmin)
        EFmax  = getarg(11, EFmax)
        EFmin = pfloat(EFmin, None)
        EFmax = pfloat(EFmax, None)
        nEF    = getintarg  (12, nEF)
        dHmin  = getarg(13, dHmin)
        dHmax  = getarg(14, dHmax)
        dHmin = pfloat(dHmin, None)
        dHmax = pfloat(dHmax, None)
        ignore_warning = getintarg(15, ignore_warning)

#        EFstep = (EFmax - EFmin) / (nEF - 1)
    if mode == 'T':
        dH_path     = getarg     ( 2, dH_path)
        CAR_path    = getarg     ( 3, CAR_path)
        CAR_path_defect = getarg  ( 4, CAR_path_defect)
        iPoint      = getintarg  ( 5, iPoint)
        Tdef        = getfloatarg( 6, Tdef)
        EFdef       = getarg     ( 7, EFdef)
        Tmin        = getfloatarg( 8, Tmin)
        Tmax        = getfloatarg( 9, Tmax)
        nT          = getintarg  (10, nT)

        ignore_warning = getarg(11, ignore_warning,)

        Tstep  = (Tmax - Tmin) / (nT - 1)

#=============================
# other functions
#=============================
def savecsv(outfile: str, header: list[str], datalist: list[list[float]]):
    """
    概要:
        データリストをCSVファイルに保存します。
    詳細説明:
        ヘッダー行とそれに続くデータ行をCSV形式で指定されたファイルに書き込みます。
        データは行ごとに整形されて出力されます。
    引数:
        :param outfile: 出力CSVファイルのパス。
        :type outfile: str
        :param header: CSVファイルのヘッダー行のリスト。
        :type header: list
        :param datalist: 書き込むデータ。
        :type datalist: list
    """
    try: 
        print("Write to [{}]".format(outfile))
        f = open(outfile, 'w')
    except:
#    except IOError:
        print("Error: Can not write to [{}]".format(outfile))
    else:
        fout = csv.writer(f, lineterminator='\n')
        fout.writerow(header)
#        fout.writerows(data)
        for i in range(0, len(datalist[0])):
            a = []
            for j in range(len(datalist)):
                a.append(datalist[j][i])
            fout.writerow(a)
        f.close()

def read_csv(infile: str, xmin: float | None = None, xmax: float | None = None, delimiter: str = ',') -> tuple[list[str], list[float], list[float]]:
    """
    概要:
        CSVファイルを読み込み、指定された範囲のデータを抽出します。
    詳細説明:
        最初の列をx軸データ、2番目の列をy軸データとして読み込みます。
        xminとxmaxが指定されている場合、x軸データがその範囲内の行のみを抽出します。
        戻り値はヘッダー行のリスト、x軸データのリスト、y軸データのリストのタプルです。
    引数:
        :param infile: 入力CSVファイルのパス。
        :type infile: str
        :param xmin: x軸データの最小値。
        :type xmin: float | None
        :param xmax: x軸データの最大値。
        :type xmax: float | None
        :param delimiter: CSVファイルの区切り文字。
        :type delimiter: str
    戻り値:
        :returns: ヘッダーリストと抽出されたx軸データリスト、y軸データリストのタプル。
        :rtype: tuple
    """
    print("xrange=", xmin, xmax)
    data = []
    try:
        infp = open(infile, "r")
        f = csv.reader(infp, delimiter = delimiter)
        header = next(f)
        print("header=", header)
        for j in range(len(header)):
            data.append([])

        for row in f:
            x = pfloat(row[0])
            if xmin is not None and xmin <= x <= xmax:
                y = pfloat(row[1])
                data[0].append(x)
                data[1].append(y)
    except:
        print("Error: Can not read [{}]".format(infile))
        exit()
    return header, data[0], data[1]

# define the DOS function
fdos = None
def DOS(E: float) -> float:
    """
    概要:
        補間関数を使用して、指定されたエネルギーにおける状態密度を返します。
    詳細説明:
        状態密度は予めVASPのDOSCARファイルから読み込まれ、fdosとして補間関数が設定されています。
        入力エネルギーが定義範囲外の場合、エラーメッセージを出力してプログラムを終了します。
    引数:
        :param E: 状態密度を評価するエネルギー (eV)。
        :type E: float
    戻り値:
        :returns: 指定されたエネルギーにおける状態密度 (states/cm^3)。
        :rtype: float
    """
    global E_norm, dos_raw
    global Enorm_min, Enorm_max, Enorm_step, nDOS
    global fdos

    if E < Enorm_min or Enorm_max < E:
        app.terminate("Error in DOS(E): Given E={} exceeds the DOS E range [{}, {}]".format(E, Enorm_min, Enorm_max), usage = usage, pause = True)
        exit()

    return fdos(E)

# Fermi-Dirac distribution function
def fe(E: float, T: float, EF: float) -> float:
    """
    概要:
        フェルミ・ディラック分布関数を計算します(電子占有確率)。
    詳細説明:
        指定されたエネルギー、温度、フェルミ準位において、
        電子がそのエネルギー準位を占有する確率を計算します。
    引数:
        :param E: エネルギー (eV)。
        :type E: float
        :param T: 温度 (K)。
        :type T: float
        :param EF: フェルミ準位 (eV)。
        :type EF: float
    戻り値:
        :returns: 電子占有確率。
        :rtype: float
    """
    global e, kB

    if T == 0.0:
        if E < EF:
            return 1.0
        elif E == EF:
            return 0.5
        else:
            return 0.0

    k = (E - EF) * e / kB / T
    if k > nexp:
        k = nexp
    if k < -nexp:
        k = -nexp
    return 1.0 / (exp(k) + 1.0)

def fh(E: float, T: float, EF: float) -> float:
    """
    概要:
        フェルミ・ディラック分布関数に基づいて正孔の占有確率を計算します。
    詳細説明:
        電子占有確率を用いて、正孔の占有確率を計算します。
    引数:
        :param E: エネルギー (eV)。
        :type E: float
        :param T: 温度 (K)。
        :type T: float
        :param EF: フェルミ準位 (eV)。
        :type EF: float
    戻り値:
        :returns: 正孔占有確率。
        :rtype: float
    """
    return 1.0 - fe(E, T, EF)

# Define the function to be integrated
def DOSfe(E: float, T: float, EF: float) -> float:
    """
    概要:
        状態密度と電子のフェルミ・ディラック分布関数の積を計算します。
    詳細説明:
        これは電子密度を計算するための積分の被積分関数となります。
    引数:
        :param E: エネルギー (eV)。
        :type E: float
        :param T: 温度 (K)。
        :type T: float
        :param EF: フェルミ準位 (eV)。
        :type EF: float
    戻り値:
        :returns: 状態密度と電子占有確率の積の値。
        :rtype: float
    """
    return DOS(E) * fe(E, T, EF)

def DOSfh(E: float, T: float, EF: float) -> float:
    """
    概要:
        状態密度と正孔のフェルミ・ディラック分布関数の積を計算します。
    詳細説明:
        これは正孔密度を計算するための積分の被積分関数となります。
    引数:
        :param E: エネルギー (eV)。
        :type E: float
        :param T: 温度 (K)。
        :type T: float
        :param EF: フェルミ準位 (eV)。
        :type EF: float
    戻り値:
        :returns: 状態密度と正孔占有確率の積の値。
        :rtype: float
    """
    return DOS(E) * fh(E, T, EF)

def integrate_trapezoid(func: callable, E0: float, E1: float, h: float) -> float:
    """
    概要:
        関数を台形則で数値積分します。
    詳細説明:
        指定された関数funcを積分範囲[E0, E1]で刻み幅hを用いて台形則で数値積分します。
    引数:
        :param func: 積分する関数。引数としてエネルギーを取り、値を返します。
        :type func: callable
        :param E0: 積分範囲の開始エネルギー (eV)。
        :type E0: float
        :param E1: 積分範囲の終了エネルギー (eV)。
        :type E1: float
        :param h: 積分刻み幅 (eV)。
        :type h: float
    戻り値:
        :returns: 計算された積分の値。
        :rtype: float
    """
    n = int((E1 - E0) / h + 1.000001)
    y = [func(E0 + i * h) for i in range(n+1)]
    
    S = 0.5 * (y[0] + y[n]) + sum(y[1:n])
    return h * S


def integrate_trapezoid(func: callable, E0: float, E1: float, h: float) -> float:
    """
    概要:
        関数を台形則で数値積分します(再定義版)。
    詳細説明:
        指定された関数funcを積分範囲[E0, E1]で刻み幅hを用いて台形則で数値積分します。
        この関数は上の同名関数を上書きします。積分刻み幅は再計算されます。
    引数:
        :param func: 積分する関数。引数としてエネルギーを取り、値を返します。
        :type func: callable
        :param E0: 積分範囲の開始エネルギー (eV)。
        :type E0: float
        :param E1: 積分範囲の終了エネルギー (eV)。
        :type E1: float
        :param h: 積分刻み幅 (eV)。
        :type h: float
    戻り値:
        :returns: 計算された積分の値。
        :rtype: float
    """
    n = int((E1 - E0) / h + 1.000001)
    h = (E1 - E0) / n
    y = [func(E0 + i * h) for i in range(n+1)]
    
    S = 0.5 * (y[0] + y[n]) + sum(y[1:n])
    return h * S

def integrate_trapezoid_by_list(y: list[float], h: float) -> float:
    """
    概要:
        データ点のリストを台形則で数値積分します。
    詳細説明:
        予め計算された関数値のリストyと刻み幅hを用いて、台形則で数値積分を行います。
    引数:
        :param y: 積分する関数値のリスト。
        :type y: list
        :param h: データ間の刻み幅。
        :type h: float
    戻り値:
        :returns: 計算された積分の値。
        :rtype: float
    """
    n = len(y)
    S = 0.5 * (y[0] + y[n-1]) + sum(y[1:n-1])
    return h * S

def convert_range2index(E0: float, E1: float) -> tuple[int, int]:
    """
    概要:
        エネルギー範囲を状態密度データ配列のインデックス範囲に変換します。
    詳細説明:
        与えられたエネルギーに対応する状態密度データ配列の開始・終了インデックスを計算します。
    引数:
        :param E0: エネルギー範囲の開始値 (eV)。
        :type E0: float
        :param E1: エネルギー範囲の終了値 (eV)。
        :type E1: float
    戻り値:
        :returns: 開始インデックスと終了インデックスのタプル。
        :rtype: tuple
    """
    global Eraw_step, Eraw_min, Eraw_max

    i0 = (int)((E0 - Enorm_min) / Enorm_step)
    i1 = (int)((E1 - Enorm_min) / Enorm_step + 0.99999)
#    print("E0, E1, i0, i1 = ", E0, E1, i0, i1)
#    exit()
    
    return i0, i1

def Ne(T: float, EF: float, E0: float, E1: float) -> float:
    """
    概要:
        指定されたエネルギー範囲における電子密度を計算します。
    詳細説明:
        温度とフェルミ準位における電子のDOSfe関数を積分範囲で数値積分することにより
        電子密度を計算します。
    引数:
        :param T: 温度 (K)。
        :type T: float
        :param EF: フェルミ準位 (eV)。
        :type EF: float
        :param E0: 積分範囲の開始エネルギー (eV)。
        :type E0: float
        :param E1: 積分範囲の終了エネルギー (eV)。
        :type E1: float
    戻り値:
        :returns: 計算された電子密度 (cm^-3)。
        :rtype: float
    """
    global Enorm_step, Enorm_min, Enorm_max
    global integrator 

    if integrator == 'interpolate':
        N = integrate_trapezoid(lambda E: DOSfe(E, T, EF), E0, E1, Estep)
    else:
        iEinteg0, iEinteg1 = convert_range2index(E0, E1)
        flist = [dos_raw[iE] * fe(Emin + iE * Estep, T, EF) for iE in range(iEinteg0, iEinteg1 + 1)]
        N = integrate_trapezoid_by_list(flist, Estep)

    return N

def Nh(T: float, EF: float, E0: float, E1: float) -> float:
    """
    概要:
        指定されたエネルギー範囲における正孔密度を計算します。
    詳細説明:
        温度とフェルミ準位における正孔のDOSfh関数を積分範囲で数値積分することにより
        正孔密度を計算します。
    引数:
        :param T: 温度 (K)。
        :type T: float
        :param EF: フェルミ準位 (eV)。
        :type EF: float
        :param E0: 積分範囲の開始エネルギー (eV)。
        :type E0: float
        :param E1: 積分範囲の終了エネルギー (eV)。
        :type E1: float
    戻り値:
        :returns: 計算された正孔密度 (cm^-3)。
        :rtype: float
    """
    global Estep, Emin
    global integrator 
    
    if integrator == 'interpolate':
        N = integrate_trapezoid(lambda E: DOSfh(E, T, EF), E0, E1, Estep)
    else:
        iEinteg0, iEinteg1 = convert_range2index(E0, E1)
        flist = [dos_raw[iE] * fh(Emin + iE * Estep, T, EF) for iE in range(iEinteg0, iEinteg1 + 1)]
        N = integrate_trapezoid_by_list(flist, Estep)

    return N


def FindBandEdges(E: list[float], DOS: list[float], Emidgap: float, DOSth: float) -> tuple[float, float]:
    """
    概要:
        状態密度データから価電子帯上端と伝導帯下端を探索します。
    詳細説明:
        状態密度データと対応するエネルギーを基に、バンドギャップ中央から
        状態密度の閾値を超える点を探索し、価電子帯上端と伝導帯下端を決定します。
    引数:
        :param E: エネルギーのリスト (eV)。
        :type E: list
        :param DOS: 状態密度のリスト。
        :type DOS: list
        :param Emidgap: バンドギャップの中央エネルギーの推定値 (eV)。
        :type Emidgap: float
        :param DOSth: 状態密度がこの値を超えたときにバンド端と判定する閾値。
        :type DOSth: float
    戻り値:
        :returns: 価電子帯上端と伝導帯下端のエネルギーのタプル。
        :rtype: tuple
    """
    EV = 1.0e30
    EC = 1.0e30
    i0 = -1
    for i in range(len(E)):
        if E[i] >= Emidgap:
            i0 = i
            break

    for i in range(i0, 0, -1):
        if DOS[i] > DOSth:
            EV = E[i+1]
            break

    for i in range(i0, len(E), +1):
        if DOS[i] > DOSth:
            EC = E[i-1]
            break

    return EV, EC

def diff(h: float, T: float, EF: float, ECmin: float, ECmax: float, EVmin: float, EVmax: float, N0: float) -> float:
    """
    概要:
        全電荷中性条件からの偏差(dQ)のフェルミ準位による数値微分を計算します(レガシー版)。
    詳細説明:
        この関数は現在使用されていないようです。tkDefectsクラスのdiffメソッドが使用されています。
        dQ関数が未定義のため、この関数は実行されません。
    引数:
        :param h: 微小なEFの摂動量 (eV)。
        :type h: float
        :param T: 電子温度 (K)。
        :type T: float
        :param EF: フェルミ準位 (eV)。
        :type EF: float
        :param ECmin: 伝導帯のエネルギー下限 (eV)。
        :type ECmin: float
        :param ECmax: 伝導帯のエネルギー上限 (eV)。
        :type ECmax: float
        :param EVmin: 価電子帯のエネルギー下限 (eV)。
        :type EVmin: float
        :param EVmax: 価電子帯のエネルギー上限 (eV)。
        :type EVmax: float
        :param N0: (現在未使用)。
        :type N0: float
    戻り値:
        :returns: dQのEFによる数値微分。
        :rtype: float
    """
    # dQ関数が定義されていないため、この関数は現状では機能しない
    # ここでは仮の戻り値として0を返すか、エラーを出すべきだが、既存コード変更禁止のためそのまま
    return 0.0 # (dQ(T, EF + h, ECmin, ECmax, EVmin, EVmax, N0) - dQ(T, EF - h, ECmin, ECmax, EVmin, EVmax, N0)) / 2.0 / h
    
def callbackfunc(obj: object):
    """
    概要:
        数値計算ソルバーの進行状況を表示するためのコールバック関数。
    詳細説明:
        ソルバーの反復回数、現在の解、関数値、およびステップサイズをコンソールに出力します。
    引数:
        :param obj: ソルバーオブジェクトのインスタンス。
        :type obj: object
    戻り値:
        :returns: 常に1を返します。
        :rtype: int
    """
    if obj.xa is None:
        print(f"Iter {obj.iter:5d}: x: {obj.x:>12.6g} f: {obj.f:>12.6g}, dx = {obj.dx:>10.4g}")
    else:
        print(f"Iter {obj.iter:5d}: x: {obj.x:>12.6g} f: {obj.f:>12.6g} in [{obj.xa:>12.6g}, {obj.xb:>12.6g}], dx = {obj.dx:>10.4g}")
    return 1


def read_files(vasp: 'tkVASP', cry: 'tkCrystal', cry_defect: 'tkCrystal', defects: 'tkDefects', POSCAR_path: str, OUTCAR_path: str, EIGENVAL_path: str, DOSCAR_path: str, fp: object = outfp) -> tuple[dict, dict, dict, dict, float, float, float, float, float]:
    """
    概要:
        VASPの計算結果ファイルを読み込み、関連情報を設定します。
    詳細説明:
        結晶構造、セル体積、フェルミ準位、バンド端、状態密度データを読み込み、
        それらをグローバル変数に設定します。また、POSCARと欠陥データのサイト数の一致をチェックします。
        戻り値はOUTCAR情報辞書、EIGENVAL情報辞書、バンド端情報辞書、状態密度情報辞書、
        価電子帯上端、伝導帯下端、バンドギャップ、HOMO準位、LUMO準位のタプルです。
    引数:
        :param vasp: tkVASPオブジェクトのインスタンス。
        :type vasp: tkVASP
        :param cry: 完全結晶のtkCrystalオブジェクトのインスタンス。
        :type cry: tkCrystal
        :param cry_defect: 欠陥結晶のtkCrystalオブジェクトのインスタンス。
        :type cry_defect: tkCrystal
        :param defects: tkDefectsオブジェクトのインスタンス。
        :type defects: tkDefects
        :param POSCAR_path: 完全結晶のPOSCARファイルのパス。
        :type POSCAR_path: str
        :param OUTCAR_path: OUTCARファイルのパス。
        :type OUTCAR_path: str
        :param EIGENVAL_path: EIGENVALファイルのパス。
        :type EIGENVAL_path: str
        :param DOSCAR_path: DOSCARファイルのパス。
        :type DOSCAR_path: str
        :param fp: 出力先のファイルポインタ。
        :type fp: object | None
    戻り値:
        :returns: VASPファイルから読み込んだ情報のタプル。
        :rtype: tuple
    """
    global EF0
    global Vcell
    global E_raw, E_norm, dos_raw, nDOS, Enorm_min, Enorm_max, Enorm_step, Emin, Emax, Estep
    global EFmin, EFmax, nEF, EFstep
    global EV, EC
    global fdos

    mprint("", fp = outfp)
    mprint("Crystal structure from [{}]".format(POSCAR_path), fp = outfp)
    a, b, c, alpha, beta, gamm = cry.LatticeParameters()
    Vcell = cry.Volume()
#    cry.PrintInf("cell")
    mprint("  cell: {:12.8f} {:12.8f} {:12.8f} A   {:10.6f} {:10.6f} {:10.6f}".format(a, b, c, alpha, beta, gamm), fp = outfp)
    mprint("  volume: {:12.6f} A^-3".format(Vcell), fp = outfp)

    ret = check_atom_sites(cry_defect, defects, fp = outfp)
#    ret = check_atom_sites(cry, defects, fp = outfp)
    if not ret:
        mprint("", fp = outfp)
        app.terminate("", usage = usage, pause = True)

    mprint("", fp = outfp)
    mprint("Read [{}]".format(OUTCAR_path), fp = outfp)
    outcarinf = vasp.read_outcar_inf(OUTCAR_path)
    EF0   = outcarinf["EF"]
    ISPIN = outcarinf["ISPIN"]
    mprint("Information in OUTCAR:", fp = outfp)
    mprint("  ISPIN: ", ISPIN, fp = outfp)
    mprint("  EF = {} eV corrected to 0".format(EF0), fp = outfp)
    if stop_spin_polarized and ISPIN == 2:
        mprint("", fp = outfp)
        mprint("Error: ISPIN=2 is not implemented", fp = outfp)
        mprint("", fp = outfp)
        app.terminate("", usage = usage, pause = True)

    mprint("", fp = outfp)
    mprint("Read EV, EC, Eg with EF0={} eV from [{}]".format(EF0, EIGENVAL_path), fp = outfp)
    eigenvalinf = vasp.read_eigenval(EIGENVAL_path, EF = 0.0)
    nk      = eigenvalinf["nk"]
    nLevels = eigenvalinf["nLevels"]
    mprint("k points in EIGENVAL:", fp = outfp)
    print("nk=", nk)
    print("nLevels=", nLevels)

    bandedgeinf = vasp.find_band_edges_from_eigenval(EF0 = EF0, eigenvalinf = eigenvalinf, ISPIN = ISPIN)
    EV    = bandedgeinf["EV"]
    EC    = bandedgeinf["EC"]
    Eg    = bandedgeinf["Eg"]
    EHOMO = bandedgeinf["EHOMO"]
    ELUMO = bandedgeinf["ELUMO"]

    mprint("", fp = outfp)
    mprint("*** Read Total DOS from [{}]".format(DOSCAR_path), fp = outfp)
# E_raw is measured from EF in OUTCAR
#    E_raw, dos_raw = np.genfromtxt(DOSCAR_path, skip_header=6, usecols=(0,1), unpack=True)
#    nDOS = len(E_raw)
    dosinf = vasp.read_doscar(DOSCAR_path, cry, IsSpinPolarized = ISPIN, IsNonCollinear = None, EF = 0.0)
    E_raw   = dosinf["E"]
    dos_raw = dosinf["TotalDOS"]
    nDOS    = dosinf["nE"]
    for i in range(len(E_raw)):
        dos_raw[i] *= 1.0 / (Vcell * 1.0e-24)
    Emin  = E_raw[0]
    Emax  = E_raw[nDOS-1]
    Estep = (E_raw[nDOS-1] - E_raw[0]) / (nDOS - 1)
    mprint("  DOS E range                     : {} - {}, {} eV step".format(Emin, Emax, Estep), fp = outfp)

# E_norm is measured from EV
    E_norm = []
    for i in range(len(E_raw)):
        E_norm.append(E_raw[i] - EV)
    Enorm_min  = E_norm[0]
    Enorm_max  = E_norm[nDOS-1]
    Enorm_step = (E_norm[nDOS-1] - E_norm[0]) / (nDOS - 1)
    Emin  -= EV
    Emax  -= EV
    mprint("  DOS E range normalized to EV = 0: {} - {}, {} eV step".format(Emin, Emax, Estep), fp = outfp)
#    fdos = interp1d(E_norm, dos_raw, kind = 'cubic')
    fdos = interp1d(E_norm, dos_raw, kind = 'linear')

    mprint("", fp = outfp)
    mprint("*** Find band edges from eigenvalinf", fp = outfp)
    mprint("Band edges from EIGENVAL       : EV={:10.6f}  EC={:10.6f}  Eg={:10.6f}  EF0={:10.6f} eV"
        .format(EV, EC, Eg, EF0), fp = outfp)
    EF0   -= EV
    EC    -= EV
    EHOMO -= EV
    ELUMO -= EV
    EV = 0.0
    mprint("Band edges normalized by EV = 0: EV={:10.6f}  EC={:10.6f}  Eg={:10.6f}  EF0={:10.6f} eV"
        .format(EV, EC, Eg, EF0), fp = outfp)

    return outcarinf, eigenvalinf, bandedgeinf, dosinf, EV, EC, Eg, EHOMO, ELUMO

def check_atom_sites(cry: 'tkCrystal', defects: 'tkDefects', fp: object = None) -> bool:
    """
    概要:
        POSCARファイルと欠陥形成エネルギーのExcelファイルで指定された原子サイト数が一致するかをチェックします。
    詳細説明:
        POSCARから原子サイト数を取得し、Excelファイルから読み込まれた欠陥データの基準サイト数と比較します。
        不一致がある場合、ユーザーに続行するか終了するかを問い合わせる警告メッセージを表示します。
    引数:
        :param cry: 結晶構造情報を持つtkCrystalオブジェクト。
        :type cry: tkCrystal
        :param defects: 欠陥情報を持つtkDefectsオブジェクト。
        :type defects: tkDefects
        :param fp: 出力先のファイルポインタ。
        :type fp: object | None
    戻り値:
        :returns: サイト数チェックを続行する場合True、終了する場合False。
        :rtype: bool
    """
# Check atom sites
    mprint("", fp = outfp)
    mprint("Check atom sites", fp = outfp)
    mprint("  from POSCAR:", fp = outfp)
#    AtomSites = cry.ExpandedAtomSiteList()
#    for atom in AtomSites:
#        label     = atom.Label()
#        pos       = atom.Position()
#        print("  {} ({}, {}, {})".format(label, *pos))
    AtomTypes = cry.AtomTypeList()
    nsites = {}
    for t in AtomTypes:
        type = t.AtomTypeOnly()
        nsites[type] = cry.count_by_type(type, mode = 'short', target = 'expanded')
    
    for type in nsites.keys():
        mprint("    {:4} {:4d}".format(type, nsites[type]))
    mprint("  from {}:".format(dH_path), fp = outfp)
    checked = {}
    is_error = ''
    site_names = list(nsites.keys())
    for d in defects.defects:
        if d.site not in checked.keys():
            mprint("    {:4} {:4f}".format(d.site, d.N0))
            site_name = d.site
            checked[site_name] = 1
            if d.site not in site_names:
                m = re.match(r'([A-Z][a-z]?)', d.site)
                if m:
#                    d.site = m.groups()[0]
                    site_name = m.groups()[0]

            if site_name in nsites.keys():
                ns = nsites[site_name]
#                ns = nsites[d.site]
            else:
                ns = 0
            
            if abs(d.N0 - ns) > 1.0e-3:
                is_error += "{} ({} != {}), ".format(d.site, d.N0, ns)

    if not ignore_warning and is_error != '':
        mprint("", fp = outfp)
        mprint("========================================================================", fp = outfp)
        mprint(f"Error: Site numbers mismatch between POSCAR and {dH_path}", fp = outfp)
        mprint(f"   {is_error}")
        mprint("========================================================================", fp = outfp)
        print("")
        print(f"Choose continue: Enter continue to disregard this error")
        print(f"Choose stop    : Enter stop to terminate this run, and correct site numbers in {dH_path}")
        print(f"                 so as to corresponds to POSCAR")
        while 1:
            print("  coninue or stop>>", end = '')
            answer = input()
            if answer == 'continue':
                break
            elif answer == 'stop':
                return False
    return True

def print_transition_level(defects: 'tkDefects') -> tuple[str, list[list[list]], list[list[list]]]:
    """
    概要:
        欠陥の遷移準位を計算し、出力ファイルに保存して表示します。
    詳細説明:
        tkDefectsオブジェクトから遷移準位のリストを取得し、
        すべての形成エンタルピー曲線と最小形成エンタルピー曲線をExcelファイルに保存します。
        戻り値は最後の欠陥名、全遷移点リスト、最小曲線上の点リストのタプルです。
    引数:
        :param defects: 欠陥情報を持つtkDefectsオブジェクト。
        :type defects: tkDefects
    戻り値:
        :returns: 欠陥名と遷移点リストのタプル。
        :rtype: tuple
    """
#==============================================
# Transition levels
#==============================================
    print("")
    print("Find transition levels")
    dHEF_all_xlsx = modify_path(dH_path, '-dH-EF-all-{}.xlsx'.format(defects.plabels[iPoint]))
    dHEF_min_xlsx = modify_path(dH_path, '-dH-EF-min-{}.xlsx'.format(defects.plabels[iPoint]))
    print("  Remove [{}]".format(dHEF_all_xlsx))
    delete_file(dHEF_all_xlsx)
    defects.save_allpts(dHEF_all_xlsx)

    print("  Remove [{}]".format(dHEF_min_xlsx))
    delete_file(dHEF_min_xlsx)
    minwb = tkExcel(dHEF_min_xlsx, 'w')
    defects.save_minpts(dHEF_min_xlsx)

    names, allpts_list, minpts_list = defects.get_transitionlevel_lists()

    for iname in range(len(names)):
        name = names[iname]
#        groupdata = defects.get_groupdata(name, iPoint)
        allpts = allpts_list[iname]
        minpts = minpts_list[iname]

        mprint("", fp = outfp)
        for ipt in range(len(allpts)):
            if abs(allpts[ipt][2] - allpts[ipt][3]) > 1.0e-3:
                mprint("  ***all point: {:6}({:4}/{:4}) EF={:10.4f} eV  dH={:10.4f} eV  iorder={} idx={}"
                              .format(name, allpts[ipt][2], allpts[ipt][3], 
                                        allpts[ipt][0], allpts[ipt][1], allpts[ipt][4], allpts[ipt][5]), fp = outfp)
        mprint("", fp = outfp)
        for ipt in range(len(minpts)):
            if abs(minpts[ipt][2] - minpts[ipt][3]) > 1.0e-3:
                mprint("  ***min point: {:6}({:4}/{:4}) EF={:10.4f} eV  dH={:10.4f} eV  iorder={} idx={}"
                                .format(name, minpts[ipt][2], minpts[ipt][3], minpts[ipt][0], minpts[ipt][1], 
                                        minpts[ipt][4], minpts[ipt][5]), fp = outfp)

    return name, allpts_list, minpts_list

def exec_T():
    """
    概要:
        温度依存性の計算とプロットを実行します。
    詳細説明:
        欠陥形成エネルギーのExcelファイルとVASP計算結果ファイルを読み込み、
        様々な温度における電子・正孔密度、欠陥密度、およびフェルミ準位を計算します。
        計算結果はExcelファイルに保存され、グラフとして可視化されます。
    """
    global mode, plot_mode, T0, EF0
    global dH_path, CAR_path, CAR_path_defect, DOSCAR_path, Vcell
    global E_raw, E_norm, dos_raw, nDOS, Enorm_min, Enorm_max, Enorm_step, Emin, Emax, Estep
    global EV, EC
    global EFmin, EFmax, nEF, EFstep
    global Einteg0
    global nmaxiter_bisection, eps_bisec
    global nmaxiter_newton, eps_newton
    global fdos
    global outfp # グローバル変数としてoutfpを宣言

    kBTe    = kB * T0 / e
    kBTedef = kB * Tdef / e
    if T0 < Tdef:
        dE = nrange * kBTedef
    else:
        dE = nrange * kBTe

    vasp = tkVASP()
    CAR_path = vasp.getdir(CAR_path)
    INCAR_path   = vasp.get_INCAR(CAR_path)
    POSCAR_path  = vasp.get_POSCAR(CAR_path)
    CONTCAR_path = vasp.get_CONTCAR(CAR_path)
    OUTCAR_path  = vasp.get_OUTCAR(CAR_path)
    EIGENVAL_path = vasp.get_VASPPath(CAR_path, 'EIGENVAL')
    DOSCAR_path   = vasp.get_VASPPath(CAR_path, 'DOSCAR')

    CAR_path_defect = vasp.getdir(CAR_path_defect)
    POSCAR_path_defect  = vasp.get_POSCAR(CAR_path_defect)

    print("")
    print("*** Read crystal structure from [{}]".format(POSCAR_path))
    cry = vasp.read_poscar(POSCAR_path)
#    cry.PrintInf("cell")
    Vcell = cry.Volume()
    mprint("  volume: {:12.6f} A^-3".format(Vcell))

    mprint("Crystal structure of defect model from [{}]".format(POSCAR_path_defect))
    cry_d = vasp.read_poscar(POSCAR_path_defect)
#    a_d, b_d, c_d, alpha_d, beta_d, gamm_d = cry_d.LatticeParameters()
    Vcell_d = cry_d.Volume()
#    cry.PrintInf("cell")
#    mprint("  cell: {:12.8f} {:12.8f} {:12.8f} A   {:10.6f} {:10.6f} {:10.6f}".format(a_d, b_d, c_d, alpha_d, beta_d, gamm_d), fp = outfp)
    mprint("  volume: {:12.6f} A^-3".format(Vcell_d))

    print("")
    print("Read defect formation enthalpies from [{}]".format(dH_path))
    defects = tkDefects()
#    defects.SetVolume(Vcell)
    defects.read_excel(dH_path, Vcell * int(Vcell_d / Vcell + 0.2))
    defects.Print()

    print("")
    outputfile    = modify_path(dH_path, '-out-{}.txt'.format(defects.plabels[iPoint]))
    print("  Remove [{}]".format(outputfile))
    delete_file(outputfile)
    T_xlsx = modify_path(dH_path, '-T-{}.xlsx'.format(defects.plabels[iPoint]))
    delete_file(T_xlsx)

    print("")
    print("*** Open output file {}".format(outputfile))
    outfp = open(outputfile, 'w')
    if outfp is None:
        app.terminate("Error: Can not write to [{}]".format(outputfile), pause = True)

    mprint("", fp = outfp)
    mprint("=======================================================", fp = outfp)
    mprint(" Calculate T dependence of semiconductor properties", fp = outfp)
    mprint("=======================================================", fp = outfp)
    mprint("mode     : ", mode, fp = outfp)

    mprint("", fp = outfp)
    mprint("Files:", fp = outfp)
    mprint("  dH input  : ", dH_path, fp = outfp)
    mprint("", fp = outfp)
    mprint("Files:", fp = outfp)
    mprint("  dH input  : ", dH_path, fp = outfp)
    mprint("  Output    : ", outputfile, fp = outfp)

    mprint("", fp = outfp)
    mprint("Defects red from [{}]".format(dH_path), fp = outfp)
    defects.Print(outfp = outfp) 
    mprint("", fp = outfp)
    mprint("To be analyzed for Point {}".format(defects.plabels[iPoint]), fp = outfp)
    mprint("T(defects): {} K".format(Tdef))

    mprint("", fp = outfp)
    mprint("Temperature range: {} - {} K, {} points".format(Tmin, Tmax, nT), fp = outfp)

    mprint("", fp = outfp)
    mprint("VASP files  :", fp = outfp)
    mprint("  CAR dir   : ", CAR_path, fp = outfp)
    mprint("  INCAR     : ", INCAR_path, fp = outfp)
    mprint("  POSCAR    : ", POSCAR_path, fp = outfp)
    mprint("  CONTCAR   : ", CONTCAR_path, fp = outfp)
    mprint("  OUTCAR    : ", OUTCAR_path, fp = outfp)
    mprint("  EIGENVAL  : ", EIGENVAL_path)
    mprint("  DOSCAR    : ", DOSCAR_path)
    mprint("  POSCAR(defect): ", POSCAR_path_defect, fp = outfp)

    outcarinf, eigenvalinf, bandedgeinf, dosinf, EV, EC, Eg, EHOMO, ELUMO = \
        read_files(vasp, cry, cry_d, defects, POSCAR_path, OUTCAR_path, EIGENVAL_path, DOSCAR_path, fp = outfp)

    EFmin = EV + dEFmin
    EFmax = EC + dEFmax
    EFstep = (EFmax - EFmin) / (nEF - 1)

    mprint("", fp = outfp)
    mprint("EF range to save to transition level files", fp = outfp)
    mprint("  EF range: {} - {}, {} eV step".format(EFmin, EFmax, EFstep), fp = outfp)

    mprint("", fp = outfp)
    mprint("T dependence", fp = outfp)
    mprint("  T range: {} - {}, {} K step".format(Tmin, Tmax, Tstep), fp = outfp)

    mprint("", fp = outfp)
    mprint("Integration configuration", fp = outfp)
    mprint("  Integration E range: {} - {} eV, or EF+-{}*kBT eV".format(Einteg0, Einteg1, nrange), fp = outfp)

    ECmin = EC - Estep
    ECmax = EC + dE
    EVmin = EV - dE
    EVmax = EV + Estep

#==============================================
# Print transition levels
#==============================================
    name, allpts_list, minpts_list = print_transition_level(defects)

#==============================================
# Start calculations
#==============================================
# At Tdef
    EFmin_ini = EV + dEFmin
    EFmax_ini = EC + dEFmax
    mprint("", fp = outfp)
    mprint("Calculate defect densities at T(defects) = {} K".format(Tdef))
    mprint("  Calculate EF,eq: Bisecton method: eps = {}, nmaxiter = {}".format(eps_bisec, nmaxiter_bisection), fp = outfp)
    mprint("                        Initial EF range: {:12.4f} - {:12.4f} eV".format(EFmin_ini, EFmax_ini))
    mprint("                   Newton method  : eps = {}, nmaxiter = {}, dump = {}".format(eps_newton, nmaxiter_newton, dump_newton), fp = outfp)
    if EFdef == 'eq':
        EFeqTdef, dQ = defects.find_EF(Tdef, Tdef, None, ECmin, ECmax, EVmin, EVmax,
                    callbackfunc = callbackfunc, 
                    eps_bisec  = eps_bisec,  nmaxiter_bisection = nmaxiter_bisection, initial = [EFmin_ini, EFmax_ini],
                    eps_newton = eps_newton, nmaxiter_newton = nmaxiter_newton, dump_newton = dump_newton,
                    h_newton = h_newton
                    )
        if EFeqTdef is None:
            app.terminate("Error in tkDefects.find_EF: Could not reach convergence for EF", pause = True)
    else:
        try:
            EFeqTdef = float(EFdef)
        except:
            app.terminate("Error in tkDefects.find_EF: EFdef [{}] must be numeral".format(EFdef), pause = True)

    mprint("", fp = outfp)
    mprint("  EF,eq({:8.3f} K)={:12.6g} eV".format(Tdef, EFeqTdef), fp = outfp)
    Z, ne, nh, Nds, NdsTdef, dGs, Qtot = defects.calculate_densities(T0, Tdef, None, EFeqTdef, ECmin, ECmax, EVmin, EVmax)
    for id in range(defects.ndefects):
        d = defects.defects[id]
        label  = "{}_{}".format(d.atom, d.site)
        labelq = "{}^{}".format(label, d.charge)
        mprint("    {:>16} (N0={:5g} Ndoped={:5g}):".format(labelq, d.N0, d.Ndoped), end = '', fp = outfp)
        mprint("       Nds={:12.4g} cm-3  dGs={:12.4g} eV".format(Nds[id], dGs[id]), fp = outfp)
    mprint("", fp = outfp)
    mprint("  Total defect densities at {} K".format(Tdef), fp = outfp)
    for key in NdsTdef.keys():
        Nsite = NdsTdef[key] * defects.V * 1.0e-24
        mprint("    {:>8}: {:12.6g} cm-3 ({:8.2f} sites)".format(key, NdsTdef[key], Nsite), fp = outfp)

# At T0, EF,eq(T0)
    EFmin_ini = EV + dEFmin
    EFmax_ini = EC + dEFmax
    mprint("", fp = outfp)
    if T0 <= Tdef:
        mprint("Calculate defect densities at T(electron) = {} K".format(T0))
        mprint("  Total defect densities frozen at {} K with EF = {:12.4f} eV".format(Tdef, EFeqTdef), fp = outfp)
        for key in NdsTdef.keys():
            Nsite = NdsTdef[key] * defects.V * 1.0e-24
            mprint("    {:>8}: {:12.6g} cm-3 ({:8.2f} sites)".format(key, NdsTdef[key], Nsite), fp = outfp)
    else:
        mprint("T(electron) is higher than T(defects). Defect densities are calculated at T(electron)")
        mprint("Calculate defect densities at T(electron) = {} K, T(defects) = {} K".format(T0, T0))
    mprint("  Calculate EF,eq: Bisecton method: eps = {}, nmaxiter = {}".format(eps_bisec, nmaxiter_bisection), fp = outfp)
    mprint("                        Initial EF range: {:12.4f} - {:12.4f} eV".format(EFmin_ini, EFmax_ini))
    mprint("                   Newton method  : eps = {}, nmaxiter = {}, dump = {}".format(eps_newton, nmaxiter_newton, dump_newton), fp = outfp)
    EFeqT0, dQ = defects.find_EF(T0, Tdef, NdsTdef, ECmin, ECmax, EVmin, EVmax,
                    callbackfunc = callbackfunc, 
                    eps_bisec  = eps_bisec,  nmaxiter_bisection = nmaxiter_bisection, initial = [EFmin_ini, EFmax_ini],
                    eps_newton = eps_newton, nmaxiter_newton = nmaxiter_newton, dump_newton = dump_newton,
                    h_newton = h_newton
                    )
    if EFeqT0 is None:
        app.terminate("Error in tkDefects.find_EF: Could not reach convergence for EF", pause = True)

#==============================================
# Start calculations
#==============================================
# At Tdef
    EFmin_ini = EV + dEFmin
    EFmax_ini = EC + dEFmax
    mprint("", fp = outfp)
    mprint("Calculate defect densities at T(electron) = T(defects) = {} K".format(Tdef))
    mprint("  Calculate EF,eq: Bisecton method: eps = {}, nmaxiter = {}".format(eps_bisec, nmaxiter_bisection), fp = outfp)
    mprint("                        Initial EF range: {:12.4f} - {:12.4f} eV".format(EFmin_ini, EFmax_ini))
    mprint("                   Newton method  : eps = {}, nmaxiter = {}, dump = {}".format(eps_newton, nmaxiter_newton, dump_newton), fp = outfp)
    if EFdef == 'eq':
        EFeqTdef, dQ = defects.find_EF(Tdef, Tdef, None, ECmin, ECmax, EVmin, EVmax,
                    callbackfunc = callbackfunc, 
                    eps_bisec  = eps_bisec,  nmaxiter_bisection = nmaxiter_bisection, initial = [EFmin_ini, EFmax_ini],
                    eps_newton = eps_newton, nmaxiter_newton = nmaxiter_newton, dump_newton = dump_newton,
                    h_newton = h_newton
                    )
        if EFeqTdef is None:
            app.terminate("Error in tkDefects.find_EF: Could not reach convergence for EF", pause = True)
    else:
        try:
            EFeqTdef = float(EFdef)
        except:
            app.terminate("Error in tkDefects.find_EF: EFdef [{}] must be numeral".format(EFdef), pause = True)

    mprint("", fp = outfp)
    mprint("  EF,eq({:8.3f} K)={:12.6g} eV".format(Tdef, EFeqTdef), fp = outfp)
    Z, ne, nh, Nds, NdsTdef, dGs, Qtot = defects.calculate_densities(T0, Tdef, None, EFeqTdef, ECmin, ECmax, EVmin, EVmax)
    for id in range(defects.ndefects):
        d = defects.defects[id]
        label  = "{}_{}".format(d.atom, d.site)
        labelq = "{}^{}".format(label, d.charge)
        mprint("    {:>16} (N0={:5g} Ndoped={:5g}):".format(labelq, d.N0, d.Ndoped), end = '', fp = outfp)
        mprint("       Nds={:12.4g} cm-3  dGs={:12.4g} eV".format(Nds[id], dGs[id]), fp = outfp)
    mprint("", fp = outfp)
    mprint("  Total defect densities at {} K".format(Tdef), fp = outfp)
    for key in NdsTdef.keys():
        Nsite = NdsTdef[key] * defects.V * 1.0e-24
        mprint("    {:>8}: {:12.6g} cm-3 ({:8.2f} sites)".format(key, NdsTdef[key], Nsite), fp = outfp)

# At T0, EF,eq(T0)
    EFmin_ini = EV + dEFmin
    EFmax_ini = EC + dEFmax
    mprint("", fp = outfp)
    mprint("Calculate defect densities at T(electron) = {} K, T(defects) = {} K".format(T0, Tdef))
    mprint("  Calculate EF,eq: Bisecton method: eps = {}, nmaxiter = {}".format(eps_bisec, nmaxiter_bisection), fp = outfp)
    mprint("                        Initial EF range: {:12.4f} - {:12.4f} eV".format(EFmin_ini, EFmax_ini))
    mprint("                   Newton method  : eps = {}, nmaxiter = {}, dump = {}".format(eps_newton, nmaxiter_newton, dump_newton), fp = outfp)
    EFeqT0, dQ = defects.find_EF(T0, Tdef, NdsTdef, ECmin, ECmax, EVmin, EVmax,
                    callbackfunc = callbackfunc, 
                    eps_bisec  = eps_bisec,  nmaxiter_bisection = nmaxiter_bisection, initial = [EFmin_ini, EFmax_ini],
                    eps_newton = eps_newton, nmaxiter_newton = nmaxiter_newton, dump_newton = dump_newton,
                    h_newton = h_newton
                    )
    if EFeqT0 is None:
        app.terminate("Error in tkDefects.find_EF: Could not reach convergence for EF", pause = True)

    mprint("", fp = outfp)
    mprint("  EF,eq({:8.3g} K)={:12.6g} eV".format(T0, EFeqT0), fp = outfp)
    Z, ne, nh, Nds, NdsTdef, dGs, Qtot = defects.calculate_densities(T0, Tdef, NdsTdef, EFeqT0, ECmin, ECmax, EVmin, EVmax)
    for id in range(defects.ndefects):
        d = defects.defects[id]
        label  = f"{d.atom}_{d.site}"
        labelq = f"{label}^{d.charge}"
        mprint(f"    {labelq:>16} (Nds(tot)={NdsTdef[labelq]:12.4g} cm-3):", end = '', fp = outfp)
#        mprint(f"    {labelq:>16} (Nds(tot)={NdsTdef[label]:12.6g} cm-3):", end = '', fp = outfp)
        mprint(f"       Nds={Nds[id]:12.4g} cm-3  dGs={dGs[id]:12.4g} eV", fp = outfp)

# At Tmin, EF,eq(Tmin)
    EFmin_ini = EV + dEFmin
    EFmax_ini = EC + dEFmax
    mprint("", fp = outfp)
    mprint("Calculate defect densities at T(min) = {} K, T(defects) = {} K".format(Tmin, Tdef))
    mprint("  Calculate EF,eq: Bisecton method: eps = {}, nmaxiter = {}".format(eps_bisec, nmaxiter_bisection), fp = outfp)
    mprint("                        Initial EF range: {:12.4f} - {:12.4f} eV".format(EFmin_ini, EFmax_ini))
    mprint("                   Newton method  : eps = {}, nmaxiter = {}, dump = {}".format(eps_newton, nmaxiter_newton, dump_newton), fp = outfp)
    EFeqTmin, dQ = defects.find_EF(Tmin, Tdef, NdsTdef, ECmin, ECmax, EVmin, EVmax,
                    callbackfunc = callbackfunc, 
                    eps_bisec  = eps_bisec,  nmaxiter_bisection = nmaxiter_bisection, initial = [EFmin_ini, EFmax_ini],
                    eps_newton = eps_newton, nmaxiter_newton = nmaxiter_newton, dump_newton = dump_newton,
                    h_newton = h_newton
                    )

    if EFeqTmin is None:
        app.terminate("Error in tkDefects.find_EF: Could not reach convergence for EF at {} K".format(EFeqTmin), pause = True)
    else:
        mprint("  EF,eq at {} K = {:12.6g} eV".format(Tmin, EFeqTmin), fp = outfp)

# At various T, EF,eq(T)
    mprint("", fp = outfp)
    mprint("Calculate defect densities at various T = {} - {} K".format(Tmin, Tmax))
    mprint("  Note if T(electron) is higher than T(defects), defects are not frozen.")
    mprint("  If T is lower than T(defects) = {} K, ".format(Tdef), fp = outfp)
    mprint("  total defect densities are frozen at {} K with EF = {:.4f} eV as follows.".format(Tdef, EFeqTdef), fp = outfp)
    for key in NdsTdef.keys():
        Nsite = NdsTdef[key] * defects.V * 1.0e-24
        mprint("    {:>8}: {:12.6g} cm-3 ({:8.2f} sites)".format(key, NdsTdef[key], Nsite), fp = outfp)

    mprint("  Calculate EF,eq: Bisecton method: eps = {}, nmaxiter = {}".format(eps_bisec, nmaxiter_bisection), fp = outfp)
    mprint("                   Newton method  : eps = {}, nmaxiter = {}, dump = {}".format(eps_newton, nmaxiter_newton, dump_newton), fp = outfp)
    EF = EFeqTmin
    xT = [Tmin + i * Tstep for i in range(nT)]
    xT1000 = []
    yEF    = []
    ydQ    = []
    yne  = []
    ynh  = []
    ynds = []
    for id in range(defects.ndefects):
        ynds.append([0.0]*nT)
    EFmin_ini = EV + dEFmin
    EFmax_ini = EC + dEFmax
    for iT in range(nT):
        T = xT[iT]
        xT1000.append(1000.0 / T)

        print("T: {} K".format(T))
        EF, dQ = defects.find_EF(T, Tdef, NdsTdef, ECmin, ECmax, EVmin, EVmax,
                    callbackfunc = callbackfunc, 
                    eps_bisec  = eps_bisec,  nmaxiter_bisection = nmaxiter_bisection, initial = [EFmin_ini, EFmax_ini],
                    eps_newton = eps_newton, nmaxiter_newton = nmaxiter_newton, dump_newton = dump_newton,
                    h_newton = h_newton
                    )
        Z, ne, nh, Nds, NdsTdef, dGs, Qtot = defects.calculate_densities(T, Tdef, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax)
        if ne == 0.0:
            ne = 1.0
        if nh == 0.0:
            nh = 1.0
        for id in range(defects.ndefects):
            ynds[id][iT] = Nds[id]

        yEF.append(EF)
        yne.append(ne)
        ynh.append(nh)
        ydQ.append(Qtot)

        EFmin_ini = EF - 0.2
        EFmax_ini = EF + 0.2

    mprint("")
    mprint("T(defect frozen): {} K".format(Tdef))
    Twb = tkExcel(T_xlsx, 'w')
    mprint("{:8}: {:12} {:12} {:8} {:12} {:8}".format("T(K)", "EF(eV)", "ne(cm^-3)", "Ea(ne)(eV)", "nh(cm^-3)", "Ea(nh)(eV)"), end = '', fp = outfp)
    Twb.Print(["T(K)", "1000/T(K^-1)", "EF(eV)", "ne(cm^-3)", "Ea(ne)(eV)", "nh(cm^-3)", "Ea(nh)(eV)"], end = '')
    for id in range(defects.ndefects):
        d = defects.defects[id]
        mprint(" {:12}".format(d.name), end = '', fp = outfp)
        Twb.Print([d.name], end = '')
    mprint("", fp = outfp)
    Twb.Print(['dQ'])
    for iT in range(nT):
        T    = xT[iT]
        EF   = yEF[iT]
        ne   = yne[iT]
        nh   = ynh[iT]
        Qtot = ydQ[iT]
        if iT == 0:
            slopene = (log(yne[1]) - log(yne[0])) / (xT1000[1] - xT1000[0])
            slopenh = (log(ynh[1]) - log(ynh[0])) / (xT1000[1] - xT1000[0])
        elif iT == nT - 1:
            slopene = (log(yne[nT-1]) - log(yne[nT-2])) / (xT1000[nT-1] - xT1000[nT-2])
            slopenh = (log(ynh[nT-1]) - log(ynh[nT-2])) / (xT1000[nT-1] - xT1000[nT-2])
        else:
            slopene = (log(yne[iT+1]) - log(yne[iT-1])) / (xT1000[iT+1] - xT1000[iT-1])
            slopenh = (log(ynh[iT+1]) - log(ynh[iT-1])) / (xT1000[iT+1] - xT1000[iT-1])
        Eane = -slopene * kB / e * 1000.0
        Eanh = -slopenh * kB / e * 1000.0

        mprint("{:8.4g}: {:12.6g} {:12.6g} {:8.3f} {:12.6g} {:8.3f}".format(T, EF, ne, Eane, nh, Eanh), end = '', fp = outfp)
        Twb.Print([T, 1000.0 / T, EF, ne, Eane, nh, Eanh], end = '')
        for id in range(defects.ndefects):
            mprint(" {:12.6g}".format(Nds[id]), end = '', fp = outfp)
            Twb.Print([Nds[id]], end = '')
        mprint(" {:12.6g}".format(Qtot), fp = outfp)
        Twb.Print([Qtot])


    mprint("", fp = outfp)
    mprint("*** Save N-EF data to [{}]".format(T_xlsx))
    Twb.Close()

#=============================
# Plot graphs
#=============================
    fig = plt.figure(figsize = (12, 8))

    axdos = fig.add_subplot(2, 2, 1)
    axEF  = fig.add_subplot(2, 2, 2)
    axdH  = axEF.twiny()
    axNT  = fig.add_subplot(2, 2, 4)
    axNiT = fig.add_subplot(2, 2, 3)

    axdos.set_title("{}, for Point {}".format(dH_path, defects.plabels[iPoint]))

    axdos.plot(E_raw,  dos_raw, label = 'DOS(EF={:.3} eV)'.format(EF0), color = 'cyan',  linewidth =0.5)
    axdos.plot(E_norm, dos_raw, label = 'DOS(EV=0)',                    color = 'black', linewidth =1.0)
    yrange = [min(dos_raw), max(dos_raw)]
    axdos.plot([EV, EV],     yrange, label = '$E_V$',      color = 'blue', linestyle = 'dashed', linewidth = 0.5)
    axdos.plot([EC, EC],     yrange, label = '$E_C$',      color = 'blue', linestyle = 'dashed', linewidth = 0.5)
    EFeqmin = min(yEF)
    EFeqmax = max(yEF)
    axdos.plot([EFeqmin, EFeqmin], yrange, label = '$E_{F,eq}$(min)', color = 'red',  linestyle = 'dashed', linewidth = 1.0)
    axdos.plot([EFeqmax, EFeqmax], yrange, label = '$E_{F,eq}$(max)', color = 'blue', linestyle = 'dashed', linewidth = 1.0)
#    axdos.set_xlim([EFmin, EFmax])
    axdos.set_xlim([min([EFmin, view_Emin]), max([EFmax, view_Emax])])
    axdos.set_ylim([0.0, None])
    axdos.set_xlabel("$E$, $E - E_V$ (eV)")
    axdos.set_ylabel("DOS (states/cm$^3$)")
    _legend = axdos.legend()
    _legend.set_draggable(True)

    axEF.plot(xT, yEF, label = '$E_F$ (eV)')
    axEF.plot([Tmin, Tmax], [EV, EV], label = '$E_V$', color = 'red',  linestyle = 'dashed', linewidth = 1.0)
    axEF.plot([Tmin, Tmax], [EC, EC], label = '$E_C$', color = 'red',  linestyle = 'dashed', linewidth = 1.0)
    axEF.set_xlabel("$T$ (K)")
    axEF.set_ylabel("$E_F - E_V$ (eV)")
#    axEF.set_ylim([-1.0e19, 1.0e19])
    axEF.set_xlim([Tmin, Tmax])
#    axEF.legend()

    names = defects.get_names()
    for iname in range(len(names)):
        color = colors[iname % len(colors)]

        name = names[iname]
        data = defects.get_groupdata(name, iPoint)
        axdH.plot([], [], label = "{}".format(name), linestyle = 'dashed', linewidth = 0.5, color = color)

        minpts = minpts_list[iname]
        for ipt in range(len(minpts) - 1):
            axdH.plot([minpts[ipt][1], minpts[ipt+1][1]], [minpts[ipt][0], minpts[ipt+1][0]],
                    linestyle = 'dashed', linewidth = 0.5, color = color, 
                    marker = 'o', markersize = 3.0, markeredgewidth = 1, markerfacecolor = 'w', markeredgecolor = color)
    axdH.set_xlabel("$\Delta$$H$ (eV)")
#    axdH.set_xlim([view_dHmax, axdH.get_xlim()[0]])
    axdH.set_xlim([-0.2, axdH.get_xlim()[0]])

    h1, l1 = axEF.get_legend_handles_labels()
    h2, l2 = axdH.get_legend_handles_labels()
    _legend = axEF.legend(h1 + h2, l1 + l2, bbox_to_anchor=(1.05, 1.0), loc='upper left', borderaxespad = 0, fontsize = legend_fontsize)
    _legend.set_draggable(True)

    for i in [0, 1]:
        if i == 0:
            ax = axNiT
            xx = xT1000
            xxrange = [1000.0 / Tmax, 1000.0 / Tmin]
            ax.set_xlabel("1000/$T$ (K$^{-1}$)")
        else:
            ax = axNT
            xx = xT
            xxrange = [Tmin, Tmax]
            ax.set_xlabel("$T$ (K)")

        ax.plot(xx, yne, label = '$N_e$', color = 'red',  linestyle = '-', linewidth = 1.0)
        ax.plot(xx, ynh, label = '$N_h$', color = 'blue', linestyle = '-', linewidth = 1.0)
        for id in range(len(ynds)):
            d = defects.defects[id]
            ymax = max(ynds[id])
            if ymax < view_Nmin:
                continue
            ax.plot(xx, ynds[id], label = "{}({})".format(d.name, d.charge), linestyle = 'dashed', linewidth = 1.0)
        yrange = ax.get_ylim()
        yrange = [yrange[0], yrange[1] * 10.0]
        ax.plot(xxrange, [EV, EV], color = 'red', linestyle = 'dashed', linewidth = 0.5)
        ax.plot(xxrange, [EC, EC], color = 'red', linestyle = 'dashed', linewidth = 0.5)
        ax.set_ylabel("$N$  (cm$^{-3}$)")
        ax.set_yscale('log')
        ax.set_xlim(xxrange)
        ax.set_ylim([view_Nmin, 1.0e23])
#        ax.legend()
        _legend = ax.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left', borderaxespad = 0, fontsize = legend_fontsize)
        _legend.set_draggable(True)


# Rearange the graph axes so that they are not overlapped
    plt.tight_layout()

    """
    print("")
    print("Close graph window to terminate")
    plt.show()
    """

    plt.pause(0.1)

    app.terminate("", usage = usage, pause = True)

    if outfp:
        outfp.close()

def exec_EF():
    """
    概要:
        フェルミ準位依存性の計算とプロットを実行します。
    詳細説明:
        欠陥形成エネルギーのExcelファイルとVASP計算結果ファイルを読み込み、
        様々なフェルミ準位における電子・正孔密度、欠陥密度、および全電荷の中性条件からの偏差を計算します。
        計算結果はExcelファイルに保存され、グラフとして可視化されます。
    """
    global mode, plot_mode, T0, EF0
    global dH_path, CAR_path, CAR_path_defect, DOSCAR_path, Vcell
    global E_raw, E_norm, dos_raw, nDOS, Enorm_min, Enorm_max, Enorm_step, Emin, Emax, Estep
    global EV, EC
    global EFmin, EFmax, nEF, EFstep
    global Einteg0
    global outfp # グローバル変数としてoutfpを宣言

    kBTe    = kB * T0 / e
    kBTedef = kB * Tdef / e
    if T0 < Tdef:
        dE = nrange * kBTedef
    else:
        dE = nrange * kBTe

    vasp = tkVASP()
    CAR_path = vasp.getdir(CAR_path)
    INCAR_path   = vasp.get_INCAR(CAR_path)
    POSCAR_path  = vasp.get_POSCAR(CAR_path)
    CONTCAR_path = vasp.get_CONTCAR(CAR_path)
    OUTCAR_path  = vasp.get_OUTCAR(CAR_path)
    EIGENVAL_path = vasp.get_VASPPath(CAR_path, 'EIGENVAL')
    DOSCAR_path   = vasp.get_VASPPath(CAR_path, 'DOSCAR')

    CAR_path_defect = vasp.getdir(CAR_path_defect)
    POSCAR_path_defect = vasp.get_POSCAR(CAR_path_defect)
    
    print("")
    print("*** Read crystal structure from [{}]".format(POSCAR_path))
    cry = vasp.read_poscar(POSCAR_path)
#   cry.PrintInf("cell")
    Vcell = cry.Volume()
    mprint("  volume: {:12.6f} A^-3".format(Vcell))

    mprint("Crystal structure of defect model from [{}]".format(POSCAR_path_defect))
    cry_d = vasp.read_poscar(POSCAR_path_defect)
#    a_d, b_d, c_d, alpha_d, beta_d, gamm_d = cry_d.LatticeParameters()
    Vcell_d = cry_d.Volume()
#    cry.PrintInf("cell")
#    mprint("  cell: {:12.8f} {:12.8f} {:12.8f} A   {:10.6f} {:10.6f} {:10.6f}".format(a_d, b_d, c_d, alpha_d, beta_d, gamm_d), fp = outfp)
    mprint("  volume: {:12.6f} A^-3".format(Vcell_d))

    print("")
    print("Read defect formation enthalpies from [{}]".format(dH_path))
    defects = tkDefects()
#    defects.SetVolume(Vcell)
    defects.read_excel(dH_path, Vcell * int(Vcell_d / Vcell + 0.2))
    defects.Print()

    outputfile    = modify_path(dH_path, '-out-{}.txt'.format(defects.plabels[iPoint]))
    NEF_xlsx      = modify_path(dH_path, '-N-EF-{}.xlsx'.format(defects.plabels[iPoint]))
    dHEF_all_xlsx = modify_path(dH_path, '-dH-EF-all-{}.xlsx'.format(defects.plabels[iPoint]))
    dHEF_min_xlsx = modify_path(dH_path, '-dH-EF-min-{}.xlsx'.format(defects.plabels[iPoint]))
    print("  Remove [{}]".format(outputfile))
    delete_file(outputfile)
    print("  Remove [{}]".format(NEF_xlsx))
    delete_file(NEF_xlsx)
    print("  Remove [{}]".format(dHEF_all_xlsx))
    delete_file(dHEF_all_xlsx)
    print("  Remove [{}]".format(dHEF_min_xlsx))
    delete_file(dHEF_min_xlsx)

    print("")
    print("*** Open output file {}".format(outputfile))
    outfp = open(outputfile, 'w')
    if outfp is None:
        app.terminate("Error: Can not write to [{}]".format(outputfile), pause = True)

    if plot_mode == 'all':
        plotall = True
    else:
        plotall = False

    mprint("", fp = outfp)
    mprint("=======================================================", fp = outfp)
    mprint(" Calculate EF dependence of semiconductor properties", fp = outfp)
    mprint("=======================================================", fp = outfp)
    mprint("mode     : ", mode, fp = outfp)
    mprint("plot mode: ", plot_mode, fp = outfp)
    mprint("plot all: ", plotall, fp = outfp)

    mprint("", fp = outfp)
    mprint("Files:", fp = outfp)
    mprint("  dH input  : ", dH_path, fp = outfp)
    mprint("", fp = outfp)
    mprint("Files:", fp = outfp)
    mprint("  dH input  : ", dH_path, fp = outfp)
    mprint("  Output    : ", outputfile, fp = outfp)
    mprint("  N-EF      : ", NEF_xlsx, fp = outfp)
    mprint("  dH-EF(all): ", dHEF_all_xlsx, fp = outfp)
    mprint("  dH-EF(min): ", dHEF_min_xlsx, fp = outfp)

    mprint("", fp = outfp)
    mprint("Defects red from [{}]".format(dH_path), fp = outfp)
    defects.Print(outfp = outfp) 
    mprint("", fp = outfp)
    mprint("To be analyzed for Point {}".format(defects.plabels[iPoint]), fp = outfp)
    mprint("T(defects): {} K".format(Tdef), fp = outfp)
    
    mprint("", fp = outfp)
    mprint("VASP files  :", fp = outfp)
    mprint("  CAR dir   : ", CAR_path, fp = outfp)
    mprint("  INCAR     : ", INCAR_path, fp = outfp)
    mprint("  POSCAR    : ", POSCAR_path, fp = outfp)
    mprint("  CONTCAR   : ", CONTCAR_path, fp = outfp)
    mprint("  OUTCAR    : ", OUTCAR_path, fp = outfp)
    mprint("  EIGENVAL  : ", EIGENVAL_path, fp = outfp)
    mprint("  DOSCAR    : ", DOSCAR_path, fp = outfp)
    mprint("  POSCAR(defect): ", POSCAR_path_defect, fp = outfp)

    outcarinf, eigenvalinf, bandedgeinf, dosinf, EV, EC, Eg, EHOMO, ELUMO = \
        read_files(vasp, cry, cry_d, defects, POSCAR_path, OUTCAR_path, EIGENVAL_path, DOSCAR_path, fp = outfp)

# EF plot range
    if EFmin is None: EFmin = EV + dEFmin
    if EFmax is None: EFmax = EC + dEFmax
    EFstep = (EFmax - EFmin) / (nEF - 1)
    mprint("", fp = outfp)
    mprint("EF dependence", fp = outfp)
    mprint("  EF range: {} - {}, {} eV step".format(EFmin, EFmax, EFstep), fp = outfp)

    mprint("", fp = outfp)
    mprint("Integration configuration", fp = outfp)
    mprint("  Integration E range: {} - {} eV, or EF+-{}*kBT eV".format(Einteg0, Einteg1, nrange), fp = outfp)

    ECmin = EC - Estep
    ECmax = EC + dE
    EVmin = EV - dE
    EVmax = EV + Estep

#==============================================
# Print transition levels
#==============================================
    name, allpts_list, minpts_list = print_transition_level(defects)
#    exit()
    
#==============================================
# Start calculations
#==============================================
# At Tdef
#    EFmin_ini = min([EV - dEFmin, EF - dEFmin])
#    EFmax_ini = max([EC + dEFmax, EF + dEFmax])
    EFmin_ini = EV + dEFmin
    EFmax_ini = EC + dEFmax
    mprint("", fp = outfp)
    mprint("Calculate defect densities at T(electron) = T(defects) = {} K".format(Tdef))
    mprint("  Calculate EF,eq: Bisecton method: eps = {}, nmaxiter = {}".format(eps_bisec, nmaxiter_bisection), fp = outfp)
    mprint("                        Initial EF range: {:12.4f} - {:12.4f} eV".format(EFmin_ini, EFmax_ini))
    mprint("                   Newton method  : eps = {}, nmaxiter = {}, dump = {}".format(eps_newton, nmaxiter_newton, dump_newton), fp = outfp)

    if EFdef == 'eq':
        EFeqTdef, dQ = defects.find_EF(Tdef, Tdef, None, ECmin, ECmax, EVmin, EVmax,
                    callbackfunc = callbackfunc, 
                    eps_bisec  = eps_bisec,  nmaxiter_bisection = nmaxiter_bisection, initial = [EFmin_ini, EFmax_ini],
                    eps_newton = eps_newton, nmaxiter_newton = nmaxiter_newton, dump_newton = dump_newton,
                    h_newton = h_newton
                    )
        if EFeqTdef is None:
            app.terminate("Error in tkDefects.find_EF: Could not reach convergence for EF", pause = True)
    else:
        try:
            EFeqTdef = float(EFdef)
        except:
            app.terminate("Error in tkDefects.find_EF: EFdef [{}] must be numeral".format(EFdef), pause = True)

    mprint("", fp = outfp)
    mprint("  EF,eq({:8.3f} K)={:12.6g} eV".format(Tdef, EFeqTdef), fp = outfp)

    Z, ne, nh, Nds, NdsTdef, dGs, Qtot = defects.calculate_densities(T0, Tdef, None, EFeqTdef, ECmin, ECmax, EVmin, EVmax)

    for id in range(defects.ndefects):
        d = defects.defects[id]
        label  = "{}_{}".format(d.atom, d.site)
        labelq = "{}^{}".format(label, d.charge)
        mprint("    {:>16} (N0={:5g} Ndoped={:5g}):".format(labelq, d.N0, d.Ndoped), end = '', fp = outfp)
        mprint("       Nds={:12.4g} cm-3  dGs={:12.4g} eV".format(Nds[id], dGs[id]), fp = outfp)
    mprint("", fp = outfp)

    mprint("  Total defect densities at {} K".format(Tdef), fp = outfp)
    for defect_name in NdsTdef.keys():
        Nsite = NdsTdef[defect_name] * defects.V * 1.0e-24
        mprint(f"    {defect_name:>8}: {NdsTdef[defect_name]:12.6g} cm-3 ({Nsite:8.2f} sites)", fp = outfp)

# At T0, EF,eq(T0)
    EFmin_ini = EV + dEFmin
    EFmax_ini = EC + dEFmax
    mprint("", fp = outfp)
    if T0 <= Tdef:
        mprint("Calculate electron distribution at T(electron) = {} K".format(T0))
        mprint("  Total defect densities are frozen at {} K with EF = {:12.4f} eV".format(Tdef, EFeqTdef), fp = outfp)
        for key in NdsTdef.keys():
            Nsite = NdsTdef[key] * defects.V * 1.0e-24
            mprint("    {:>8}: {:12.6g} cm-3 ({:8.2f} sites)".format(key, NdsTdef[key], Nsite), fp = outfp)
    else:
        mprint("T(electron) is higher than T(defects). Defect densities are calculated at T(electron)")
        mprint("Calculate electron distribution at T(electron) = {} K, T(defects) = {} K".format(T0, T0))

    mprint("  Calculate EF,eq: Bisecton method: eps = {}, nmaxiter = {}".format(eps_bisec, nmaxiter_bisection), fp = outfp)
    mprint("                        Initial EF range: {:12.4f} - {:12.4f} eV".format(EFmin_ini, EFmax_ini))
    mprint("                   Newton method  : eps = {}, nmaxiter = {}, dump = {}".format(eps_newton, nmaxiter_newton, dump_newton), fp = outfp)

    EFeqT0, dQ = defects.find_EF(T0, Tdef, NdsTdef, ECmin, ECmax, EVmin, EVmax,
                    callbackfunc = callbackfunc, 
                    eps_bisec  = eps_bisec,  nmaxiter_bisection = nmaxiter_bisection, initial = [EFmin_ini, EFmax_ini],
                    eps_newton = eps_newton, nmaxiter_newton = nmaxiter_newton, dump_newton = dump_newton,
                    h_newton = h_newton
                    )
    if EFeqT0 is None:
        app.terminate("Error in tkDefects.find_EF: Could not reach convergence for EF", pause = True)

    mprint("", fp = outfp)
    mprint("  EF,eq({:8.3g} K)={:12.6g} eV".format(T0, EFeqT0), fp = outfp)
    Z, ne, nh, Nds, NdsTdef, dGs, Qtot = defects.calculate_densities(T0, Tdef, NdsTdef, EFeqT0, ECmin, ECmax, EVmin, EVmax)
    for id in range(defects.ndefects):
        d = defects.defects[id]
        label  = "{}_{}".format(d.atom, d.site)
        labelq = "{}^{}".format(label, d.charge)
        mprint(f"    {labelq:>16} (Nds(tot)={NdsTdef[labelq]:12.6g} cm-3):", end = '', fp = outfp)
#        mprint(f"    {labelq:>16} (Nds(tot)={NdsTdef[label]:12.6g} cm-3):", end = '', fp = outfp)
        mprint(f"       Nds={Nds[id]:12.6g} cm-3  dGs={dGs[id]:12.4g} eV", fp = outfp)
    
# At T0, various EF
    outwb = tkExcel(NEF_xlsx, 'w')

    xEF = [EFmin + i * EFstep for i in range(nEF)]
    mprint("", fp = outfp)
    mprint("Calculate defect and carrier densities at {} K at {} point as functions of EF:"
            .format(T0, defects.plabels[iPoint]), fp = outfp)
    if T0 <= Tdef:
        mprint("  T(electron) is lower than T(defects).", fp = outfp)
        mprint("  Total defect densities frozen at {} K with EF = {:12.4f} eV".format(Tdef, EFeqTdef), fp = outfp)
        for key in NdsTdef.keys():
            Nsite = NdsTdef[key] * defects.V * 1.0e-24
            mprint("    {:>8}: {:12.6g} cm-3 ({:8.2f} sites)".format(key, NdsTdef[key], Nsite), fp = outfp)
    else:
        mprint("T(electron) is higher than T(defects). Defects are not frozen")
        mprint("Calculate defect densities at T(electron) = {} K, T(defects) = {} K".format(T0, T0))

    mprint("{:8}: {:12} {:12}".format("EF(eV)", "ne(cm^-3)", "nh(cm^-3)"), end = '', fp = outfp)
    outwb.Print(["EF(eV)", "ne(cm^-3)", "nh(cm^-3)"], end = '')
    for id in range(defects.ndefects):
        d = defects.defects[id]
        label = "{}({})".format(d.name, d.charge)
        mprint(" {:12}".format(label), end = '', fp = outfp)
        outwb.Print(["{}".format(label)], end = '')
    mprint(" {:12}".format('dQ'), fp = outfp)
    outwb.Print(['dQ'])

    ydQ    = []
    ylogdQ = []
    xEF_dQp  = []
    xEF_dQm  = []
    ylogdQp = []
    ylogdQm = []
    yne  = []
    ynh  = []
    ynds = []
    for id in range(defects.ndefects):
        ynds.append([0.0]*nEF)
    
    for i in range(nEF):
        EF = xEF[i]
        Z, ne, nh, Nds, NdsTdef, dGs, Qtot = defects.calculate_densities(T0, Tdef, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax)
        for id in range(defects.ndefects):
            ynds[id][i] = Nds[id]

        ydQ.append(Qtot)
        yne.append(ne)
        ynh.append(nh)

# Make an array for plotting log(|dQ|) - EF
        if Qtot > 0.0:
            logdQ = log(Qtot) / log10
            ylogdQ.append(logdQ)
            xEF_dQp.append(EF)
            ylogdQp.append(logdQ)
        elif Qtot < 0.0:
            logdQ = -log(-Qtot) / log10
            ylogdQ.append(logdQ)
            xEF_dQm.append(EF)
            ylogdQm.append(logdQ)
        else:
            ylogdQ.append(0.0)
            xEF_dQp.append(0.0)
            ylogdQp.append(0.0)

        mprint(f"{EF:8.4g}: {ne:12.4g} {nh:12.4g}", end = '', fp = outfp)
        outwb.Print([EF, ne, nh], end = '')
        for id in range(defects.ndefects):
            mprint(" {:12.4g}".format(Nds[id]), end = '', fp = outfp)
            outwb.Print([Nds[id]], end = '')
        mprint(" {:12.4g}".format(Qtot), fp = outfp)
        outwb.Print([Qtot])

    outwb.Close()

#    mprint("", fp = outfp)
#    mprint("*** Open [{}] to save N-EF data".format(NEF_xlsx))
#    for i in range(len(ylogdQ)):
#        print(f"{xEF[i]:10.4g}  {ydQ[i]:10.4g}  {ylogdQ[i]:10.4g}")

    mprint("", fp = outfp)
    mprint("*** Open [{}] to save dH-EF data (all)".format(dHEF_all_xlsx))
    mprint("*** Open [{}] to save dH-EF data (min)".format(dHEF_min_xlsx))
    allwb = tkExcel(dHEF_all_xlsx, 'w')
    minwb = tkExcel(dHEF_min_xlsx, 'w')
    for id in range(defects.ndefects):
        d = defects.defects[id]
        name = d.name
        q    = d.charge
        allwb.Print([name, q])
        allwb.Print(["EF(eV)", "dH(eV)"])
        allwb.Print([EFmin, d.dH0s[iPoint] + q * EFmin])
        allwb.Print([EFmax, d.dH0s[iPoint] + q * EFmax])
    allwb.Close()
    minwb.Close()

#=============================
# Plot graphs
#=============================
    fig = plt.figure(figsize = (12, 8))
    plot_event = tkPlotEvent(plt)

    axdos = fig.add_subplot(2, 2, 1)
    axdQ  = fig.add_subplot(2, 2, 3)
    axdH  = fig.add_subplot(2, 2, 2)
    axdN  = fig.add_subplot(2, 2, 4)

# plot_event初期化
    root = get_window_from_plt(plt)
    plot_event = tkPlotEvent(plt)
    plot_event.prepare_annotation()
    plot_event.prepare_move_text(fig)
    plot_event.prepare_popup_menu(fig, parent = root)
    popup_menu = plot_event.popup_menu.menu

#    selector0 = RangeSelector('', ax1, color = 'green', print_level = 0)
#    selector1 = RangeSelector('', ax2, color = 'blue', print_level = 0)

    vars = tkParams() #cparams.copy()
    vars.caller = "EF"
    vars.plot_event = plot_event
#    vars.selector0 = selector0
#    vars.selector1 = selector1

    if len(app.config.plugin_dir) >= 1 \
            and (app.config.plugin_dir[0] != "/" and app.config.plugin_dir[0] != "\\" and app.config.plugin_dir[2] != "\\"):
        plugin_dir = app.replace_path(None, template = ["{dirname}", app.config.plugin_dir])
    else:
        plugin_dir = app.config.plugin_dir

    print()
    print(f"Read plugins from : {plugin_dir}")
    module_names, modules = app.load_modules(plugin_dir, "*.py", target = "popup_menu", sort = True, is_print = True)
    for m in modules:
        if hasattr(m, "add_popup_menu"):
            m.add_popup_menu(popup_menu, app = app, vars = vars, parent = root)

    axdos.set_title("{}, for Point {}".format(dH_path, defects.plabels[iPoint]))
    axdos.plot(E_raw,  dos_raw, label = 'DOS(EF={:.3} eV)'.format(EF0), picker = True, color = 'cyan', linewidth =0.3)
    axdos.plot(E_norm, dos_raw, label = 'DOS(EV=0)',                    picker = True, color = 'black', linewidth =1.0)
    yrange = [min(dos_raw), max(dos_raw)]
    axdos.plot([EV, EV],     yrange, label = '$E_V$',      color = 'blue', linestyle = 'dashed', linewidth = 0.5)
    axdos.plot([EC, EC],     yrange, label = '$E_C$',      color = 'blue', linestyle = 'dashed', linewidth = 0.5)
    axdos.plot([EFeqT0,   EFeqT0],   yrange, label = '$E_{F,eq}$(T0)',   color = 'red',  linestyle = 'dashed', linewidth = 1.0)
    axdos.plot([EFeqTdef, EFeqTdef], yrange, label = '$E_{F,eq}$(Tdef)', color = 'green',  linestyle = 'dashed', linewidth = 1.0)
#    axdos.set_xlim([EFmin, EFmax])
    xmin = min([EFmin, view_Emin])
    xmax = max([EFmax, view_Emax])
    axdos.set_xlim([xmin, xmax])
    ymin, ymax = minmax_xy(x = E_raw, y = dos_raw, xmin = xmin, xmax = xmax)
    axdos.set_ylim([0.0, ymax])
    axdos.set_xlabel("$E$, $E - E_V$ (eV)")
    axdos.set_ylabel("DOS (states/cm$^3$)")
    _legend = axdos.legend()
    _legend.set_draggable(True)

#    axdQ.plot(xEF, ylogdQ, label = 'log$_{10}$|$\Delta$$Q$|', picker = True, marker = 'o', markersize = 2.0)
    axdQ.plot(xEF_dQp, ylogdQp, label = 'log$_{10}$|$\Delta$$Q$| (dQ >= 0)', picker = True, marker = 'o', markersize = 2.0, markeredgecolor = 'blue', markerfacecolor = 'blue')
    axdQ.plot(xEF_dQm, ylogdQm, label = 'log$_{10}$|$\Delta$$Q$| (dQ < 0)',  picker = True, marker = 'o', markersize = 2.0, markeredgecolor = 'red', markerfacecolor = 'red')
    axdQ.plot(axdQ.get_xlim(), [0.0, 0.0], color = 'red', linestyle = '-', linewidth = 1.0)
    yrange = axdQ.get_ylim()
    axdQ.plot([EV, EV],     yrange, label = '$E_V$',      color = 'blue', linestyle = 'dashed', linewidth = 0.5)
    axdQ.plot([EC, EC],     yrange, label = '$E_C$',      color = 'blue', linestyle = 'dashed', linewidth = 0.5)
    axdQ.plot([EFeqT0,   EFeqT0],   yrange, label = '$E_{F,eq}$(T0)',   color = 'red',  linestyle = 'dashed', linewidth = 1.0)
    axdQ.plot([EFeqTdef, EFeqTdef], yrange, label = '$E_{F,eq}$(Tdef)', color = 'green',  linestyle = 'dashed', linewidth = 1.0)
    axdQ.set_xlabel("$E_F - E_V$ (eV)")
    axdQ.set_ylabel("log$_{10}$ |$\Delta$$Q$ / $e$|")
    axdQ.set_xlim([EFmin, EFmax])
#    axdQ.set_ylim([-1.0e19, 1.0e19])
    _legend = axdQ.legend()
    _legend.set_draggable(True)
    
    alpha = 0.8
    names = defects.get_names()
    ymin = 1.0e300
    ymax = -1.0e300
    for iname, name in enumerate(names):
        groupdata = defects.get_groupdata(name, iPoint)
        allpts, minpts = defects.find_all_transitions(name, groupdata, EFmin, EFmax)
        for pts in minpts:
            ymin = min([ymin, pts[0], pts[1]])
            ymax = max([ymax, pts[0], pts[1]])

    yrange = [ymin, min(ymax, view_dHmax)]
    if dHmin is not None: yrange[0] = dHmin
    if dHmax is not None: yrange[1] = dHmax
    def plot_dH_EF(axdH):
        axdH.set_title('$E_{F,eq}$: ' + "{:8.3g} eV(Tdef={} K)".format(EFeqTdef, Tdef)
                                      + "  {:8.3g} eV(T0={} K)".format(EFeqT0, T0))

        for iname, name in enumerate(names):
            color = colors[iname % len(colors)]
            args    = {'color': color, 'markersize': 3.0, 'markeredgewidth': 0, 'markerfacecolor': color }

            groupdata = defects.get_groupdata(name, iPoint)
            ngroupdata = len(groupdata)

            minwb.Print([name])
            minwb.Print(["EF(eV)", "dH(eV)", "q(-)", "q(+)"])

            allpts, minpts = defects.find_all_transitions(name, groupdata, EFmin, EFmax)
        
            min_x = []
            min_y = []
            for ipt in range(len(minpts)):
                min_x.append(minpts[ipt][0])
                min_y.append(minpts[ipt][1])

            label = name
            largs = { 'label': label }
            line, = axdH.plot(min_x, min_y, linestyle = '-', marker = 'o', **args, **largs)
#        axdH.plot(min_x, min_y, picker = True, linestyle = '-', marker = 'o', **args, **largs)

            format = "{label}: line#{iline} data#{idata}:\n EF - EV={x_list:8.3g} eV  dH={y_list:8.3g} eV"
            plot_event.annotation.add_line(label, axdH, axdH, min_x, min_y, line,
                    inf_list = {"x_label": "EF - EV (eV)", "y_label": "dH (eV)", "x_list": min_x, "y_list": min_y,},
                    annotation_format = format, inf_format = format)

            plot_event.move_text.add_annotation(axdH, axdH, x_list = min_x, y_list = min_y,
                     frac = None, ylim = yrange,
                     text = name, fontsize = 10, ha = 'right', va = 'center', color = color, alpha = alpha, fc = "w", ec = "none")

            for id in range(0, ngroupdata):
                d = groupdata[id]
                name = d["name"]
                dH0 = d["dH0"]
                q   = d["charge"]
                if plotall:
                    label = f"{name} (dashed line)"
                    line, = axdH.plot([EFmin, EFmax], [dH0 + q * EFmin, dH0 + q * EFmax], linestyle = 'dashed', linewidth = 0.5, **args)
                    plot_event.annotation.add_line(label, axdH, axdH, min_x, min_y, line,
                            inf_list = {"x_label": "EF - EV (eV)", "y_label": "dH (eV)", "x_list": min_x, "y_list": min_y,},
                            annotation_format = format, inf_format = format)

        axdH.plot([EV, EV],     yrange, color = 'blue', linestyle = 'dashed', linewidth = 0.5)
        axdH.plot([EC, EC],     yrange, color = 'blue', linestyle = 'dashed', linewidth = 0.5)
        axdH.plot([EFeqT0,   EFeqT0],   yrange, label = '$E_{F,eq}$(T0)',   color = 'red',  linestyle = 'dashed', linewidth = 1.0)
        axdH.plot([EFeqTdef, EFeqTdef], yrange, label = '$E_{F,eq}$(Tdef)', color = 'green',  linestyle = 'dashed', linewidth = 1.0)
        axdH.set_xlim([EFmin, EFmax])
        axdH.set_ylim(yrange)
        axdH.set_xlabel("$E_F - E_V$ (eV)")
        axdH.set_ylabel("$\Delta$$H$ (eV)")
        _legend = axdH.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left', borderaxespad = 0, fontsize = legend_fontsize)
        _legend.set_draggable(True)

    plot_dH_EF(axdH)


    format = "{label}: line#{iline} data#{idata}:\n EF - EV={x_list:8.3g} eV  {x_label}={y_list:8.3g} 1/cm3"
    label = '$N_e$'
    line, = axdN.plot(xEF, yne, label = label, picker = True, color = 'red',  linestyle = '-', linewidth = 1.0)
    plot_event.annotation.add_line(label, axdN, axdN, xEF, yne, line,
                inf_list = {"x_label": "EF - EV (eV)", "y_label": label, "x_list": xEF, "y_list": yne},
                annotation_format = format, inf_format = format)
    plot_event.move_text.add_annotation(axdN, axdN, x_list = xEF, y_list = yne,
                frac = None, xlim = [EFmin, EFmax], ylim = [view_Nmin, 1.0e23],
                text = label, fontsize = 10, ha = 'right', va = 'center', color = 'red', alpha = alpha, fc = "w", ec = "none")
    label = '$N_h$'
    line, axdN.plot(xEF, ynh, label = label, picker = True, color = 'blue', linestyle = '-', linewidth = 1.0)
    plot_event.annotation.add_line(label, axdN, axdN, xEF, ynh, line,
                inf_list = {"x_label": "EF - EV (eV)", "y_label": label, "x_list": xEF, "y_list": ynh},
                annotation_format = format, inf_format = format)
    plot_event.move_text.add_annotation(axdN, axdN, x_list = xEF, y_list = ynh,
                frac = None, xlim = [EFmin, EFmax], ylim = [view_Nmin, 1.0e23],
                text = label, fontsize = 10, ha = 'right', va = 'center', color = 'blue', alpha = alpha, fc = "w", ec = "none")
    ncolor = len(colors)
    for id in range(len(ynds)):
        d = defects.defects[id]
        ymax = max(ynds[id])
        if ymax < view_Nmin:
            continue

        color = colors[id % ncolor]
        label = f"{d.name}({d.charge})"
        line, axdN.plot(xEF, ynds[id], label = label, picker = True, linestyle = 'dashed', linewidth = 1.0, color = color)

        plot_event.annotation.add_line(label, axdN, axdN, xEF, ynds[id], line,
                inf_list = {"x_label": "EF - EV (eV)", "y_label": label, "x_list": xEF, "y_list": ynds[id]},
                annotation_format = format, inf_format = format)

        plot_event.move_text.add_annotation(axdN, axdN, x_list = xEF, y_list = ynds[id], 
                frac = None, xlim = [EFmin, EFmax], ylim = [view_Nmin, 1.0e23],
                text = label, fontsize = 10, ha = 'right', va = 'center', color = color, alpha = alpha, fc = "w", ec = "none")

    yrange = axdN.get_ylim()
    yrange = [yrange[0], yrange[1] * 10.0]
    axdN.plot([EV, EV],     yrange, color = 'blue', linestyle = 'dashed', linewidth = 0.5)
    axdN.plot([EC, EC],     yrange, color = 'blue', linestyle = 'dashed', linewidth = 0.5)
    axdN.plot([EFeqT0,   EFeqT0],   yrange, label = '$E_{F,eq}$(T0)',   color = 'red',  linestyle = 'dashed', linewidth = 1.0)
    axdN.plot([EFeqTdef, EFeqTdef], yrange, label = '$E_{F,eq}$(Tdef)', color = 'green',  linestyle = 'dashed', linewidth = 1.0)
    axdN.set_xlabel("$E_F - E_V$ (eV)")
    axdN.set_ylabel("$N$  (cm$^{-3}$)")
    axdN.set_yscale('log')
    axdN.set_xlim([EFmin, EFmax])
    axdN.set_ylim([view_Nmin, 1.0e23])
#    axdN.legend()
    _legend = axdN.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left', borderaxespad = 0, fontsize = legend_fontsize)
    _legend.set_draggable(True)

# Rearange the graph axes so that they are not overlapped
    plt.tight_layout()

    def rescale(ax: object, EFmin: float, EFmax: float, dHmin: float, dHmax: float):
        """
        概要:
            グラフの軸のスケールを再設定します。
        引数:
            :param ax: スケールを再設定するAxesオブジェクト。
            :type ax: object
            :param EFmin: EF軸の最小値 (eV)。
            :type EFmin: float
            :param EFmax: EF軸の最大値 (eV)。
            :type EFmax: float
            :param dHmin: dH軸の最小値 (eV)。
            :type dHmin: float
            :param dHmax: dH軸の最大値 (eV)。
            :type dHmax: float
        """
#        ax.cla()
#        plot_dH_EF(ax)
        ax.set_xlim([EFmin, EFmax])
        ax.set_ylim([dHmin, dHmax])
#        plt.draw()
        plt.pause(1.0e-4)
    
    vars.scale_callback = rescale
    vars.dH_axis = axdH
    EFlim = axdH.get_xlim()
    dHlim = axdH.get_ylim()
    vars.EF_min = EFlim[0]
    vars.EF_max = EFlim[1]
    vars.dH_min = dHlim[0]
    vars.dH_max = dHlim[1]
    plot_event.register_annotation_event(fig, activate = False, print_level = 0)
    plot_event.register_move_text_event(fig, activate = False)
    plot_event.register_popup_menu_event()
#    plot_event.register_pick(fig) # callback = lambda event: plot_event.onclick(event))


    plt.pause(1.0e-4)

    app.terminate("", pause = True) #usage = usage, pause = True)

    if outfp:
        outfp.close()

def main():
    """
    概要:
        スクリプトのメイン実行関数。
    詳細説明:
        アプリケーションの設定を読み込み、コマンドライン引数を解析してグローバル変数を更新します。
        ログファイルを初期化し、指定されたモードに応じてフェルミ準位依存性または
        温度依存性の計算とプロットを実行します。不正なモードが指定された場合はエラーで終了します。
    """
    app.config.plugin_dir = app.replace_path(None, template = ["{dirname}", "plugin/vasp_defect"])
    appconfig_path, config = app.read_app_config(print_level = 1)

    updatevars()

    vasp = tkVASP()
    base_path = vasp.getdir(CAR_path)
    logfile = app.replace_path(dH_path, template = ["{dirname}", "{filebody}-defect-out.txt"])
#    logfile = os.path.join(base_path, 'vasp_defect-out.txt')
    print("")
    print(f"Open logfile [{logfile}]")
    app.redirect(targets = ["stdout", logfile], mode = 'w')

    if mode == 'EF':
        exec_EF()
    elif mode == 'T':
        exec_T()
    else:
        app.terminate("Error: Invalid mode [{}]".format(mode), usage = usage, pause = True)


if __name__ == '__main__':
    main()