"""
概要: VASPのDOSCARファイルに基づきキャリア密度とフェルミ準位を計算するスクリプト。
詳細説明: このスクリプトは、VASP計算から得られたDOSCAR、OUTCAR、EIGENVAL、POSCARファイル、
および欠陥形成エネルギーのExcelファイルを利用して、半導体の電子・正孔キャリア密度、
欠陥密度、およびフェルミ準位の温度依存性またはフェルミ準位依存性を計算・解析します。
特に、温度変化に伴うフェルミ準位の自己無撞着計算や、欠陥の電荷状態による形成エネルギーの変化、
遷移準位の特定、およびこれらの特性のグラフ表示をサポートします。
関連リンク: :doc:`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: dict: 欠陥のプロパティを格納するキーワード引数辞書。
- "labels": list[str]: ラベルのリスト。
- "plabels": list[str]: プロットラベルのリスト。
- "names": list[str]: 欠陥名のリスト。
- "atom_label": str: 原子ラベル。
- "site_label": str: サイトラベル。
- "charge_label": str: 電荷ラベル。
- "entropy_label": str: エントロピーラベル。
- "N0_label": str: N0ラベル。
- "Ndoped_label": str: Ndopedラベル。
- "dH0_labels": list[str]: 形成エンタルピーラベルのリスト。
- "atom": str: 原子タイプ。
- "site": str: サイトタイプ。
- "name": str: 欠陥名。
- "charge": int: 欠陥の電荷。
- "entropy": float: 欠陥のエントロピー項 (S/kB)。
- "N0": float: 欠陥サイトの基準数。
- "Ndoped": float: ドープされた欠陥の濃度。
- "dH0s": list[float]: 各ポイントでの欠陥形成エンタルピーのリスト。
"""
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: str: 欠陥サイトの名前。
"""
return self._site
[ドキュメント]
class tkDefects():
"""
概要: 複数の欠陥情報を管理するクラス。
詳細説明: Excelファイルから欠陥形成エネルギーを読み込み、セル体積を考慮して
欠陥密度、キャリア密度、フェルミ準位の計算を行います。
また、欠陥の遷移準位の探索と保存も担当します。
"""
def __init__(self, **args):
"""
概要: tkDefectsクラスのコンストラクタ。
詳細説明: 欠陥に関する各種リストやパスを初期化します。
:param args: dict: 欠陥のプロパティを格納するキーワード引数辞書。
- "labels": list[str]: ラベルのリスト。
- "plabels": list[str]: プロットラベルのリスト。
- "atoms": list[str]: 原子タイプのリスト。
- "sites": list[str]: サイトタイプのリスト。
- "names": list[str]: 欠陥名のリスト。
"""
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: list[str]: ユニークな欠陥名のリスト。
"""
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:
"""
概要: 特定の欠陥の形成エネルギー（$\Delta$E）を計算します。
詳細説明: 指定されたポイント、欠陥ID、温度、フェルミ準位に基づいて、
欠陥の形成エンタルピーと電荷によるエネルギー項を合計します。
:param iPoint: int: 計算に使用する形成エンタルピーのインデックス。
:param idefect: int: 欠陥リスト内の欠陥のインデックス。
:param T: float: 温度 (K)。
:param EF: float: フェルミ準位 (eV)。
:returns: float: 計算された欠陥形成エネルギー ($\Delta$E) (eV)。
"""
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) とフェルミ準位で、全欠陥の密度を計算します。
詳細説明: 各欠陥サイトの分配関数 (ZS) を考慮し、個々の欠陥の電荷状態に応じた密度と、
全欠陥による総電荷 (Qtot) を算出します。
:param Tdef: float: 欠陥が凍結している温度 (K)。
:param EF: float: フェルミ準位 (eV)。
:param ECmin: float: 伝導帯のエネルギー下限 (eV)。
:param ECmax: float: 伝導帯のエネルギー上限 (eV)。
:param EVmin: float: 価電子帯のエネルギー下限 (eV)。
:param EVmax: float: 価電子帯のエネルギー上限 (eV)。
:returns: tuple[dict, list, dict, list, float]:
- ZS: dict: 各サイトの分配関数。
- Nds: list[float]: 各欠陥状態の密度 (cm^-3)。
- NdsTdef: dict: Tdefにおける総欠陥密度を電荷状態とサイト名でキー付けした辞書 (cm^-3)。
- dEs: list[float]: 各欠陥状態の形成エネルギー (eV)。
- Qtot: float: 全欠陥による総電荷 (単位電荷)。
"""
# 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]:
"""
概要: 電子温度 (Te) と凍結された欠陥密度 (NdsTdef) を用いて、電荷を帯びた欠陥密度を計算します。
詳細説明: まず、各サイトの全欠陥密度を合計し、電子温度での各電荷状態の分配関数 (ZSA) を計算します。
その後、各電荷状態の欠陥密度と、全欠陥による総電荷 (Qtot) を算出します。
:param Te: float: 電子温度 (K)。
:param NdsTdef: dict: 欠陥凍結温度Tdefで計算された総欠陥密度の辞書 (cm^-3)。キーは"{atom}_{site}^{charge}"。
:param EF: float: フェルミ準位 (eV)。
:param ECmin: float: 伝導帯のエネルギー下限 (eV)。
:param ECmax: float: 伝導帯のエネルギー上限 (eV)。
:param EVmin: float: 価電子帯のエネルギー下限 (eV)。
:param EVmax: float: 価電子帯のエネルギー上限 (eV)。
:returns: tuple[dict, list, dict, list, float]:
- ZSA: dict: 各サイトの（電子温度における）分配関数。
- Nds: list[float]: 各欠陥状態の密度 (cm^-3)。
- NdsTdef: dict: Tdefにおける総欠陥密度を電荷状態とサイト名でキー付けした辞書 (cm^-3)。
- dEs: list[float]: 各欠陥状態の形成エネルギー (eV)。
- Qtot: float: 全欠陥による総電荷 (単位電荷)。
"""
# 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: float: 電子温度 (K)。
:param Tdef: float: 欠陥が凍結している温度 (K)。
:param NdsTdef: dict | None: 欠陥凍結温度Tdefで計算された総欠陥密度の辞書。Noneの場合、未凍結として計算。
:param EF: float: フェルミ準位 (eV)。
:param ECmin: float: 伝導帯のエネルギー下限 (eV)。
:param ECmax: float: 伝導帯のエネルギー上限 (eV)。
:param EVmin: float: 価電子帯のエネルギー下限 (eV)。
:param EVmax: float: 価電子帯のエネルギー上限 (eV)。
:returns: tuple[dict, list, dict, list, float]:
- Z: dict: 各サイトの分配関数（またはZSA）。
- Nds: list[float]: 各欠陥状態の密度 (cm^-3)。
- NdsTdef: dict: Tdefにおける総欠陥密度を電荷状態とサイト名でキー付けした辞書 (cm^-3)。
- dGs: list[float]: 各欠陥状態の形成エネルギー (eV)。
- Qtot: float: 全欠陥による総電荷 (単位電荷)。
"""
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]:
"""
概要: 電子、正孔、および欠陥のキャリア密度と総電荷を計算します。
詳細説明: 電子温度 (T) が欠陥凍結温度 (Tdef) 以下の場合、欠陥密度はTdefで凍結されます。
TがTdefより高い場合、欠陥密度はTで再計算されます。
最終的に電子密度 (ne)、正孔密度 (nh)、各欠陥密度 (Nds)、および総電荷 (Qtot) を返します。
:param T: float: 電子温度 (K)。
:param Tdef: float: 欠陥が凍結している温度 (K)。
:param NdsTdef: dict | None: 欠陥凍結温度Tdefで計算された総欠陥密度の辞書。Noneの場合、未凍結として計算。
:param EF: float: フェルミ準位 (eV)。
:param ECmin: float: 伝導帯のエネルギー下限 (eV)。
:param ECmax: float: 伝導帯のエネルギー上限 (eV)。
:param EVmin: float: 価電子帯のエネルギー下限 (eV)。
:param EVmax: float: 価電子帯のエネルギー上限 (eV)。
:returns: tuple[dict, float, float, list, dict, list, float]:
- Z: dict: 各サイトの分配関数。
- ne: float: 電子密度 (cm^-3)。
- nh: float: 正孔密度 (cm^-3)。
- Nds: list[float]: 各欠陥状態の密度 (cm^-3)。
- NdsTdef: dict: Tdefにおける総欠陥密度を電荷状態とサイト名でキー付けした辞書 (cm^-3)。
- dGs: list[float]: 各欠陥状態の形成エネルギー (eV)。
- Qtot: float: 全電荷 (単位電荷)。
"""
# 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:
"""
概要: 全電荷の中性条件からの偏差 ($\Delta$Q) を計算します。
詳細説明: calculate_densitiesメソッドを呼び出して、電子、正孔、欠陥の密度から総電荷を計算し、
その値を返します。この関数はフェルミ準位探索のための目的関数として使用されます。
:param T: float: 電子温度 (K)。
:param Tdef: float: 欠陥が凍結している温度 (K)。
:param NdsTdef: dict | None: 欠陥凍結温度Tdefで計算された総欠陥密度の辞書。
:param EF: float: フェルミ準位 (eV)。
:param ECmin: float: 伝導帯のエネルギー下限 (eV)。
:param ECmax: float: 伝導帯のエネルギー上限 (eV)。
:param EVmin: float: 価電子帯のエネルギー下限 (eV)。
:param EVmax: float: 価電子帯のエネルギー上限 (eV)。
:returns: float: 全電荷 ($\Delta$Q) (単位電荷)。
"""
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:
"""
概要: 全電荷の中性条件からの偏差 ($\Delta$Q) のフェルミ準位 (EF) による数値微分を計算します。
詳細説明: 中心差分法を用いてdQ/dEFを計算します。これはニュートン法のヤコビアンとして使用されます。
:param h: float: 微小なEFの摂動量 (eV)。
:param T: float: 電子温度 (K)。
:param Tdef: float: 欠陥が凍結している温度 (K)。
:param NdsTdef: dict | None: 欠陥凍結温度Tdefで計算された総欠陥密度の辞書。
:param EF: float: フェルミ準位 (eV)。
:param ECmin: float: 伝導帯のエネルギー下限 (eV)。
:param ECmax: float: 伝導帯のエネルギー上限 (eV)。
:param EVmin: float: 価電子帯のエネルギー下限 (eV)。
:param EVmax: float: 価電子帯のエネルギー上限 (eV)。
:returns: float: $\Delta$QのEFによる数値微分 ($\partial Q / \partial E_F$)。
"""
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]:
"""
概要: 全電荷中性条件を満たすフェルミ準位 (EF) を探索します。
詳細説明: 二分法とニュートン法を組み合わせて、全電荷 ($\Delta$Q) がゼロとなるEFを探索します。
最初に広い範囲で二分法を実行し、その後、より精密なニュートン法を適用します。
:param T: float: 電子温度 (K)。
:param Tdef: float: 欠陥が凍結している温度 (K)。
:param NdsTdef: dict | None: 欠陥凍結温度Tdefで計算された総欠陥密度の辞書。
:param ECmin: float: 伝導帯のエネルギー下限 (eV)。
:param ECmax: float: 伝導帯のエネルギー上限 (eV)。
:param EVmin: float: 価電子帯のエネルギー下限 (eV)。
:param EVmax: float: 価電子帯のエネルギー上限 (eV)。
:param callbackfunc: callable | None: 探索の各ステップで呼び出されるコールバック関数。
:param eps_bisec: float: 二分法の収束判定閾値。
:param nmaxiter_bisection: int: 二分法の最大反復回数。
:param initial: list[float]: 二分法の初期探索範囲 [EF_min, EF_max] (eV)。
:param eps_newton: float: ニュートン法の収束判定閾値。
:param nmaxiter_newton: int: ニュートン法の最大反復回数。
:param dump_newton: float: ニュートン法のダンプ係数。
:param h_newton: float: ニュートン法の数値微分のための微小量。
:param is_print: bool: 途中経過をコンソールに出力するかどうか。
:returns: tuple[float | None, float | None]: 収束したフェルミ準位 (EF) と、そのときの$\Delta$Qの値。
収束しない場合は(None, None)を返します。
"""
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]:
"""
概要: 特定の欠陥状態 (idx0) から他のすべての欠陥状態への遷移準位EFを計算します。
詳細説明: 欠陥形成エネルギーの交点から、異なる電荷状態間の遷移準位を特定します。
電荷が同じ場合は非常に大きな値 (1.0e100) を返します。
:param idx0: int: 基準となる欠陥状態のインデックス。
:param groupdata: list[dict]: 特定の欠陥名に属する全欠陥状態のデータリスト。
各辞書は'charge' (int) と 'dH0' (float) を含む。
:returns: list[float]: 各欠陥状態への遷移準位EFのリスト (eV)。
"""
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]:
"""
概要: 現在のフェルミ準位 (prevEF) より大きく、指定範囲内の次の遷移準位EFを探索します。
詳細説明: `find_all_transition_EF` を使用してすべての遷移準位を計算し、
現在のフェルミ準位よりも大きく、かつ探索範囲 (EFmin, EFmax) 内で最小の遷移準位を特定します。
:param idx0: int: 基準となる欠陥状態のインデックス。
:param prevEF: float: 現在のフェルミ準位 (eV)。
:param EFmin: float: フェルミ準位の最小探索値 (eV)。
:param EFmax: float: フェルミ準位の最大探索値 (eV)。
:param groupdata: list[dict]: 特定の欠陥名に属する全欠陥状態のデータリスト。
:returns: tuple[int | None, float]: 次の遷移準位のインデックスとEF値。見つからない場合は (None, EFmax) を返します。
"""
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]:
"""
概要: 特定のフェルミ準位 (EF) における、ある欠陥状態の形成エンタルピーの順位を計算します。
詳細説明: 同じ欠陥名に属するすべての電荷状態について、指定されたEFにおける形成エンタルピーを計算し、
入力された(EF, dH)のペアがその中で何番目に低い形成エンタルピーを持つか (iorder) を決定します。
また、同じEFで最も低い形成エンタルピーを持つ欠陥状態のインデックス (idxmin) とそのdH (dHmin) も返します。
:param groupdata: list[dict]: 特定の欠陥名に属する全欠陥状態のデータリスト。
:param idx: int: 現在の欠陥状態のインデックス。
:param EF: float: フェルミ準位 (eV)。
:param dH: float: 現在の欠陥状態の形成エンタルピー (eV)。
:returns: tuple[int, int, float]:
- iorder: int: 指定された(EF, dH)のペアがグループ内で何番目に低い形成エンタルピーを持つか (0-indexed)。
- idxmin: int: 同じEFで最も低い形成エンタルピーを持つ欠陥状態のインデックス。
- dHmin: float: 同じEFで最も低い形成エンタルピーの値 (eV)。
"""
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まで変化させながら、各電荷状態の形成エンタルピーを計算し、
形成エンタルピーが最小となる電荷状態が変化する遷移準位を特定します。
'allpts'はすべての交点、'minpts'は最小形成エンタルピー曲線上の主要な点（遷移準位点と区間の端点）を収集します。
:param name: str: 欠陥名。
:param groupdata: list[dict]: 特定の欠陥名に属する全欠陥状態のデータリスト。
:param EFmin: float: フェルミ準位の最小探索値 (eV)。
:param EFmax: float: フェルミ準位の最大探索値 (eV)。
:returns: tuple[list[list], list[list]]:
- allpts: list[list]: すべての遷移点 (EF, dH, q_prev, q_next, iorder, idxmin) のリスト。
- minpts: list[list]: 最小形成エンタルピー曲線上の主要な点 (EF, dH, q_prev, q_next, iorder, idxmin) のリスト。
"""
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]]]:
"""
概要: 全てのユニークな欠陥名に対して、遷移準位のリストを生成します。
詳細説明: 各欠陥名について`find_all_transitions`を呼び出し、
すべての形成エンタルピー曲線上の点と、最小曲線上の点を取得します。
:returns: tuple[list[str], list[list[list]], list[list[list]]]:
- names: list[str]: ユニークな欠陥名のリスト。
- allpts_list: list[list[list]]: 各欠陥名の全遷移点のリストのリスト。
- minpts_list: list[list[list]]: 各欠陥名の最小形成エンタルピー曲線上の点のリストのリスト。
"""
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ファイルに保存します。
詳細説明: 各欠陥の電荷状態について、フェルミ準位の範囲 (EFminからEFmax) における
形成エンタルピーの線形関係を計算し、Excelシートに書き出します。
:param path: str: 出力Excelファイルのパス。
:returns: int | None: ファイル書き込みが成功すれば1、失敗すれば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ファイルに保存します。
詳細説明: `find_all_transitions`で特定された、各欠陥タイプにおいてフェルミ準位に対して
形成エンタルピーが最小となる曲線上の遷移準位点とその情報をExcelシートに書き出します。
:param path: str: 出力Excelファイルのパス。
:returns: int | None: ファイル書き込みが成功すれば1、失敗すれば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]:
"""
概要: 指定された欠陥名に属する全電荷状態のデータを取得します。
詳細説明: 全ての欠陥から、指定された`name`と一致する欠陥（異なる電荷状態を含む）を抽出し、
それらの情報を辞書形式でリストにまとめます。リストは電荷の降順にソートされます。
:param name: str: 抽出する欠陥名。
:param iPoint: int: 使用する形成エンタルピーのインデックス。
:returns: list[dict]: 抽出された欠陥データのリスト。各辞書は'defects', 'atom', 'site', 'name', 'charge', 'dH0'を含む。
"""
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: float: ユニットセルの体積 (A^3)。
:returns: float: 設定されたユニットセルの体積 (A^3)。
"""
self.V = V
return self.V
[ドキュメント]
def Volume(self) -> float | None:
"""
概要: ユニットセルの体積を取得します。
:returns: float | None: 設定されたユニットセルの体積 (A^3)。設定されていない場合はNone。
"""
return self.V
[ドキュメント]
def add(self, defect: 'tkDefect'):
"""
概要: 欠陥オブジェクトをリストに追加します。
:param defect: tkDefect: 追加するtkDefectオブジェクト。
"""
self.defects.append(defect)
[ドキュメント]
def get(self, idx: int):
"""
概要: 指定されたインデックスの欠陥オブジェクトを取得します。
:param idx: int: 取得する欠陥のインデックス。
:returns: tkDefect | None: 指定されたインデックスの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):
"""
概要: 登録されている欠陥情報を標準出力または指定されたファイルポインタに出力します。
詳細説明: 欠陥の数、ポイント数、セル体積、および各欠陥の詳細情報（原子、サイト、電荷、エントロピー、N0、Ndoped、dH0s）を整形して表示します。
:param PrintLabels: bool: ラベル情報 (labels, plabels, names) を出力するかどうか。
:param outfp: object | None: 出力先のファイルポインタ。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ファイルを読み込み、tkDefectsオブジェクトを構築します。
詳細説明: 指定されたExcelファイルから欠陥の原子、サイト、電荷、エントロピー、基準サイト数、
ドープ濃度、および各ポイントでの形成エンタルピーを抽出します。
ファイル形式のバージョン（'Version'列の有無）を自動判別し、適切にデータをパースします。
セル体積Vが与えられた場合、それを設定します。
:param infile: str: 欠陥形成エンタルピーを含むExcelファイルのパス。
:param V: float | None: ユニットセルの体積 (A^3)。与えられない場合はNone。
:returns: tkDefects: 読み込まれたデータで構築された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):
"""
概要: スクリプトのコマンドライン引数と使用方法を表示します。
詳細説明: プログラムの実行に必要な引数とその例を示します。
特にモード ('EF' または 'T') に応じた引数の指定方法を説明します。
:param app: tkApplication | None: tkApplicationのインスタンス。現在のところ使用されていません。
"""
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` を解析し、実行モード ('EF'または'T') に応じて、
欠陥データパス、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: str: 出力CSVファイルのパス。
:param header: list[str]: CSVファイルのヘッダー行のリスト。
:param datalist: list[list[float]]: 書き込むデータ。`datalist[j][i]`がi行j列のデータとなるように渡されます。
"""
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軸データがその範囲内の行のみを抽出します。
:param infile: str: 入力CSVファイルのパス。
:param xmin: float | None: x軸データの最小値。この値以上のデータのみが抽出されます。Noneの場合は無視。
:param xmax: float | None: x軸データの最大値。この値以下のデータのみが抽出されます。Noneの場合は無視。
:param delimiter: str: CSVファイルの区切り文字。
:returns: tuple[list[str], list[float], list[float]]:
- header: list[str]: CSVファイルのヘッダー行のリスト。
- data0: list[float]: 抽出されたx軸データのリスト。
- data1: list[float]: 抽出されたy軸データのリスト。
"""
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:
"""
概要: 補間関数 (fdos) を使用して、指定されたエネルギー (E) における状態密度 (DOS) を返します。
詳細説明: DOSは予めVASPのDOSCARファイルから読み込まれ、`fdos`として補間関数が設定されています。
入力エネルギーがDOSの定義範囲外の場合、エラーメッセージを出力してプログラムを終了します。
:param E: float: 状態密度を評価するエネルギー (eV)。
:returns: float: 指定されたエネルギーにおける状態密度 (states/cm^3)。
"""
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:
"""
概要: フェルミ・ディラック分布関数を計算します (電子占有確率)。
詳細説明: 指定されたエネルギー、温度、フェルミ準位において、
電子がそのエネルギー準位を占有する確率を計算します。
温度が0Kの場合、ステップ関数として動作します。
:param E: float: エネルギー (eV)。
:param T: float: 温度 (K)。
:param EF: float: フェルミ準位 (eV)。
:returns: 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:
"""
概要: フェルミ・ディラック分布関数に基づいて正孔の占有確率を計算します。
詳細説明: `fe(E, T, EF)` (電子占有確率) を用いて、正孔の占有確率 (1 - `fe`) を計算します。
:param E: float: エネルギー (eV)。
:param T: float: 温度 (K)。
:param EF: float: フェルミ準位 (eV)。
:returns: float: 正孔占有確率。
"""
return 1.0 - fe(E, T, EF)
# Define the function to be integrated
[ドキュメント]
def DOSfe(E: float, T: float, EF: float) -> float:
"""
概要: 状態密度 (DOS) と電子のフェルミ・ディラック分布関数の積を計算します。
詳細説明: これは電子密度を計算するための積分の被積分関数となります。
:param E: float: エネルギー (eV)。
:param T: float: 温度 (K)。
:param EF: float: フェルミ準位 (eV)。
:returns: float: DOS(E) * fe(E, T, EF) の値。
"""
return DOS(E) * fe(E, T, EF)
[ドキュメント]
def DOSfh(E: float, T: float, EF: float) -> float:
"""
概要: 状態密度 (DOS) と正孔のフェルミ・ディラック分布関数の積を計算します。
詳細説明: これは正孔密度を計算するための積分の被積分関数となります。
:param E: float: エネルギー (eV)。
:param T: float: 温度 (K)。
:param EF: float: フェルミ準位 (eV)。
:returns: float: DOS(E) * fh(E, T, EF) の値。
"""
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: callable: 積分する関数。引数としてエネルギー (E) を取り、値を返します。
:param E0: float: 積分範囲の開始エネルギー (eV)。
:param E1: float: 積分範囲の終了エネルギー (eV)。
:param h: float: 積分刻み幅 (eV)。
:returns: 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` を用いて台形則で数値積分します。
この関数は上の同名関数を上書きします。積分刻み幅 `h` は `(E1 - E0) / n` で再計算されます。
:param func: callable: 積分する関数。引数としてエネルギー (E) を取り、値を返します。
:param E0: float: 積分範囲の開始エネルギー (eV)。
:param E1: float: 積分範囲の終了エネルギー (eV)。
:param h: float: 積分刻み幅 (eV)。
:returns: 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: list[float]: 積分する関数値のリスト。
:param h: float: データ間の刻み幅。
:returns: 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]:
"""
概要: エネルギー範囲をDOSデータ配列のインデックス範囲に変換します。
詳細説明: グローバル変数 `Enorm_min` と `Enorm_step` を使用して、
与えられたエネルギー `E0` と `E1` に対応するDOSデータ配列の開始・終了インデックスを計算します。
:param E0: float: エネルギー範囲の開始値 (eV)。
:param E1: float: エネルギー範囲の終了値 (eV)。
:returns: tuple[int, int]: 開始インデックスと終了インデックス。
"""
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:
"""
概要: 指定されたエネルギー範囲における電子密度を計算します。
詳細説明: 温度 `T` とフェルミ準位 `EF` における電子のDOSfe関数を、
`E0` から `E1` までの範囲で数値積分することにより電子密度を計算します。
積分方法はグローバル変数 `integrator` ('interpolate'または'raw') によって選択されます。
:param T: float: 温度 (K)。
:param EF: float: フェルミ準位 (eV)。
:param E0: float: 積分範囲の開始エネルギー (eV)。
:param E1: float: 積分範囲の終了エネルギー (eV)。
:returns: float: 計算された電子密度 (cm^-3)。
"""
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:
"""
概要: 指定されたエネルギー範囲における正孔密度を計算します。
詳細説明: 温度 `T` とフェルミ準位 `EF` における正孔のDOSfh関数を、
`E0` から `E1` までの範囲で数値積分することにより正孔密度を計算します。
積分方法はグローバル変数 `integrator` ('interpolate'または'raw') によって選択されます。
:param T: float: 温度 (K)。
:param EF: float: フェルミ準位 (eV)。
:param E0: float: 積分範囲の開始エネルギー (eV)。
:param E1: float: 積分範囲の終了エネルギー (eV)。
:returns: float: 計算された正孔密度 (cm^-3)。
"""
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]:
"""
概要: DOSデータから価電子帯上端 (EV) と伝導帯下端 (EC) を探索します。
詳細説明: DOSデータ `DOS` と対応するエネルギー `E` を基に、
バンドギャップ中央 `Emidgap` からDOS閾値 `DOSth` を超える点を探索し、
EVとECを決定します。
:param E: list[float]: エネルギーのリスト (eV)。
:param DOS: list[float]: 状態密度のリスト。
:param Emidgap: float: バンドギャップの中央エネルギーの推定値 (eV)。
:param DOSth: float: DOSがこの値を超えたときにバンド端と判定する閾値。
:returns: tuple[float, float]: 価電子帯上端 (EV) と伝導帯下端 (EC) のエネルギー (eV)。
"""
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:
"""
概要: 全電荷中性条件からの偏差 ($\Delta$Q) のフェルミ準位 (EF) による数値微分を計算します (レガシー版)。
詳細説明: この関数は現在使用されていないようです。`tkDefects` クラスの `diff` メソッドが使用されています。
`dQ` 関数が未定義のため、この関数は実行されません。
:param h: float: 微小なEFの摂動量 (eV)。
:param T: float: 電子温度 (K)。
:param EF: float: フェルミ準位 (eV)。
:param ECmin: float: 伝導帯のエネルギー下限 (eV)。
:param ECmax: float: 伝導帯のエネルギー上限 (eV)。
:param EVmin: float: 価電子帯のエネルギー下限 (eV)。
:param EVmax: float: 価電子帯のエネルギー上限 (eV)。
:param N0: float: (現在未使用)。
:returns: float: $\Delta$QのEFによる数値微分 ($\partial Q / \partial E_F$)。
"""
# 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):
"""
概要: 数値計算ソルバーの進行状況を表示するためのコールバック関数。
詳細説明: ソルバーの反復回数、現在の解、関数値、およびステップサイズをコンソールに出力します。
主に`tkequation.Equation`クラスのデバッグや進捗表示に使用されます。
:param obj: object: `tkequation.Equation` オブジェクトのインスタンス。
`iter`, `x`, `f`, `xa`, `xb`, `dx` 属性を持つことが期待されます。
:returns: int: 常に1を返します（コールバックの成功を示す）。
"""
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, DOSCAR）を読み込み、関連情報を設定します。
詳細説明: 結晶構造、セル体積、OUTCARからのフェルミ準位 (EF0) とISPIN情報、
EIGENVALからのバンド端 (EV, EC, Eg) とHOMO/LUMO準位、
DOSCARからの状態密度 (DOS) データを読み込み、
それらをグローバル変数に設定し、DOSデータは`EV=0`を基準に正規化されます。
また、POSCARと欠陥データ (`dH_path`) のサイト数の一致をチェックします。
:param vasp: tkVASP: tkVASPオブジェクトのインスタンス。
:param cry: tkCrystal: 完全結晶のtkCrystalオブジェクトのインスタンス。
:param cry_defect: tkCrystal: 欠陥結晶のtkCrystalオブジェクトのインスタンス。
:param defects: tkDefects: tkDefectsオブジェクトのインスタンス。
:param POSCAR_path: str: 完全結晶のPOSCARファイルのパス。
:param OUTCAR_path: str: OUTCARファイルのパス。
:param EIGENVAL_path: str: EIGENVALファイルのパス。
:param DOSCAR_path: str: DOSCARファイルのパス。
:param fp: object | None: 出力先のファイルポインタ。Noneの場合、標準出力に出力。
:returns: tuple[dict, dict, dict, dict, float, float, float, float, float]:
- outcarinf: dict: OUTCARから読み込んだ情報。
- eigenvalinf: dict: EIGENVALから読み込んだ情報。
- bandedgeinf: dict: バンド端情報。
- dosinf: dict: DOS情報。
- EV: float: 価電子帯上端 (EV) (eV、EV=0基準)。
- EC: float: 伝導帯下端 (EC) (eV、EV=0基準)。
- Eg: float: バンドギャップ (eV)。
- EHOMO: float: HOMO準位 (eV、EV=0基準)。
- ELUMO: float: LUMO準位 (eV、EV=0基準)。
"""
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ファイルから読み込まれた欠陥データ (`dH_path`) の`N0` (基準サイト数) と比較します。
不一致がある場合、ユーザーに続行するか終了するかを問い合わせる警告メッセージを表示します。
:param cry: tkCrystal: 結晶構造情報を持つtkCrystalオブジェクト。
:param defects: tkDefects: 欠陥情報を持つtkDefectsオブジェクト。
:param fp: object | None: 出力先のファイルポインタ。Noneの場合、標準出力に出力。
:returns: bool: サイト数チェックを続行する場合True、終了する場合False。
"""
# 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: 欠陥情報を持つtkDefectsオブジェクト。
:returns: tuple[str, list[list[list]], list[list[list]]]:
- name: str: 最後の欠陥名。
- allpts_list: list[list[list]]: 各欠陥名の全遷移点のリストのリスト。
- minpts_list: list[list[list]]: 各欠陥名の最小形成エンタルピー曲線上の点のリストのリスト。
"""
#==============================================
# 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ファイルに保存され、グラフとして可視化されます。
特に、欠陥密度が凍結される温度 (Tdef) を考慮した計算を行います。
"""
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():
"""
概要: フェルミ準位 (EF) 依存性の計算とプロットを実行します。
詳細説明: 欠陥形成エネルギーのExcelファイルとVASP計算結果ファイルを読み込み、
様々なEFにおける電子・正孔密度、欠陥密度、および全電荷の中性条件からの偏差を計算します。
計算結果はExcelファイルに保存され、グラフとして可視化されます。
特に、欠陥密度が凍結される温度 (Tdef) を考慮した計算を行います。
"""
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: matplotlib.axes.Axes: スケールを再設定するAxesオブジェクト。
:param EFmin: float: EF軸の最小値 (eV)。
:param EFmax: float: EF軸の最大値 (eV)。
:param dHmin: float: dH軸の最小値 (eV)。
:param dHmax: float: dH軸の最大値 (eV)。
"""
# 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():
"""
概要: スクリプトのメイン実行関数。
詳細説明: アプリケーションの設定を読み込み、コマンドライン引数を解析してグローバル変数を更新します。
ログファイルを初期化し、`mode` (フェルミ準位依存性または温度依存性) に応じて
`exec_EF` または `exec_T` 関数を呼び出します。
不正なモードが指定された場合はエラーで終了します。
"""
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()