"""
概要:
    VASP計算結果から欠陥モデルのバンドフィリング補正を推定します。

詳細説明:
    本スクリプトは、欠陥モデルの形成エネルギー計算において必要となるバンドフィリング補正 dE_BF を算出します。
    理想結晶と欠陥結晶のVASP計算結果 (OUTCAR, DOSCAR, EIGENVAL) を比較し、
    dEVBMによるバンドアラインメント後の状態占有率から電子と正孔によるエネルギー補正量を評価します。
    また、計算結果のDOSとバンド構造をグラフとして可視化し、サマリーファイルに結果を保存します。

関連リンク:
    vasp_correction_bandfilling_usage
"""
import os
import sys
import shutil
import glob
import csv
import re
import numpy as np
from numpy import exp, log, sin, cos, tan, arcsin, arccos, arctan, pi
from scipy.interpolate import interp1d
from pprint import pprint
from matplotlib import pyplot as plt


from tklib.tkfile import tkFile
from tklib.tkinifile import tkIniFile
import tklib.tkre as tkre
from tklib.tkutils import save_csv
from tklib.tkapplication import tkApplication
from tklib.tkutils import IsDir, IsFile, SplitFilePath
from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
from tklib.tksci.tksci import Reduce01, Round
from tklib.tksci.tkconvolution import convolution, convolve_func
from tklib.tksci.tkconvolution import convolve_xydata
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.tkatomtype import tkAtomType
from tklib.tkcrystal.tkvasp import tkVASP
from tklib.tkcrystal.tkbandstructure import plot_band
from tklib.tksci import tkequation
import tklib.tkcsv


#================================
# global parameters
#================================
debug = 0

validmodes = ['BF']
mode = 'BF'

CAR_dir1 = '.'
CAR_dir2 = '.'
summaryfile = None

# dEVBM obtained by VBM correction
dEVBM = ''

# FWHM of Gaussian smearing function for DOS
WG_DOS = 0.1 # eV


# Energy to search edge energies
EF0  = 0.1
# occupancy threthold to separate HOMO and LUMO
occ_th = 0.5

# Critera DOS value to search band edges
DOSth = 1.0e-5

# Max number of electrons occupies one state
Nemax = 1.0

# Plot configuration
band_marker_size       = 4
band_marker_edge_width = 0.5

Emin = -10.0  # eV
Emax =  10.0  # eV

#colors = ['#000000', '#ff0000', '#00aa00', '#0000ff', '#aaaa00', '#ff00ff', '#00ffff', '#aa0000', '#00aa00', '#0000aa']
#colors = ['k', 'r', 'g', 'b', 'y', 'm', 'c']
figsize = (10, 8)
fontsize        = 16
titlefontsize   = 12
labelfontsize   = 12
legend_fontsize = 12


app = tkApplication()


#=============================
# Treat argments
#=============================
def usage():
    """
    概要:
        スクリプトの正しい使用方法を標準出力に表示します。

    詳細説明:
        スクリプトの実行に必要なコマンドライン引数のフォーマットと、
        各引数の意味について説明します。

    戻り値:
        :returns: なし
        :rtype: None
    """
    global mode
    
    if mode not in validmodes:
        mode = validmodes[0]

    print("")
    print("Usage:")
    print("  (a)  python {} mode CAR_dir(ideal) CAR_dir(defect) dEVBM WG Emin Emax EF0".format(sys.argv[0]))
    print("         mode: {}".format(validmodes))
    print("     ex: python {} {} {} {} {} {} {} {} {}".format(sys.argv[0], mode, CAR_dir1, CAR_dir2, dEVBM, WG_DOS, Emin, Emax, EF0))

