import os import sys import shutil import glob import csv import re import numpy as np from numpy import exp, log, sin, cos, tan, arcsin, arccos, arctan, pi from scipy.interpolate import interp1d from pprint import pprint from matplotlib import pyplot as plt from tklib.tkfile import tkFile from tklib.tkinifile import tkIniFile import tklib.tkre as tkre from tklib.tkutils import save_csv from tklib.tkapplication import tkApplication from tklib.tkutils import IsDir, IsFile, SplitFilePath from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg from tklib.tksci.tksci import Reduce01, Round from tklib.tksci.tkconvolution import convolution, convolve_func from tklib.tksci.tkconvolution import convolve_xydata from tklib.tksci.tkmatrix import make_matrix1, make_matrix2, make_matrix3 from tklib.tkcrystal.tkcif import tkCIF, tkCIFData from tklib.tkcrystal.tkcrystal import tkCrystal from tklib.tkcrystal.tkatomtype import tkAtomType from tklib.tkcrystal.tkvasp import tkVASP from tklib.tkcrystal.tkbandstructure import plot_band from tklib.tksci import tkequation import tklib.tkcsv """ Estimate band filling correction for defect model """ #================================ # global parameters #================================ debug = 0 validmodes = ['BF'] mode = 'BF' CAR_dir1 = '.' CAR_dir2 = '.' summaryfile = None # dEVBM obtained by VBM correction dEVBM = '' # FWHM of Gaussian smearing function for DOS WG_DOS = 0.1 # eV # Energy to search edge energies EF0 = 0.1 # occupancy threthold to separate HOMO and LUMO occ_th = 0.5 # Critera DOS value to search band edges DOSth = 1.0e-5 # Max number of electrons occupies one state Nemax = 1.0 # Plot configuration band_marker_size = 4 band_marker_edge_width = 0.5 Emin = -10.0 # eV Emax = 10.0 # eV #colors = ['#000000', '#ff0000', '#00aa00', '#0000ff', '#aaaa00', '#ff00ff', '#00ffff', '#aa0000', '#00aa00', '#0000aa'] #colors = ['k', 'r', 'g', 'b', 'y', 'm', 'c'] figsize = (10, 8) fontsize = 16 titlefontsize = 12 labelfontsize = 12 legend_fontsize = 12 app = tkApplication() #============================= # Treat argments #============================= def usage(): global mode if mode not in validmodes: mode = validmodes[0] print("") print("Usage:") print(" (a) python {} mode CAR_dir(ideal) CAR_dir(defect) dEVBM WG Emin Emax EF0".format(sys.argv[0])) print(" mode: {}".format(validmodes)) print(" ex: python {} {} {} {} {} {} {} {} {}".format(sys.argv[0], mode, CAR_dir1, CAR_dir2, dEVBM, WG_DOS, Emin, Emax, EF0)) def updatevars(): global mode global CAR_dir1, CAR_dir2 global EF0, dEVBM global Emin, Emax, WG_DOS mode = getarg (1, mode) CAR_dir1 = getarg (2, CAR_dir1) CAR_dir2 = getarg (3, CAR_dir2) dEVBM = getarg (4, dEVBM) WG_DOS = getfloatarg(5, WG_DOS) Emin = getfloatarg(6, Emin) Emax = getfloatarg(7, Emax) EF0 = getfloatarg(8, EF0) if mode == 'BF': pass else: terminate("Error: Invalide mode [{}]".format(mode), usage = usage) #def save_csv(path, headerlist, datalist, is_print = 0): def BF_correction(mode, CAR_path1, CAR_path2): global dEVBM, EF0 debug = 0 vasp1 = tkVASP() vasp2 = tkVASP() base_path1 = vasp1.getdir(CAR_path1) base_path2 = vasp2.getdir(CAR_path2) print("") print(f"Summary file: {summaryfile}") INCAR_path1 = vasp1.get_INCAR(base_path1) POSCAR_path1 = vasp1.get_POSCAR(base_path1) OUTCAR_path1 = vasp1.get_OUTCAR(base_path1) DOSCAR_path1 = vasp1.get_VASPPath(base_path1, 'DOSCAR') EIGENVAL_path1 = vasp1.get_VASPPath(base_path1, 'EIGENVAL') INCAR_path2 = vasp2.get_INCAR(base_path2) POSCAR_path2 = vasp2.get_POSCAR(base_path2) OUTCAR_path2 = vasp2.get_OUTCAR(base_path2) DOSCAR_path2 = vasp2.get_VASPPath(base_path2, 'DOSCAR') EIGENVAL_path2 = vasp2.get_VASPPath(base_path2, 'EIGENVAL') print("CAR dir(ideal) : ", CAR_path1) print(" INCAR : ", INCAR_path1) print(" POSCAR : ", POSCAR_path1) print(" OUTCAR : ", OUTCAR_path1) print(" DOSCAR : ", DOSCAR_path1) print(" EIGENVAL: ", EIGENVAL_path1) print("CAR dir(defect): ", CAR_path2) print(" INCAR : ", INCAR_path2) print(" POSCAR : ", POSCAR_path2) print(" OUTCAR : ", OUTCAR_path2) print(" DOSCAR : ", DOSCAR_path2) print(" EIGENVAL: ", EIGENVAL_path2) print("") print("dEVBM: {} eV".format(dEVBM)) print("") print("DOS plot") print(f" FWHM of Gauss function for convolution: {WG_DOS} eV") # print("") # print("EF0: {} eV".format(EF0)) dEVBM = pfloat(dEVBM, '') if dEVBM == '': app.terminate("Error: dEVBM value must be given", pause = True) print("") print("Ideal crystal model:") print("Read crystal structure from [{}]".format(POSCAR_path1)) cry1 = vasp1.read_poscar(POSCAR_path1) if cry1 is None: terminate("Error: Can not read [{}]".format(POSCAR_path1), usage = usage) a1, b1, c1, alpha1, beta1, gamm1 = cry1.LatticeParameters() Vcell1 = cry1.Volume() types1 = cry1.AtomTypeList() ntypes1 = len(types1) sites1 = cry1.ExpandedAtomSiteList() nsites1 = len(sites1) atomnamelist1 = cry1.get_atom_name_list(mode = 'all', NameOnly = True) poslist1 = cry1.get_position_list(mode = 'all', IsReduce01 = True) cry1.PrintInf("cell") print("") print("Read OUTCAR from [{}]".format(OUTCAR_path1)) outcarinf1 = vasp1.read_outcar_inf(OUTCAR_path1) if outcarinf1 is None: terminate("Error: Can not read [{}]".format(OUTCAR_path1), usage = usage) EF1 = vasp1.outcarinf["EF"] ISPIN1 = outcarinf1["ISPIN"] # ETOT1 = outcarinf1["TOTEN"] ETOT1 = outcarinf1["TOTEN_zero_sigma"] print(f"ISPIN={ISPIN1}") print(f"ETOT(sigma->0)={ETOT1}") print("") print("Defect crystal model:") print("Read from [{}]".format(POSCAR_path2)) cry2 = vasp2.read_poscar(POSCAR_path2) if cry2 is None: terminate("Error: Can not read [{}]".format(POSCAR_path2), usage = usage) a2, b2, c2, alpha2, beta2, gamm2 = cry2.LatticeParameters() Vcell2 = cry2.Volume() types2 = cry2.AtomTypeList() ntypes2 = len(types2) sites2 = cry2.ExpandedAtomSiteList() nsites2 = len(sites2) atomnamelist2 = cry2.get_atom_name_list(mode = 'all', NameOnly = True) poslist2 = cry2.get_position_list(mode = 'all', IsReduce01 = True) cry2.PrintInf("cell") print("") print("Read OUTCAR from [{}]".format(OUTCAR_path2)) outcarinf2 = vasp2.read_outcar_inf(OUTCAR_path2) if outcarinf2 is None: terminate("Error: Can not read [{}]".format(OUTCAR_path2), usage = usage) EF2 = vasp2.outcarinf["EF"] ISPIN2 = outcarinf2["ISPIN"] # ETOT2 = outcarinf2["TOTEN"] ETOT2 = outcarinf2["TOTEN_zero_sigma"] print(f"ISPIN={ISPIN2}") print(f"ETOT(sigma->0)={ETOT2}") # Estimate supercell size print("") nx = int(a2 / a1 + 0.2) ny = int(b2 / b1 + 0.2) nz = int(c2 / c1 + 0.2) print("Supercell multipliers of the defect crystal with respect to the ideal crystal") print(" (nx, ny, nz) = ({}, {}, {})".format(nx, ny, nz)) print("") print(f"Ideal crystal model: Read DOSCAR from [{DOSCAR_path1}]") doscarinf1 = vasp1.read_doscar(DOSCAR_path1, EF = EF1, normalize_E = False, unit = '/cm3') nE1 = doscarinf1["nE"] E1 = doscarinf1["E"] if nE1 == 0: Edmin1 = 0.0 Edmax1 = 0.0 else: Edmin1 = min(E1) Edmax1 = max(E1) tDOSup1 = doscarinf1["TotalDOSup"] Neup1 = doscarinf1["Neup"] tDOSdn1 = doscarinf1["TotalDOSdn"] Nedn1 = doscarinf1["Nedn"] if nE1 > 0: print(f"EF={EF1} eV") print(f"DOS E range: {Edmin1:10.6g} - {Edmax1:10.6g} eV, {nE1} points") print(f" Energy is not normalized") else: print("") app.terminate(f"Error: Zero nE for DOS(E) in {DOSCAR_path1}", pause = True) print(f"Read EIGENVAL from [{DOSCAR_path1}]") eigenvalinf1 = vasp1.read_eigenval(EIGENVAL_path1, EF = EF1, normalize_E = False) nk1 = eigenvalinf1["nk"] nLevels1 = eigenvalinf1["nLevels"] bandedgeinf1a = vasp1.find_band_edges_from_eigenval(EF0 = EF0, eigenvalinf = eigenvalinf1, ISPIN = ISPIN1, occ_th = occ_th) bandedgeinf1b = vasp1.gbandedges(base_path1, OUTCAR_path1, EIGENVAL_path1, EF1) print("k points in EIGENVAL:") print("nk=", nk1) print("nLevels=", nLevels1) """ for i in range(nk1): el = eigenvalinf1['EList'][i] kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = el print(" ({:8.4f} {:8.4f} {:8.4f}) w={:8.4g} {:8.4f} {:12.4f}".format(kx, ky, kz, wk, dk, ktot)) """ EV1a = bandedgeinf1a["EV"] EC1a = bandedgeinf1a["EC"] Eg1a = bandedgeinf1a["Eg"] EV1 = bandedgeinf1b["EVBM"] EC1 = bandedgeinf1b["ECBM"] Eg1 = bandedgeinf1b["Eg"] EHOMO1 = bandedgeinf1a["EHOMO"] ELUMO1 = bandedgeinf1a["ELUMO"] print(f"Band edge from EIGENVAL:") print(f"EF0={EF0:10.6f} eV") print(f" find_band_edges: EV={EV1a:10.6f} EC={EC1a:10.6f} Eg={Eg1a:10.6f} eV") print(f" HOMO:{EHOMO1:10.6f} eV") print(f" LUMO:{ELUMO1:10.6f} eV") print(f" gbandedges() : EV={EV1:10.6f} EC={EC1:10.6f} Eg={Eg1:10.6f} eV") print("") print(f"Defect crystal model: Read DOSCAR from [{DOSCAR_path2}]") doscarinf2 = vasp2.read_doscar(DOSCAR_path2, EF = EF2, normalize_E = False, unit = '/cm3') nE2 = doscarinf2["nE"] E2 = doscarinf2["E"] if nE2 == 0: Edmin2 = 0.0 Edmax2 = 0.0 else: Edmin2 = min(E2) Edmax2 = max(E2) tDOSup2 = doscarinf2["TotalDOSup"] Neup2 = doscarinf2["Neup"] tDOSdn2 = doscarinf2["TotalDOSdn"] Nedn2 = doscarinf2["Nedn"] if nE2 > 0: print(f"EF={EF2} eV") print(f"DOS E range: {Edmin2:10.6g} - {Edmax2:10.6g} eV, {nE2} points") print(f" Energy is not normalized") print(f"Read EIGENVAL from [{DOSCAR_path2}]") eigenvalinf2 = vasp2.read_eigenval(EIGENVAL_path2, EF = EF2, normalize_E = False) nk2 = eigenvalinf2["nk"] nLevels2 = eigenvalinf2["nLevels"] bandedgeinf2a = vasp2.find_band_edges_from_eigenval(EF0 = EF0, eigenvalinf = eigenvalinf2, ISPIN = ISPIN2, occ_th = occ_th) bandedgeinf2b = vasp2.gbandedges(base_path2, OUTCAR_path2, EIGENVAL_path2, EF2) print("k points in EIGENVAL:") print("nk=", nk2) print("nLevels=", nLevels2) for i in range(nk2): el = eigenvalinf2['EList'][i] kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = el print(" ({:8.4f} {:8.4f} {:8.4f}) w={:8.4g} {:8.4f} {:12.4f}".format(kx, ky, kz, wk, dk, ktot)) EV2a = bandedgeinf1a["EV"] EC2a = bandedgeinf1a["EC"] Eg2a = bandedgeinf1a["Eg"] EV2 = bandedgeinf1b["EVBM"] EC2 = bandedgeinf1b["ECBM"] Eg2 = bandedgeinf1b["Eg"] EHOMO2 = bandedgeinf1a["EHOMO"] ELUMO2 = bandedgeinf1a["ELUMO"] print(f"Band edge from EIGENVAL:") print(f"EF0={EF0:10.6f} eV") print(f" find_band_edges: EV={EV2a:10.6f} EC={EC2a:10.6f} Eg={Eg2a:10.6f} eV") print(f" HOMO:{EHOMO2:10.6f} eV") print(f" LUMO:{ELUMO2:10.6f} eV") print(f" gbandedges() : EV={EV2:10.6f} EC={EC2:10.6f} Eg={Eg2:10.6f} eV") print("") print("DOS plot range") print(f" Ideal crystal model:") print(f" EF={EF1} eV") print(f" E range: {Edmin1:10.6g} - {Edmax1:10.6g} eV, {nE1} points") """ print(f" DOS(E) is scaled for the size of the defect model by {nx}x{ny}x{nz}") k = nx * ny * nz * 1.0 for i in range(len(E1)): tDOSup1[i] *= k Neup1[i] *= k if ISPIN1 == 2: tDOSdn2[i] *= k Nedn1[i] *= k """ print(f" Defect crystal model:") print(f" Original EF={EF2} eV") print(f" Original E range: {Edmin2:10.6g} - {Edmax2:10.6g} eV, {nE2} points") EF2 -= dEVBM Edmin2 -= dEVBM Edmax2 -= dEVBM for i in range(len(E2)): E2[i] -= dEVBM for ik in range(nk2): el = eigenvalinf2['EList'][ik] kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = [*el] Eups[ik] -= dEVBM if ISPIN2 == 2: Edns[ik] -= dEVBM print(f" Adjusted EF={EF2} eV") print(f" Adjusted E range: {Edmin2:10.6g} - {Edmax2:10.6g} eV, {nE2} points") print("") print("Calculate band filling energy") totalw = 0.0 totalNe = 0.0 dEe = 0.0 dEh = 0.0 for ik in range(nk2): el = eigenvalinf2['EList'][ik] kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = [*el] totalw += wk print(" k[{:3d}]: ({:8.4f} {:8.4f} {:8.4f}) w={:8.4g}".format(ik, kx, ky, kz, wk)) for iE in range(nLevels2): Eup = Eups[iE] neup = occups[iE] totalNe += neup * wk if Eup <= EV1 and abs(neup - Nemax) > 1.0e-6: dEh += wk * (Nemax - neup) * (Eup - EV1) print(f" dEh: up iE={iE:3d} E={Eup:10.6g} eV dNh={Nemax - neup:10.4g} dEh={dEh:12.6g} eV") if Eup >= EC1 and neup > 1.0e-6: dEe += wk * neup * (Eup - EC1) print(f" dEe: up iE={iE:3d} E={Eup:10.6g} eV dNe={neup:10.4g} dEe={dEe:12.6g} eV") if ISPIN2 == 2: Edn = Edns[iE] nedn = occdns[iE] totalNe += nedn * wk if Edn <= EV1 and abs(nedn - Nemax) > 1.0e-6: dEh += wk * (Nemax - nedn) * (Edn - EV1) print(f" dEh: dn iE={iE:3d} E={Edn:10.6g} eV dNh={Nemax - nedn:10.4g} dEh={dEh:12.6g} eV") if Edn >= EC1 and nedn > 1.0e-6: dEe += wk * nedn * (Edn - EC1) print(f" dEe: dn iE={iE:3d} E={Edn:10.6g} eV dNe={nedn:10.4g} dEe={dEe:12.6g} eV") print("") print("From defect crystal model:") print(" Total weight: ", totalw) print(" NELECT : ", outcarinf2["NELECT"]) dEtot = dEh + dEe if ISPIN2 == 1: totalNe *= 2.0 dEh *= 2.0 dEe *= 2.0 dEtot *= 2.0 print(" Count spin degeneracy of 2:") print(" Total Ne : ", totalNe) print(" dEh = {:12.6g} eV".format(dEh)) print(" dEe = {:12.6g} eV".format(dEe)) print(" dEtot = {:12.6g} eV".format(dEtot)) print("") print(f"Save summary to {summaryfile}") ini = tkIniFile() ini.write_from_scratch(summaryfile, 'results', Car_dir1 = base_path1, Car_dir2 = base_path2, Etot1 = ETOT1 * nx * ny * nz, Etot2 = ETOT2, totalNe = totalNe, dEh = dEh, dEe = dEe, dEtot = dEtot ) #============================= # Plot graphs #============================= print("") fig = plt.figure(figsize = figsize) ax1 = fig.add_subplot(2, 2, 1) ax2 = fig.add_subplot(2, 2, 3) axband = fig.add_subplot(1, 2, 2) print("E1=", E1) tDOSup_s1 = convolution(E1, tDOSup1, WG_DOS, func_type = 'gauss') ax1.plot(E1, tDOSup_s1, label = 'up(ideal)', linestyle = '-', linewidth = 0.5, color = 'black') if ISPIN1 == 2: tDOSdn_s1 = convolution(E, tDOSdn1, WG_DOS, func_type = 'gauss') ax1.plot(E1, tDOSdn_s1, label = 'dn(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'black') ylim = ax1.get_ylim() ax1.plot([EV1, EV1], ylim, label = '$E_V$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'red') ax1.plot([EC1, EC1], ylim, label = '$E_C$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'red') ax1.plot([EF1, EF1], ylim, label = '$E_F$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'blue') # ax1.plot([EF0, EF0], ylim, label = '$E_{F0}$', linestyle = 'dashed', linewidth = 0.5, color = 'blue') # ax1.plot([EHOMO1, EHOMO1], ylim, label = '$E_{HOMO}$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'green') # ax1.plot([ELUMO1, ELUMO1], ylim, label = '$E_{LUMO}$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'green') # ax1.set_xlabel("$E$ (eV)", fontsize = fontsize) ax1.set_ylabel("DOS (states/cm$^3$/eV)", fontsize = fontsize) # ax1.set_xlim([Emin, Emax]) ax1.tick_params(labelsize = fontsize) tDOSup_s2 = convolution(E2, tDOSup2, WG_DOS, func_type = 'gauss') ax1.plot(E2, tDOSup_s2, label = 'up(defect)', linestyle = '-', linewidth = 1.0, color = 'red') if ISPIN2 == 2: tDOSdn_s2 = convolution(E, tDOSdn2, WG_DOS, func_type = 'gauss') ax1.plot(E2, tDOSdn_s2, label = 'dn(defect)', linestyle = 'dashed', linewidth = 1.0, color = 'red') ax1.plot([EF2, EF2], ylim, label = '$E_F$(defect)', linestyle = 'dashed', linewidth = 1.0, color = 'blue') ax1.tick_params(labelsize = fontsize) ax1.set_title(f"{CAR_path1}\n{CAR_path2}", fontsize = titlefontsize) ax1.legend(fontsize = legend_fontsize) ax2.plot(E2, Neup2, label = 'up(defect)', linestyle = '-', linewidth = 0.5, color = 'black') if ISPIN2 == 2: ax2.plot(E2, Nedn2, label = 'dn(defect)', linestyle = '-', linewidth = 0.5, color = 'red') ylim = ax2.get_ylim() ax2.plot([EV1, EV1], ylim, label = '$E_V$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'red') ax2.plot([EC1, EC1], ylim, label = '$E_C$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'red') # ax2.plot([EF0, EF0], ylim, label = '$E_{F0}$', linestyle = 'dashed', linewidth = 0.5, color = 'blue') # ax2.plot([EHOMO1, EHOMO1], ylim, label = '$E_{HOMO}$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'green') # ax2.plot([ELUMO1, ELUMO1], ylim, label = '$E_{LUMO}$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'green') ax2.plot([EF1, EF1], ylim, label = '$E_F$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'blue') ax2.plot([EF2, EF2], ylim, label = '$E_F$(defect)', linestyle = 'dashed', linewidth = 1.0, color = 'blue') ax2.set_xlabel("$E$ (eV)", fontsize = fontsize) ax2.set_ylabel("$N_e$ (states/unit cell)", fontsize = fontsize) # ax2.set_xlim([Emin, Emax]) ax2.legend(fontsize = legend_fontsize) ax2.tick_params(labelsize = fontsize) xk = [] yEup = [] yNup = [] yEdn = [] yNdn = [] for i in range(nk2): el = eigenvalinf2['EList'][i] kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = [*el] xk.append(ktot) yEup.append(Eups) yNup.append(occups) if ISPIN2 == 2: yEdn.append(Edns) yNdn.append(occdns) plot_band(plt, axband, ISPIN2, xk, yEup, yEdn, [Emin, Emax], occups = yNup, occdns = yNdn, ktotallist = None, ktotal_namelist = None, EV = EV1, EC = EC1, EF = EF2, EHOMO = None, ELUMO = None, EFlabel = '$E_{F}$(defect)', markersize = band_marker_size, markeredgewidth = band_marker_edge_width) axband.set_title(f"dEh={dEh:12.6g} eV\ndEe={dEe:12.6g} eV\ndEtot={dEtot:12.6g} eV", fontsize = titlefontsize) plt.tight_layout() # Rearange the graph axes so that they are not overlapped plt.tight_layout() """ plt.show() print("") print("Close graph window to terminate") """ plt.pause(0.1) app.terminate("", pause = True) def main(): global mode global cifpath, poscarpath global summaryfile updatevars() vasp = tkVASP() base_path = vasp.getdir(CAR_dir2) logfile = os.path.join(base_path, 'BF_correction-out.txt') summaryfile = os.path.join(base_path, 'BF_correction-summary.prm') print("") print(f"Open logfile [{logfile}]") app.redirect(targets = ["stdout", logfile], mode = 'w') print("") print("=============== Estimate band filling correction for defect model calculated by VASP ============") print("") print("mode: ", mode) if mode == 'BF': BF_correction(mode, CAR_dir1, CAR_dir2) else: terminate(f"Error: Invalide mode [{mode}]", usage = usage) if __name__ == "__main__": main()