"""
vasp_correction_bandfilling.py: VASP計算結果から欠陥モデルのバンドフィリング補正を推定します。
詳細説明:
本スクリプトは、欠陥モデルの形成エネルギー計算において必要となるバンドフィリング補正 (dE_BF) を算出します。
理想結晶と欠陥結晶のVASP計算結果 (OUTCAR, DOSCAR, EIGENVAL) を比較し、
dEVBMによるバンドアラインメント後の状態占有率から電子と正孔によるエネルギー補正量を評価します。
また、計算結果のDOSとバンド構造をグラフとして可視化し、サマリーファイルに結果を保存します。
関連リンク:
:doc:`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: 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: 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: str: 実行モード。現在 'BF' のみがサポートされています。
:param CAR_path1: str: 理想結晶のVASP計算結果ディレクトリへのパス。
:param CAR_path2: str: 欠陥結晶のVASP計算結果ディレクトリへのパス。
:returns: 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: 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()