def updatevars():
    """
    概要:
        コマンドライン引数に基づいてグローバル変数を更新します。

    詳細説明:
        sys.argvから引数を読み込み、mode、CAR_dir1、CAR_dir2、dEVBM、WG_DOS、Emin、Emax、EF0
        などのグローバルパラメータに値を設定します。不正なモードが指定された場合はエラーで終了します。

    戻り値:
        :returns: なし
        :rtype: None
    """
    global mode
    global CAR_dir1, CAR_dir2
    global EF0, dEVBM
    global Emin, Emax, WG_DOS

    mode     = getarg     (1, mode)
    CAR_dir1 = getarg     (2, CAR_dir1)
    CAR_dir2 = getarg     (3, CAR_dir2)
    dEVBM    = getarg     (4, dEVBM)
    WG_DOS   = getfloatarg(5, WG_DOS)
    Emin     = getfloatarg(6, Emin)
    Emax     = getfloatarg(7, Emax)
    EF0      = getfloatarg(8, EF0)
    if mode == 'BF':
        pass
    else:
        terminate("Error: Invalide mode [{}]".format(mode), usage = usage)


#def save_csv(path, headerlist, datalist, is_print = 0):

def BF_correction(mode, CAR_path1, CAR_path2):
    """
    概要:
        バンドフィリング補正を計算し、結果をファイルに保存し、グラフをプロットします。

    詳細説明:
        理想結晶と欠陥結晶のVASP計算結果ファイル（POSCAR, OUTCAR, DOSCAR, EIGENVAL）を読み込み、
        それぞれのバンド構造、DOS、フェルミ準位などの情報を取得します。
        dEVBMを用いて欠陥結晶のエネルギー準位を補正した後、バンドフィリング補正量（dEh, dEe, dEtot）を算出します。
        計算結果はサマリーファイルに保存され、DOSとバンド構造のグラフが生成されます。

    引数:
        :param mode: 実行モード。現在 BF のみがサポートされています。
        :type mode: str
        :param CAR_path1: 理想結晶のVASP計算結果ディレクトリへのパス。
        :type CAR_path1: str
        :param CAR_path2: 欠陥結晶のVASP計算結果ディレクトリへのパス。
        :type CAR_path2: str

    戻り値:
        :returns: なし
        :rtype: None
    """
    global dEVBM, EF0

    debug = 0
    
    vasp1 = tkVASP()
    vasp2 = tkVASP()

    base_path1 = vasp1.getdir(CAR_path1)
    base_path2 = vasp2.getdir(CAR_path2)

    print("")
    print(f"Summary file: {summaryfile}")
    INCAR_path1    = vasp1.get_INCAR(base_path1)
    POSCAR_path1   = vasp1.get_POSCAR(base_path1)
    OUTCAR_path1   = vasp1.get_OUTCAR(base_path1)
    DOSCAR_path1   = vasp1.get_VASPPath(base_path1, 'DOSCAR')
    EIGENVAL_path1 = vasp1.get_VASPPath(base_path1, 'EIGENVAL')
    INCAR_path2    = vasp2.get_INCAR(base_path2)
    POSCAR_path2   = vasp2.get_POSCAR(base_path2)
    OUTCAR_path2   = vasp2.get_OUTCAR(base_path2)
    DOSCAR_path2   = vasp2.get_VASPPath(base_path2, 'DOSCAR')
    EIGENVAL_path2 = vasp2.get_VASPPath(base_path2, 'EIGENVAL')
    print("CAR dir(ideal) : ", CAR_path1)
    print("  INCAR   : ", INCAR_path1)
    print("  POSCAR  : ", POSCAR_path1)
    print("  OUTCAR  : ", OUTCAR_path1)
    print("  DOSCAR  : ", DOSCAR_path1)
    print("  EIGENVAL: ", EIGENVAL_path1)
    print("CAR dir(defect): ", CAR_path2)
    print("  INCAR   : ", INCAR_path2)
    print("  POSCAR  : ", POSCAR_path2)
    print("  OUTCAR  : ", OUTCAR_path2)
    print("  DOSCAR  : ", DOSCAR_path2)
    print("  EIGENVAL: ", EIGENVAL_path2)

    print("")
    print("dEVBM: {} eV".format(dEVBM))
    print("")
    print("DOS plot")
    print(f"  FWHM of Gauss function for convolution: {WG_DOS} eV")
