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 from tklib.tkexcel import tkExcel from tklib.tkgraphic.tkplotevent import tkPlotEvent from tklib.tkgui.tksimple_gui import get_window_from_plt, CustomDialog_with_config """ Calculate carrier densities and Fermi level based on DOSCAR of VASP. """ # 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^-1"; 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): 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): return self._site class tkDefects(): def __init__(self, **args): 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): names = {} for i in range(self.ndefects): d = self.defects[i] names[d.name] = 1 return list(names.keys()) def cal_dE(self, iPoint, idefect, T, EF): 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, EF, ECmin, ECmax, EVmin, EVmax): # 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 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, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax): # 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, Tdef, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax): 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, Tdef, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax): # 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, Tdef, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax): Z, ne, nh, Nds, NdsTdef, dGs, Qtot = self.calculate_densities(T, Tdef, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax) return Qtot def diff(self, h, T, Tdef, NdsTdef, EF, ECmin, ECmax, EVmin, EVmax): 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, Tdef, NdsTdef, ECmin, ECmax, EVmin, EVmax, callbackfunc = None, eps_bisec = 1.0e-1, nmaxiter_bisection = 300, initial = [-1.0, 4.0], eps_newton = 1.0e-6, nmaxiter_newton = 300, dump_newton = 0.3, h_newton = 1.0e-6, is_print = False): 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, groupdata): 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, prevEF, EFmin, EFmax, groupdata): 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, idx, EF, dH): 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, groupdata, EFmin, EFmax): 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 iorder, idxmin, dH1 = self.find_order(groupdata, ngroupdata-1, EFmax, 1.0e100) allpts.append([EFmax, dH1, q0, groupdata[idxmin]["charge"], 0, idxmin]) minpts.append([EFmax, dH1, q0, groupdata[idxmin]["charge"], 0, idxmin]) # if name == 'O_Ge1': # if name == 'V_O1': # print("allpts=", allpts) # print("minpts=", minpts) # exit() return allpts, minpts def get_transitionlevel_lists(self): 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): 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): 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, iPoint): 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): self.V = V return self.V def Volume(self): return self.V def add(self, defect): self.defects.append(defect) def get(self, idx): try: return self.defects[idx] except: print("") print("Warning in tkDefects.get: No defect for index=", idx) print("") return None def Print(self, PrintLabels = False, outfp = 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, V): 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 = None): global CAR_path if CAR_path == '': CAR_path = '.' print("") print("Usage: Variables in () are optional") print(" Input files: INCAR, POSCAR, POTCAR, OUTCAR, DOSCAR, EIGENVAL of perfect crystal model") print(" Excel (.xlsx) file of defect formation energies") print(" python {} mode (other args)".format(sys.argv[0], )) print(" NOTE: If T(electron) <= T(defect) defect densities are frozen at T(defect)") print(" If T(electron) > T(defect) defect densities are calculated at T(electron)") print(" (i) python {} EF [min|all] dH_path CAR_path CAR_path(defect) iPoint T(electron) T(defect) EF(defect) EFmin EFmax nEF".format(sys.argv[0])) print(" ex: python {} EF min {} {} {} {} {} {} {} {} {} {}" .format(sys.argv[0], dH_path, CAR_path, CAR_path_defect, iPoint, T0, Tdef, EFdef, EFmin, EFmax, nEF)) print(" ex: python {} EF all {} {} {} {} {} {} {} {} {} {}" .format(sys.argv[0], dH_path, CAR_path, CAR_path_defect, iPoint, T0, Tdef, EFdef, EFmin, EFmax, nEF)) print(" (ii) python {} T dH_path, CAR_path iPoint T(defect) Tmin Tmax nT".format(sys.argv[0])) print(" ex: python {} T {} {} {} {} {} {} {} {} {}" .format(sys.argv[0], dH_path, CAR_path, CAR_path_defect, iPoint, Tdef, EFdef, Tmin, Tmax, nT)) def updatevars(): 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, header, datalist): 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, xmin = None, xmax = None, delimiter = ','): 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): 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, T, EF): 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, T, EF): return 1.0 - fe(E, T, EF) # Define the function to be integrated def DOSfe(E, T, EF): return DOS(E) * fe(E, T, EF) def DOSfh(E, T, EF): return DOS(E) * fh(E, T, EF) def integrate_trapezoid(func, E0, E1, h): 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(func, E0, E1, h): 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, h): n = len(y) S = 0.5 * (y[0] + y[n-1]) + sum(y[1:n-1]) return h * S def convert_range2index(E0, E1): 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 # Calculate the number of electrons from E0 to E1 with the Fermi level EF at T def Ne(T, EF, E0, E1): 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, EF, E0, E1): 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, DOS, Emidgap, DOSth): 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, T, EF, ECmin, ECmax, EVmin, EVmax, N0): return (dQ(T, EF + h, ECmin, ECmax, EVmin, EVmax, N0) - dQ(T, EF - h, ECmin, ECmax, EVmin, EVmax, N0)) / 2.0 / h def callbackfunc(obj): 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, cry, cry_defect, defects, POSCAR_path, OUTCAR_path, EIGENVAL_path, DOSCAR_path, fp = outfp): 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, defects, fp = None): # 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): #============================================== # 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(): 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 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(" {:>16} (Nds(tot)={:12.4g} cm-3):".format(labelq, NdsTdef[labelq]), end = '', fp = outfp) mprint(" Nds={:12.4g} cm-3 dGs={:12.4g} eV".format(Nds[id], dGs[id]), 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(): 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 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, EFmin, EFmax, dHmin, dHmax): # ax.cla() # plot_dH_EF(ax) ax.set_xlim([EFmin, EFmax]) ax.set_ylim([dHmin, dHmax]) # plt.draw() plt.pause(1.0e-4) vars.scale_callback = rescale vars.dH_axis = axdH EFlim = axdH.get_xlim() dHlim = axdH.get_ylim() vars.EF_min = EFlim[0] vars.EF_max = EFlim[1] vars.dH_min = dHlim[0] vars.dH_max = dHlim[1] plot_event.register_annotation_event(fig, activate = False, print_level = 0) plot_event.register_move_text_event(fig, activate = False) plot_event.register_popup_menu_event() # plot_event.register_pick(fig) # callback = lambda event: plot_event.onclick(event)) plt.pause(1.0e-4) app.terminate("", pause = True) #usage = usage, pause = True) if outfp: outfp.close() def main(): app.config.plugin_dir = app.replace_path(None, template = ["{dirname}", "plugin/vasp_defect"]) appconfig_path, config = app.read_app_config(print_level = 1) updatevars() vasp = tkVASP() base_path = vasp.getdir(CAR_path) logfile = app.replace_path(dH_path, template = ["{dirname}", "{filebody}-defect-out.txt"]) # logfile = os.path.join(base_path, 'vasp_defect-out.txt') print("") print(f"Open logfile [{logfile}]") app.redirect(targets = ["stdout", logfile], mode = 'w') if mode == 'EF': exec_EF() elif mode == 'T': exec_T() else: app.terminate("Error: Invalid mode [{}]".format(mode), usage = usage, pause = True) if __name__ == '__main__': main()