import os import sys from math import exp, sqrt import numpy as np from math import log, exp 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.tkapplication import tkApplication 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 """ 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 # exponent threshold nexp = 100.0 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") 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_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): d = self.defects[id] ZS[d.site] = 1.0 for id in range(self.ndefects): d = self.defects[id] dG = self.cal_dG(iPoint, id, Tdef, EF) # dH = d.dH0s[iPoint] + d.charge * EF Kdef = dG / kBTedef if Kdef > nexp: Kdef = nexp elif Kdef <= -nexp: Kdef = -nexp ZS[d.site] += exp(-Kdef) # Calculate defect densities Nds = [] dGs = [] Qtot = 0.0 for id in range(self.ndefects): d = self.defects[id] dG = self.cal_dG(iPoint, id, Tdef, EF) # dH = d.dH0s[iPoint] + d.charge * 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) 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) try: NdsTdef[label] += Nds[id] except: NdsTdef[label] = Nds[id] 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 ZSA = {} 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) dG = self.cal_dG(iPoint, id, Te, EF) # dH = d.dH0s[iPoint] + d.charge * 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) dGs = make_matrix1(self.ndefects) Qtot = 0.0 for id in range(self.ndefects): d = self.defects[id] label = "{}_{}".format(d.atom, d.site) dG = self.cal_dG(iPoint, id, Te, EF) # dH = d.dH0s[iPoint] + d.charge * EF # print("k=", kBTe, dH, d.N0, d.N0 * exp(-dH / kBTedef), self.V) Ke = dG / kBTe if Ke > nexp: Ke = nexp elif Ke <= -nexp: Ke = -nexp n = NdsTdef[label] * exp(-Ke) / ZSA[label] # cm^-3 Nds[id] = n dGs[id] = dG Qtot += d.charge * n 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): 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): 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): ndata = len(groupdata) dH2 = [] for id2 in range(ndata): d2 = groupdata[id2] q2 = d2["charge"] H2 = d2["dH0"] dH2.append(H2 + q2 * EF) # 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=", dH2) iorder = 0 for id2 in range(ndata): if dH2[id2] < dH: print(f" current: idx={idx} ({EF}, {dH}): for id2={id2}: dH2[id2]={dH2[id2]} < dH={dH}") iorder += 1 idxmin = idx dHmin = dH for id2 in range(ndata): if dH2[id2] < dHmin: dHmin = dH2[id2] idxmin = id2 print(f" at ({EF}, {dH}): iorder={iorder:2} for idx={idx} idxmin={idxmin:2} dHmin={dHmin:8.3g}") # iorder: (EF, dH)が同じグループで何番目か # idxmin, dHmin: 同じグループでEFにおいて最小のdHをもつidx,dH return iorder, idxmin, dHmin def find_all_transitions(self, groupdata, EFmin, EFmax): eps = 1.0e-6 ngroupdata = len(groupdata) allpts = [] minpts = [] idxmin = 0 d0 = groupdata[idxmin] q0 = d0["charge"] H0 = d0["dH0"] print() print(f"name[0]: {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: ({EFmin}, {dHmin}) q={q0} idx=0") prevEF = EFmin while 1: nextidx, nextEF = self.find_next_transition_EF(idxmin, prevEF, EFmin, EFmax, groupdata) if nextidx is None: break 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: iorder={iorder} ({nextEF}, {dHmin}) idx={nextidx} next_q={q1}") allpts.append([nextEF, dH1, q0, q1, iorder, idxmin]) if iorder == 0: minpts.append([nextEF, dH1, q0, q1, iorder, idxmin]) prevEF = nextEF 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]) 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(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(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} {:>5} {:^6} {:^8} {:^12} ".format("Defect", "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(" {:6}: {:5} {:6.2g} {:8.4g} {:12.4g} ".format(d.name, d.charge, d.entropy, d.N0, d.Ndoped), 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_defct, iPoint, T0, Tdef, EFdef, EFmin, EFmax, nEF)) print(" ex: python {} EF all {} {} {} {} {} {} {} {} {} {}" .format(sys.argv[0], dH_path, CAR_path, CAR_path_defct, 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_defct, 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 global Tmin, Tmax, nT, Tstep global iPoint 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 = getfloatarg(10, EFmin) EFmax = getfloatarg(11, EFmax) nEF = getintarg (12, nEF) 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) 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, 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, 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 = '' for d in defects.defects: if d.site not in checked.keys(): mprint(" {:4} {:4f}".format(d.site, d.N0)) checked[d.site] = 1 try: ns = nsites[d.site] except: ns = 0 if abs(d.N0 - ns) > 1.0e-3: is_error += "{} ({} != {}), ".format(d.site, d.N0, ns) if is_error != '': mprint("", fp = outfp) mprint("========================================================================", fp = outfp) mprint("Error: Site numbers mismatch between POSCAR and {}".format(dH_path), fp = outfp) mprint(" {}".format(is_error)) mprint("========================================================================", fp = outfp) print("") print("Choose continue: Enter continue to disregard this error") print("Choose stop : Enter stop to terminate this run, and correct site numbers in {} 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, 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(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) 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 = "{}_{}".format(d.atom, d.site) labelq = "{}^{}".format(label, d.charge) mprint(" {:>16} (Nds(tot)={:12.4g} cm-3):".format(labelq, NdsTdef[label]), 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$)") axdos.legend() 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]]) h1, l1 = axEF.get_legend_handles_labels() h2, l2 = axdH.get_legend_handles_labels() axEF.legend(h1 + h2, l1 + l2, bbox_to_anchor=(1.05, 1.0), loc='upper left', borderaxespad = 0, fontsize = legend_fontsize) 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() ax.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left', borderaxespad = 0, fontsize = legend_fontsize) # 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, defects, POSCAR_path, OUTCAR_path, EIGENVAL_path, DOSCAR_path, fp = outfp) # EF plot range EFmin = EV + dEFmin 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 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) 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(" {:>16} (Nds(tot)={:12.4g} cm-3):".format(labelq, NdsTdef[label]), end = '', fp = outfp) mprint(" Nds={:12.4g} cm-3 dGs={:12.4g} eV".format(Nds[id], dGs[id]), 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}") #============================= # 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) axdos.set_title("{}, for Point {}".format(dH_path, defects.plabels[iPoint])) axdH.set_title('$E_{F,eq}$: ' + "{:8.3g} eV(Tdef={} K)".format(EFeqTdef, Tdef) + " {:8.3g} eV(T0={} K)".format(EFeqT0, T0)) 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$)") axdos.legend() # 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]) axdQ.legend() 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]) names = defects.get_names() for iname in range(len(names)): color = colors[iname % len(colors)] args = {'color': color, 'markersize': 3.0, 'markeredgewidth': 0, 'markerfacecolor': color } name = names[iname] 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(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]) largs = { 'label': "{}".format(name) } axdH.plot(min_x, min_y, picker = True, linestyle = '-', marker = 'o', **args, **largs) for id in range(0, ngroupdata): d = groupdata[id] name = d["name"] dH0 = d["dH0"] q = d["charge"] if plotall: axdH.plot([EFmin, EFmax], [dH0 + q * EFmin, dH0 + q * EFmax], linestyle = 'dashed', linewidth = 0.5, **args) for ipt in range(len(minpts) - 1): axdH.plot([minpts[ipt][0], minpts[ipt+1][0]], [minpts[ipt][1], minpts[ipt+1][1]], linestyle = '-', linewidth = 1.0, **args) allwb.Close() minwb.Close() yrange = axdH.get_ylim() yrange = [yrange[0], min(yrange[1], view_dHmax)] # yrange = [yrange[0], view_dHmax] # axdH.plot(axdH.get_xlim(), [0.0, 0.0], color = 'red', linestyle = 'dashed', linewidth = 0.5) 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)") # axdH.legend() axdH.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left', borderaxespad = 0, fontsize = legend_fontsize) axdN.plot(xEF, yne, label = '$N_e$', picker = True, color = 'red', linestyle = '-', linewidth = 1.0) axdN.plot(xEF, ynh, label = '$N_h$', picker = True, 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 axdN.plot(xEF, ynds[id], label = "{}({})".format(d.name, d.charge), picker = True, linestyle = 'dashed', linewidth = 1.0) 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() axdN.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left', borderaxespad = 0, fontsize = legend_fontsize) # Rearange the graph axes so that they are not overlapped plt.tight_layout() plot_event.register_pick(fig) # callback = lambda event: plot_event.onclick(event)) """ plt.show() print("") print("Close graph window to terminate") """ plt.pause(0.1) app.terminate("", usage = usage, pause = True) if outfp: outfp.close() def main(): 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()