#    print("")
#    print("EF0: {} eV".format(EF0))
    
    dEVBM = pfloat(dEVBM, '')
    if dEVBM == '':
        app.terminate("Error: dEVBM value must be given", pause = True)

    print("")
    print("Ideal crystal model:")
    print("Read crystal structure from [{}]".format(POSCAR_path1))
    cry1 = vasp1.read_poscar(POSCAR_path1)
    if cry1 is None:
        terminate("Error: Can not read [{}]".format(POSCAR_path1), usage = usage)

    a1, b1, c1, alpha1, beta1, gamm1 = cry1.LatticeParameters()
    Vcell1        = cry1.Volume()
    types1        = cry1.AtomTypeList()
    ntypes1       = len(types1)
    sites1        = cry1.ExpandedAtomSiteList()
    nsites1       = len(sites1)
    atomnamelist1 = cry1.get_atom_name_list(mode = 'all', NameOnly = True)
    poslist1      = cry1.get_position_list(mode = 'all', IsReduce01 = True)
    cry1.PrintInf("cell")

    print("")
    print("Read OUTCAR from [{}]".format(OUTCAR_path1))
    outcarinf1 = vasp1.read_outcar_inf(OUTCAR_path1)
    if outcarinf1 is None:
        terminate("Error: Can not read [{}]".format(OUTCAR_path1), usage = usage)
    EF1 = vasp1.outcarinf["EF"]
    ISPIN1 = outcarinf1["ISPIN"]
#    ETOT1  = outcarinf1["TOTEN"]
    ETOT1  = outcarinf1["TOTEN_zero_sigma"]
    print(f"ISPIN={ISPIN1}")
    print(f"ETOT(sigma->0)={ETOT1}")

    print("")
    print("Defect crystal model:")
    print("Read from [{}]".format(POSCAR_path2))
    cry2 = vasp2.read_poscar(POSCAR_path2)
    if cry2 is None:
        terminate("Error: Can not read [{}]".format(POSCAR_path2), usage = usage)

    a2, b2, c2, alpha2, beta2, gamm2 = cry2.LatticeParameters()
    Vcell2        = cry2.Volume()
    types2        = cry2.AtomTypeList()
    ntypes2       = len(types2)
    sites2        = cry2.ExpandedAtomSiteList()
    nsites2       = len(sites2)
    atomnamelist2 = cry2.get_atom_name_list(mode = 'all', NameOnly = True)
    poslist2      = cry2.get_position_list(mode = 'all', IsReduce01 = True)
    cry2.PrintInf("cell")

    print("")
    print("Read OUTCAR from [{}]".format(OUTCAR_path2))
    outcarinf2 = vasp2.read_outcar_inf(OUTCAR_path2)
    if outcarinf2 is None:
        terminate("Error: Can not read [{}]".format(OUTCAR_path2), usage = usage)
    EF2 = vasp2.outcarinf["EF"]
    ISPIN2 = outcarinf2["ISPIN"]
#    ETOT2  = outcarinf2["TOTEN"]
    ETOT2  = outcarinf2["TOTEN_zero_sigma"]
    print(f"ISPIN={ISPIN2}")
    print(f"ETOT(sigma->0)={ETOT2}")

# Estimate supercell size
    print("")
    nx = int(a2 / a1 + 0.2)
    ny = int(b2 / b1 + 0.2)
    nz = int(c2 / c1 + 0.2)
    print("Supercell multipliers of the defect crystal with respect to the ideal crystal")
    print("  (nx, ny, nz) = ({}, {}, {})".format(nx, ny, nz))

    print("")
    print(f"Ideal crystal model: Read DOSCAR from [{DOSCAR_path1}]")
    doscarinf1 = vasp1.read_doscar(DOSCAR_path1, EF = EF1, normalize_E = False, unit = '/cm3')
    nE1 = doscarinf1["nE"]
    E1  = doscarinf1["E"]
    if nE1 == 0:
        Edmin1 = 0.0
        Edmax1 = 0.0
    else:
        Edmin1  = min(E1)
        Edmax1  = max(E1)
    tDOSup1 = doscarinf1["TotalDOSup"]
    Neup1   = doscarinf1["Neup"]
    tDOSdn1 = doscarinf1["TotalDOSdn"]
    Nedn1   = doscarinf1["Nedn"]
    if nE1 > 0:
        print(f"EF={EF1} eV")
        print(f"DOS E range: {Edmin1:10.6g} - {Edmax1:10.6g} eV, {nE1} points")
        print(f"  Energy is not normalized")
    else:
        print("")
        app.terminate(f"Error: Zero nE for DOS(E) in {DOSCAR_path1}", pause = True)

    print(f"Read EIGENVAL from [{DOSCAR_path1}]")
    eigenvalinf1 = vasp1.read_eigenval(EIGENVAL_path1, EF = EF1, normalize_E = False)
    nk1      = eigenvalinf1["nk"]
    nLevels1 = eigenvalinf1["nLevels"]
    bandedgeinf1a = vasp1.find_band_edges_from_eigenval(EF0 = EF0, eigenvalinf = eigenvalinf1, ISPIN = ISPIN1, occ_th = occ_th)
    bandedgeinf1b = vasp1.gbandedges(base_path1, OUTCAR_path1, EIGENVAL_path1, EF1)
    print("k points in EIGENVAL:")
    print("nk=", nk1)
    print("nLevels=", nLevels1)
    """
    for i in range(nk1):
        el = eigenvalinf1['EList'][i]
        kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = el
        print("  ({:8.4f} {:8.4f} {:8.4f}) w={:8.4g}  {:8.4f} {:12.4f}".format(kx, ky, kz, wk, dk, ktot))
    """

    EV1a   = bandedgeinf1a["EV"]
    EC1a   = bandedgeinf1a["EC"]
    Eg1a   = bandedgeinf1a["Eg"]
    EV1    = bandedgeinf1b["EVBM"]
    EC1    = bandedgeinf1b["ECBM"]
    Eg1    = bandedgeinf1b["Eg"]
    EHOMO1 = bandedgeinf1a["EHOMO"]
    ELUMO1 = bandedgeinf1a["ELUMO"]
    print(f"Band edge from EIGENVAL:")
    print(f"EF0={EF0:10.6f} eV")
    print(f"  find_band_edges: EV={EV1a:10.6f}  EC={EC1a:10.6f}  Eg={Eg1a:10.6f} eV")
    print(f"             HOMO:{EHOMO1:10.6f} eV")
    print(f"             LUMO:{ELUMO1:10.6f} eV")
    print(f"  gbandedges()   : EV={EV1:10.6f}  EC={EC1:10.6f}  Eg={Eg1:10.6f} eV")

    print("")
    print(f"Defect crystal model: Read DOSCAR from [{DOSCAR_path2}]")
    doscarinf2 = vasp2.read_doscar(DOSCAR_path2, EF = EF2, normalize_E = False, unit = '/cm3')
    nE2 = doscarinf2["nE"]
    E2  = doscarinf2["E"]
    if nE2 == 0:
        Edmin2 = 0.0
        Edmax2 = 0.0
    else:
        Edmin2  = min(E2)
        Edmax2  = max(E2)
    tDOSup2 = doscarinf2["TotalDOSup"]
    Neup2   = doscarinf2["Neup"]
    tDOSdn2 = doscarinf2["TotalDOSdn"]
    Nedn2   = doscarinf2["Nedn"]
    if nE2 > 0:
        print(f"EF={EF2} eV")
        print(f"DOS E range: {Edmin2:10.6g} - {Edmax2:10.6g} eV, {nE2} points")
        print(f"  Energy is not normalized")

    print(f"Read EIGENVAL from [{DOSCAR_path2}]")
    eigenvalinf2 = vasp2.read_eigenval(EIGENVAL_path2, EF = EF2, normalize_E = False)
    nk2      = eigenvalinf2["nk"]
    nLevels2 = eigenvalinf2["nLevels"]
    bandedgeinf2a = vasp2.find_band_edges_from_eigenval(EF0 = EF0, eigenvalinf = eigenvalinf2, ISPIN = ISPIN2, occ_th = occ_th)
    bandedgeinf2b = vasp2.gbandedges(base_path2, OUTCAR_path2, EIGENVAL_path2, EF2)
    print("k points in EIGENVAL:")
    print("nk=", nk2)
    print("nLevels=", nLevels2)
    for i in range(nk2):
        el = eigenvalinf2['EList'][i]
        kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = el
        print("  ({:8.4f} {:8.4f} {:8.4f}) w={:8.4g}  {:8.4f} {:12.4f}".format(kx, ky, kz, wk, dk, ktot))

    EV2a   = bandedgeinf1a["EV"]
    EC2a   = bandedgeinf1a["EC"]
    Eg2a   = bandedgeinf1a["Eg"]
    EV2    = bandedgeinf1b["EVBM"]
    EC2    = bandedgeinf1b["ECBM"]
    Eg2    = bandedgeinf1b["Eg"]
    EHOMO2 = bandedgeinf1a["EHOMO"]
    ELUMO2 = bandedgeinf1a["ELUMO"]
    print(f"Band edge from EIGENVAL:")
    print(f"EF0={EF0:10.6f} eV")
    print(f"  find_band_edges: EV={EV2a:10.6f}  EC={EC2a:10.6f}  Eg={Eg2a:10.6f} eV")
    print(f"             HOMO:{EHOMO2:10.6f} eV")
    print(f"             LUMO:{ELUMO2:10.6f} eV")
    print(f"  gbandedges()   : EV={EV2:10.6f}  EC={EC2:10.6f}  Eg={Eg2:10.6f} eV")

    print("")
    print("DOS plot range")
    print(f"  Ideal crystal model:")
    print(f"    EF={EF1} eV")
    print(f"    E range: {Edmin1:10.6g} - {Edmax1:10.6g} eV, {nE1} points")
    """
    print(f"    DOS(E) is scaled for the size of the defect model by {nx}x{ny}x{nz}")
    k = nx * ny * nz * 1.0
    for i in range(len(E1)):
        tDOSup1[i] *= k
        Neup1[i]   *= k
        if ISPIN1 == 2:
            tDOSdn2[i] *= k
            Nedn1[i]   *= k
    """

    print(f"  Defect crystal model:")
    print(f"    Original EF={EF2} eV")
    print(f"    Original E range: {Edmin2:10.6g} - {Edmax2:10.6g} eV, {nE2} points")
    EF2 -= dEVBM
    Edmin2 -= dEVBM
    Edmax2 -= dEVBM
    for i in range(len(E2)):
        E2[i] -= dEVBM
    for ik in range(nk2):
        el = eigenvalinf2['EList'][ik]
        kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = [*el]
        Eups[ik] -= dEVBM
        if ISPIN2 == 2:
            Edns[ik] -= dEVBM
    print(f"    Adjusted EF={EF2} eV")
    print(f"    Adjusted E range: {Edmin2:10.6g} - {Edmax2:10.6g} eV, {nE2} points")

    print("")
    print("Calculate band filling energy")
    totalw  = 0.0
    totalNe = 0.0
    dEe = 0.0
    dEh = 0.0
    for ik in range(nk2):
        el = eigenvalinf2['EList'][ik]
        kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = [*el]
        totalw += wk
        print("  k[{:3d}]: ({:8.4f} {:8.4f} {:8.4f})  w={:8.4g}".format(ik, kx, ky, kz, wk))
        for iE in range(nLevels2):
            Eup  = Eups[iE]
            neup = occups[iE]
            totalNe += neup * wk

            if Eup <= EV1 and abs(neup - Nemax) > 1.0e-6:
                dEh += wk * (Nemax - neup) * (Eup - EV1)
                print(f"    dEh: up iE={iE:3d} E={Eup:10.6g} eV dNh={Nemax - neup:10.4g} dEh={dEh:12.6g} eV")
            if Eup >= EC1 and neup > 1.0e-6:
                dEe += wk * neup * (Eup - EC1)
                print(f"    dEe: up iE={iE:3d} E={Eup:10.6g} eV dNe={neup:10.4g} dEe={dEe:12.6g} eV")
            if ISPIN2 == 2:
                Edn  = Edns[iE]
                nedn = occdns[iE]
                totalNe += nedn * wk

                if Edn <= EV1 and abs(nedn - Nemax) > 1.0e-6:
                    dEh += wk * (Nemax - nedn) * (Edn - EV1)
                    print(f"    dEh: dn iE={iE:3d} E={Edn:10.6g} eV dNh={Nemax - nedn:10.4g} dEh={dEh:12.6g} eV")
                if Edn >= EC1 and nedn > 1.0e-6:
                    dEe += wk * nedn * (Edn - EC1)
                    print(f"    dEe: dn iE={iE:3d} E={Edn:10.6g} eV dNe={nedn:10.4g} dEe={dEe:12.6g} eV")

    print("")
    print("From defect crystal model:")
    print("  Total weight: ", totalw)
    print("  NELECT      : ", outcarinf2["NELECT"])
    dEtot = dEh + dEe
    if ISPIN2 == 1:
        totalNe *= 2.0
        dEh   *= 2.0
        dEe   *= 2.0
        dEtot *= 2.0
    print("  Count spin degeneracy of 2:")
    print("    Total Ne    : ", totalNe)
    print("    dEh   = {:12.6g} eV".format(dEh))
    print("    dEe   = {:12.6g} eV".format(dEe))
    print("    dEtot = {:12.6g} eV".format(dEtot))

    print("")
    print(f"Save summary to {summaryfile}")
    ini = tkIniFile()
    ini.write_from_scratch(summaryfile, 'results',
                Car_dir1 = base_path1,
                Car_dir2 = base_path2,
                Etot1 = ETOT1 * nx * ny * nz,
                Etot2 = ETOT2,
                totalNe = totalNe,
                dEh = dEh,
                dEe = dEe,
                dEtot = dEtot
                )


#=============================
# Plot graphs
#=============================
    print("")

    fig = plt.figure(figsize = figsize)
    ax1    = fig.add_subplot(2, 2, 1)
    ax2    = fig.add_subplot(2, 2, 3)
    axband = fig.add_subplot(1, 2, 2)

    print("E1=", E1)
    tDOSup_s1 = convolution(E1, tDOSup1, WG_DOS, func_type = 'gauss')
    ax1.plot(E1, tDOSup_s1, label = 'up(ideal)', linestyle = '-', linewidth = 0.5, color = 'black')
    if ISPIN1 == 2:
        tDOSdn_s1 = convolution(E, tDOSdn1, WG_DOS, func_type = 'gauss')
        ax1.plot(E1, tDOSdn_s1, label = 'dn(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'black')
    ylim = ax1.get_ylim()
    ax1.plot([EV1, EV1],       ylim, label = '$E_V$(ideal)',      linestyle = 'dashed', linewidth = 0.5, color = 'red')
    ax1.plot([EC1, EC1],       ylim, label = '$E_C$(ideal)',      linestyle = 'dashed', linewidth = 0.5, color = 'red')
    ax1.plot([EF1, EF1],       ylim, label = '$E_F$(ideal)',      linestyle = 'dashed', linewidth = 0.5, color = 'blue')
#    ax1.plot([EF0, EF0],     ylim, label = '$E_{F0}$',   linestyle = 'dashed', linewidth = 0.5, color = 'blue')
#    ax1.plot([EHOMO1, EHOMO1], ylim, label = '$E_{HOMO}$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'green')
#    ax1.plot([ELUMO1, ELUMO1], ylim, label = '$E_{LUMO}$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'green')
#    ax1.set_xlabel("$E$ (eV)",    fontsize = fontsize)
    ax1.set_ylabel("DOS (states/cm$^3$/eV)", fontsize = fontsize)
#    ax1.set_xlim([Emin, Emax])
    ax1.tick_params(labelsize = fontsize)

    tDOSup_s2 = convolution(E2, tDOSup2, WG_DOS, func_type = 'gauss')
    ax1.plot(E2, tDOSup_s2, label = 'up(defect)', linestyle = '-', linewidth = 1.0, color = 'red')
    if ISPIN2 == 2:
        tDOSdn_s2 = convolution(E, tDOSdn2, WG_DOS, func_type = 'gauss')
        ax1.plot(E2, tDOSdn_s2, label = 'dn(defect)', linestyle = 'dashed', linewidth = 1.0, color = 'red')
    ax1.plot([EF2, EF2],     ylim, label = '$E_F$(defect)',      linestyle = 'dashed', linewidth = 1.0, color = 'blue')
    ax1.tick_params(labelsize = fontsize)
    ax1.set_title(f"{CAR_path1}\n{CAR_path2}", fontsize = titlefontsize)
    ax1.legend(fontsize = legend_fontsize)

    ax2.plot(E2, Neup2, label = 'up(defect)', linestyle = '-', linewidth = 0.5, color = 'black')
    if ISPIN2 == 2:
        ax2.plot(E2, Nedn2, label = 'dn(defect)', linestyle = '-', linewidth = 0.5, color = 'red')
    ylim = ax2.get_ylim()
    ax2.plot([EV1, EV1],   ylim,     label = '$E_V$(ideal)',      linestyle = 'dashed', linewidth = 0.5, color = 'red')
    ax2.plot([EC1, EC1],   ylim,     label = '$E_C$(ideal)',      linestyle = 'dashed', linewidth = 0.5, color = 'red')
#    ax2.plot([EF0, EF0], ylim,     label = '$E_{F0}$',   linestyle = 'dashed', linewidth = 0.5, color = 'blue')
#    ax2.plot([EHOMO1, EHOMO1], ylim, label = '$E_{HOMO}$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'green')
#    ax2.plot([ELUMO1, ELUMO1], ylim, label = '$E_{LUMO}$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'green')
    ax2.plot([EF1, EF1],   ylim, label = '$E_F$(ideal)',      linestyle = 'dashed', linewidth = 0.5, color = 'blue')
    ax2.plot([EF2, EF2],   ylim, label = '$E_F$(defect)',      linestyle = 'dashed', linewidth = 1.0, color = 'blue')
    ax2.set_xlabel("$E$ (eV)",    fontsize = fontsize)
    ax2.set_ylabel("$N_e$ (states/unit cell)", fontsize = fontsize)
#    ax2.set_xlim([Emin, Emax])
    ax2.legend(fontsize = legend_fontsize)
    ax2.tick_params(labelsize = fontsize)

    xk   = []
    yEup = []
    yNup = []
    yEdn = []
    yNdn = []
    for i in range(nk2):
        el = eigenvalinf2['EList'][i]
        kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = [*el]
        xk.append(ktot)
        yEup.append(Eups)
        yNup.append(occups)
        if ISPIN2 == 2:
            yEdn.append(Edns)
            yNdn.append(occdns)

    plot_band(plt, axband, ISPIN2, xk, yEup, yEdn, [Emin, Emax], occups = yNup, occdns = yNdn, 
                ktotallist = None, ktotal_namelist = None, 
                EV = EV1, EC = EC1, EF = EF2, EHOMO = None, ELUMO = None,
                EFlabel = '$E_{F}$(defect)',
                markersize = band_marker_size, markeredgewidth = band_marker_edge_width)
    axband.set_title(f"dEh={dEh:12.6g} eV\ndEe={dEe:12.6g} eV\ndEtot={dEtot:12.6g} eV", fontsize = titlefontsize)

    plt.tight_layout()


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

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

    app.terminate("", pause = True)


def main():
    """
    概要:
        スクリプトのメインエントリポイント。引数の処理とバンドフィリング補正の実行を行います。

    詳細説明:
        コマンドライン引数をパースしてグローバル変数を更新し、ログファイルをセットアップします。
        指定されたモードに応じてBF_correction関数を呼び出し、バンドフィリング補正計算を実行します。

    戻り値:
        :returns: なし
        :rtype: None
    """
    global mode
    global cifpath, poscarpath
    global summaryfile
    
    updatevars()

    vasp = tkVASP()
    base_path = vasp.getdir(CAR_dir2)
    logfile     = os.path.join(base_path, 'BF_correction-out.txt')
    summaryfile = os.path.join(base_path, 'BF_correction-summary.prm')
    print("")
    print(f"Open logfile [{logfile}]")
    app.redirect(targets = ["stdout", logfile], mode = 'w')

    print("")
    print("=============== Estimate band filling correction for defect model calculated by VASP ============")
    print("")
    print("mode: ", mode)

    if mode == 'BF':
        BF_correction(mode, CAR_dir1, CAR_dir2)
    else:
        terminate(f"Error: Invalide mode [{mode}]", usage = usage)


if __name__ == "__main__":
    main()