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.tkapplication import tkApplication from tklib.tkutils import save_csv 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.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 the correction of charged defects and dipole interaction by VASP """ #================================ # global parameters #================================ debug = 0 validmodes = ['Z'] mode = 'Z' CAR_IDIPOL_dir = '.' CAR_DIEL_el_dir = '.' CAR_DIEL_ion_dir = '.' CAR_charge_dir = '.' summaryfile = None # Energy to search edge energies EF0 = 0.1 # Max number of electrons occupies one state Nemax = 1.0 #=================================== # figure configuration #=================================== band_marker_size = 4 band_marker_edge_width = 0.5 Emin = -10.0 # eV Emax = 5.0 # eV figsize = (6, 8) #figsize = (6, 4) #colors = ['#000000', '#ff0000', '#00aa00', '#0000ff', '#aaaa00', '#ff00ff', '#00ffff', '#aa0000', '#00aa00', '#0000aa'] #colors = ['k', 'r', 'g', 'b', 'y', 'm', 'c'] fontsize = 16 labelfontsize = 12 legend_fontsize = 8 app = tkApplication() #============================= # Treat argments #============================= def usage(app = None): print("") print("Usage:") print(" (a) python {} CAR_IDIPOL_dir CAR_DIEL_el_dir CAR_DIEL_ion_dir CAR_charge_dir".format(sys.argv[0])) print(" (b) python {} CAR_IDIPOL_dir CAR_DIEL_el_dir epsilon CAR_charge_dir".format(sys.argv[0])) print(" ex: python {} {} {} {} {}".format(sys.argv[0], CAR_IDIPOL_dir, CAR_DIEL_el_dir, CAR_DIEL_ion_dir, CAR_charge_dir)) def updatevars(): global mode global CAR_IDIPOL_dir, CAR_DIEL_el_dir, CAR_DIEL_ion_dir, CAR_charge_dir # mode = getarg (1, mode) CAR_IDIPOL_dir = getarg(1, CAR_IDIPOL_dir) CAR_DIEL_el_dir = getarg(2, CAR_DIEL_el_dir) CAR_DIEL_ion_dir = getarg(3, CAR_DIEL_ion_dir) CAR_charge_dir = getarg(4, CAR_charge_dir) if mode == 'Z': pass else: app.terminate(f"Error: Invalide mode [{mode}]", usage = usage, pause = True) #def save_csv(path, headerlist, datalist, is_print = 0): def read_dipole_inf(outcarpath): fp = tkFile(outcarpath, 'r') if not fp or fp.fp is None: return None inf = {} line = fp.SkipTo('DIPCOR: dipole corrections for dipol') if line: # direction 1 min pos 1, direction 2 min pos 1, direction 3 min pos 7, fp.ReadLine() # dipolmoment 0.000000 -0.000000 -6.034717 electrons x Angstroem # Tr[quadrupol] -5.863427 line = fp.ReadLine() # print("l=", line) aa = line.split() inf["dipolemoment"] = [pfloat(aa[1]), pfloat(aa[2]), pfloat(aa[3])] aa = fp.ReadLine().split() inf["Tr[quadrupol]"] = pfloat(aa[1]) # fp.ReadLine() # energy correction for charged system 0.000000 eV line = fp.ReadLine() s = tkre.Search(r'([+\-]?[+\-\d\.eEdD]+)\s*eV', line) inf["Ecorr(charge)"] = pfloat(s[1]) # dipol+quadrupol energy correction -14.435503 eV line = fp.ReadLine() s = tkre.Search(r'([+\-]?[+\-\d\.eEdD]+)\s*eV', line) # print("s=", s) inf["Ecorr(di+quadrupol)"] = pfloat(s[1]) # added-field ion interaction 0.000000 eV (added to PSCEN) line = fp.ReadLine() s = tkre.Search(r'([+\-]?[+\-\d\.eEdD]+)\s*eV', line) inf["Ecorr(with field)"] = pfloat(s[1]) fp.Close() return inf def read_dielectric_inf_ion(outcarpath): fp = tkFile(outcarpath, 'r') if not fp or fp.fp is None: return None inf = {} while 1: eps1 = make_matrix2(3, 3) eps2 = make_matrix2(3, 3) line = fp.SkipTo(r"MACROSCOPIC STATIC DIELECTRIC TENSOR \(incl") if not line: break # ------------------------------------------------------ # 7.930930 0.000000 0.000000 # 0.000000 7.930923 0.000000 # -0.000000 0.000000 7.351176 # ------------------------------------------------------ line = fp.ReadLine() aa = make_matrix1(3) aa[0] = fp.ReadLine().split() aa[1] = fp.ReadLine().split() aa[2] = fp.ReadLine().split() for i in range(3): eps1[i][0] = pfloat(aa[i][0]) eps1[i][1] = pfloat(aa[i][1]) eps1[i][2] = pfloat(aa[i][2]) inf["epsilon(static/local field)"] = eps1 line = fp.SkipTo("MACROSCOPIC STATIC DIELECTRIC TENSOR IONIC") if not line: break line = fp.ReadLine() aa[0] = fp.ReadLine().split() aa[1] = fp.ReadLine().split() aa[2] = fp.ReadLine().split() for i in range(3): eps2[i][0] = pfloat(aa[i][0]) eps2[i][1] = pfloat(aa[i][1]) eps2[i][2] = pfloat(aa[i][2]) inf["epsilon(static/ion)"] = eps2 return inf def read_dielectric_inf_el(outcarpath): fp = tkFile(outcarpath, 'r') if not fp or fp.fp is None: return None inf = {} while 1: line = fp.SkipTo(r"REAL DIELECTRIC FUNCTION") if not line: break # E(ev) X Y Z XY YZ ZX # -------------------------------------------------------------------------------------------------------------- # 0.000000 2.861856 2.864431 2.893541 0.000807 0.000000 0.000000 # 0.058843 2.861930 2.864506 2.893617 0.000807 0.000000 0.000000 # 0.117686 2.862151 2.864728 2.893843 0.000808 0.000000 0.000000 fp.ReadLine() fp.ReadLine() line = fp.ReadLine() aa = line.split() eps_r = make_matrix2(3, 3) eps_r[0][0] = pfloat(aa[1]) # X eps_r[1][1] = pfloat(aa[2]) # Y eps_r[2][2] = pfloat(aa[3]) # Z eps_r[0][1] = pfloat(aa[4]) # XY eps_r[1][0] = pfloat(aa[4]) # XY eps_r[0][2] = pfloat(aa[6]) # ZX eps_r[2][0] = pfloat(aa[6]) # ZX eps_r[1][2] = pfloat(aa[5]) # YZ eps_r[2][1] = pfloat(aa[5]) # YZ inf["eps_opt_r"] = eps_r break return inf def correct_charged_defects(mode, CAR_IDIPOL_dir, CAR_DIEL_el_dir, CAR_DIEL_ion_dir, CAR_charge_dir): CAR_dir1 = CAR_IDIPOL_dir CAR_dir2 = CAR_DIEL_el_dir CAR_dir3 = CAR_DIEL_ion_dir CAR_dir4 = CAR_charge_dir vasp1 = tkVASP() vasp2 = tkVASP() vasp3 = tkVASP() vasp4 = tkVASP() base_path1 = vasp1.getdir(CAR_dir1) base_path2 = vasp2.getdir(CAR_dir2) base_path3 = vasp2.getdir(CAR_dir3) base_path4 = vasp2.getdir(CAR_dir4) print("") print(f"Summary file: {summaryfile}") INCAR_path1 = vasp1.get_INCAR(base_path1) POSCAR_path1 = vasp1.get_POSCAR(base_path1) POTCAR_path1 = vasp1.get_POTCAR(base_path1) OUTCAR_path1 = vasp1.get_OUTCAR(base_path1) print("Defect CAR dir(IDIPOL=4) : ", CAR_dir1) print(" INCAR : ", INCAR_path1) print(" POSCAR : ", POSCAR_path1) print(" POTCAR : ", POTCAR_path1) print(" OUTCAR : ", OUTCAR_path1) INCAR_path2 = vasp2.get_INCAR(base_path2) POSCAR_path2 = vasp2.get_POSCAR(base_path2) OUTCAR_path2 = vasp2.get_OUTCAR(base_path2) print("CAR dir(electron dielectric constant) : ", CAR_dir2) print(" INCAR : ", INCAR_path2) print(" POSCAR : ", POSCAR_path2) print(" OUTCAR : ", OUTCAR_path2) INCAR_path3 = vasp3.get_INCAR(base_path3) POSCAR_path3 = vasp3.get_POSCAR(base_path3) OUTCAR_path3 = vasp3.get_OUTCAR(base_path3) print("CAR dir(ionic dielectric constant) : ", CAR_dir3) print(" INCAR : ", INCAR_path3) print(" POSCAR : ", POSCAR_path3) print(" OUTCAR : ", OUTCAR_path3) INCAR_path4 = vasp4.get_INCAR(base_path4) POSCAR_path4 = vasp4.get_POSCAR(base_path4) POTCAR_path4 = vasp4.get_POTCAR(base_path4) OUTCAR_path4 = vasp4.get_OUTCAR(base_path4) print("CAR dir(defect charge) : ", CAR_dir4) print(" INCAR : ", INCAR_path4) print(" POSCAR : ", POSCAR_path4) print(" POTCAR : ", POTCAR_path4) print(" OUTCAR : ", OUTCAR_path4) print("") print("##########################################") print("# Defect model (IDIPOL=4)") print("##########################################") print("") print("*** Read INCAR from [{}]".format(INCAR_path1)) incarinf1 = vasp1.read_incar_inf(INCAR_path1) print("") print("*** Read POTCAR from [{}]".format(POTCAR_path1)) potcarinf1 = vasp1.read_potcar_inf(POTCAR_path1) nPPTypes = potcarinf1["nPPTypes"] print("Pseudo potentials:") print(" nPPTypes=", nPPTypes) for i in range(nPPTypes): ppname = potcarinf1["PP{}".format(i)] atomname = potcarinf1["Atom{}".format(i)] ne = potcarinf1["Ne{}".format(i)] print(" {:2d}: {}".format(i, ppname)) print(" {} ne={}".format(atomname, ne)) """ print("*** Read crystal structure from [{}]".format(POSCAR_path1)) cry1 = vasp1.read_poscar(POSCAR_path1) if cry1 is None: app.terminate(f"Error: Can not read [{POSCAR_path1}]", usage = usage, pause = True) a1, b1, c1, alpha1, beta1, gamm1 = cry1.LatticeParameters() Vcell = cry1.Volume() cry1.PrintInf("cell") print(" Vcell={} A^3".format(Vcell)) L = Vcell**(1.0/3.0) print(" L={:10.6g} A".format(L)) poscarinf1 = vasp1.poscarinf nAtomTypes = poscarinf1["nAtomTypes"] print(" nAtomTypes=", nAtomTypes) print(" AtomTypes =", poscarinf1["AtomTypes"]) print(" nAtomTypes=", poscarinf1["nAtomTypes"]) print(" nAtoms =", poscarinf1["nAtoms"]) """ print("") print("*** Read OUTCAR from [{}]".format(OUTCAR_path1)) outcarinf1 = vasp1.read_outcar_inf(OUTCAR_path1) ISPIN = outcarinf1["ISPIN"] NELECT = outcarinf1["NELECT"] IDIPOL = outcarinf1.get("IDIPOL", None) print("Information in OUTCAR:") print(" ISPIN : ", ISPIN) print(" IDIPOL: ", IDIPOL) print(" NELECT: ", NELECT) if IDIPOL is None or IDIPOL != 4: app.terminate(f"Error: IDIPOL in {CAR_dir1} must be 4", usage = usage, pause = True) print("") print("##########################################") print("# Ideal model (electronic dielectric constant)") print("##########################################") if IsFile(CAR_DIEL_el_dir) or IsDir(CAR_DIEL_el_dir): print("*** Read electronic dielectric constants from [{}]".format(OUTCAR_path2)) dielinf_el = read_dielectric_inf_el(OUTCAR_path2) print("") print("*** Read INCAR from [{}]".format(INCAR_path2)) incarinf2 = vasp2.read_incar_inf(INCAR_path2) print("*** Read OUTCAR from [{}]".format(OUTCAR_path2)) outcarinf2 = vasp3.read_outcar_inf(OUTCAR_path2) ISPIN = outcarinf2["ISPIN"] LOPTICS = outcarinf2.get("LOPTICS", None) print("Information in OUTCAR:") print(" ISPIN : ", ISPIN) print(" LOPTICS: ", LOPTICS) if LOPTICS is None or 't' not in LOPTICS.lower(): app.terminate(f"Error: LOPTICS in {CAR_dir2} must be .True.", usage = usage, pause = True) eps_el = dielinf_el["eps_opt_r"] print("e(el): |{:10.4g} {:10.4g} {:10.4g}|".format(eps_el[0][0], eps_el[0][1], eps_el[0][2])) print(" |{:10.4g} {:10.4g} {:10.4g}|".format(eps_el[1][0], eps_el[1][1], eps_el[1][2])) print(" |{:10.4g} {:10.4g} {:10.4g}|".format(eps_el[2][0], eps_el[2][1], eps_el[2][2])) eps_av_el = (eps_el[0][0] + eps_el[1][1] + eps_el[2][2]) / 3.0 print(" trace average: {:10.6g}".format(eps_av_el)) else: print("*** Read dielectric constant from argument [{}]".format(CAR_DIEL_el_dir)) eps_av_el = pfloat(CAR_DIEL_el_dir) print(" e(el): {:10.6g}".format(eps_av_el)) print("") print("##########################################") print("# Ideal model (ionic dielectric constant)") print("##########################################") if IsFile(CAR_DIEL_ion_dir) or IsDir(CAR_DIEL_ion_dir): print("*** Read ionic dielectric constants from [{}]".format(OUTCAR_path3)) dielinf = read_dielectric_inf_ion(OUTCAR_path3) print("") print("*** Read INCAR from [{}]".format(INCAR_path3)) incarinf3 = vasp3.read_incar_inf(INCAR_path3) print("*** Read OUTCAR from [{}]".format(OUTCAR_path3)) outcarinf3 = vasp3.read_outcar_inf(OUTCAR_path3) ISPIN = outcarinf3["ISPIN"] LEPSILON = outcarinf3.get("LEPSILON", None) print("Information in OUTCAR:") print(" ISPIN : ", ISPIN) print(" LEPSILON: ", LEPSILON) if LEPSILON is None or 't' not in LEPSILON.lower(): app.terminate(f"Error: LEPSILON in {CAR_dir3} must be .True.", usage = usage, pause = True) epstot = dielinf["epsilon(static/local field)"] epsion = dielinf["epsilon(static/ion)"] print("e(tot): |{:10.4g} {:10.4g} {:10.4g}|".format(epstot[0][0], epstot[0][1], epstot[0][2])) print(" |{:10.4g} {:10.4g} {:10.4g}|".format(epstot[1][0], epstot[1][1], epstot[1][2])) print(" |{:10.4g} {:10.4g} {:10.4g}|".format(epstot[2][0], epstot[2][1], epstot[2][2])) epstotav = (epstot[0][0] + epstot[1][1] + epstot[2][2]) / 3.0 print(" trace average: {:10.6g}".format(epstotav)) print("e(ion): |{:10.4g} {:10.4g} {:10.4g}|".format(epsion[0][0], epsion[0][1], epsion[0][2])) print(" |{:10.4g} {:10.4g} {:10.4g}|".format(epsion[1][0], epsion[1][1], epsion[1][2])) print(" |{:10.4g} {:10.4g} {:10.4g}|".format(epsion[2][0], epsion[2][1], epsion[2][2])) epsionav = (epsion[0][0] + epsion[1][1] + epsion[2][2]) / 3.0 print(" trace average: {:10.6g}".format(epsionav)) else: print("*** Read dielectric constant from argument [{}]".format(CAR_DIEL_ion_dir)) epsionav = pfloat(CAR_DIEL_ion_dir) print(" eps(tot): {:10.6g}".format(epsionav)) eps_tot = epsionav + eps_av_el print("") print("*** Read IDIPOL information from [{}]".format(OUTCAR_path1)) dipoleinf = read_dipole_inf(OUTCAR_path1) # print("d=", dipoleinf) print(" dipolemoment:", dipoleinf["dipolemoment"]) print(" Tr[quadrupol] :", dipoleinf["Tr[quadrupol]"]) print(" Ecorr(di+quadrupol): {:12.6g} eV".format(dipoleinf["Ecorr(di+quadrupol)"])) print(" Ecorr(with field) : {:12.6g} eV".format(dipoleinf["Ecorr(with field)"])) print(" Ecorr(charge)=q^2 * aM / L : {:12.6g} eV".format(dipoleinf["Ecorr(charge)"])) print("") print("##########################################") print("# Defect model (defect charge)") print("##########################################") print("") print("*** Read crystal structure from [{}]".format(POSCAR_path4)) cry4 = vasp4.read_poscar(POSCAR_path4) if cry4 is None: app.terminate(f"Error: Can not read [{POSCAR_path4}]", usage = usage, pause = True) a4, b4, c4, alpha4, beta4, gamm4 = cry4.LatticeParameters() Vcell = cry4.Volume() cry4.PrintInf("cell") print(" Vcell={} A^3".format(Vcell)) L = Vcell**(1.0/3.0) print(" L={:10.6g} A".format(L)) poscarinf4 = vasp4.poscarinf nAtomTypes = poscarinf4["nAtomTypes"] print(" nAtomTypes=", nAtomTypes) print(" AtomTypes =", poscarinf4["AtomTypes"]) print(" nAtomTypes=", poscarinf4["nAtomTypes"]) print(" nAtoms =", poscarinf4["nAtoms"]) print("") print("*** Read POTCAR from [{}]".format(POTCAR_path4)) potcarinf4 = vasp4.read_potcar_inf(POTCAR_path4) nPPTypes = potcarinf4["nPPTypes"] print("Pseudo potentials:") print(" nPPTypes=", nPPTypes) print("") print("*** Read OUTCAR from [{}]".format(OUTCAR_path4)) outcarinf4 = vasp4.read_outcar_inf(OUTCAR_path4) ISPIN = outcarinf4["ISPIN"] NELECT = outcarinf4["NELECT"] print("Information in OUTCAR:") print(" ISPIN : ", ISPIN) print(" NELECT: ", NELECT) nVe = 0.0 for i in range(nAtomTypes): key = "Ne{}".format(i) ne = potcarinf4[key] nAtoms = poscarinf4["nAtoms"][i] nVe += nAtoms * ne # print(" n=", i, ne, nAtoms, nVe) print(" # of total valcence electrons: ", nVe) q = nVe - NELECT print(" NELECT =", NELECT) print(" nVe =", nVe) print(" q =", q) print(" eps(tot)=", eps_tot, "e0") print(" L =", L, " A") print("") print("Summary:") print("Electronic dielectric constant (average): {:10.6g}".format(eps_av_el)) print("Ionic dielectric constant (average) : {:10.6g}".format(epsionav)) print("Total dielectric constant (average) : {:10.6g}".format(eps_tot)) if abs(q) <= 1.0e-3: aM = dipoleinf["Ecorr(charge)"] * L # print(" q^2 * aM : {:12.6g}".format(aM)) # dEcorr = 1.0 / 3.0 * q * q * dipoleinf["Ecorr(charge)"] / eps_tot dEcoor = 0.0 print(f" q={q} is very small, regarded as zero") print(" (1/3) * q^2 * aM / L / eps(total): {:12.6g} eV".format(dEcorr)) # dEcorr = 2.0 / 3.0 * q * q * dipoleinf["Ecorr(charge)"] / epstotav # print(" (2/3) * q^2 * aM / L / eps(tot,av): {:12.6g} eV".format(dEcorr)) else: aM = dipoleinf["Ecorr(charge)"] * L / q / q print(" aM : {:12.6g} eV*A".format(aM)) dEcorr = 1.0 / 3.0 * q * q * dipoleinf["Ecorr(charge)"] / eps_tot print(" (1/3) * q^2 * aM / L / eps(total): {:12.6g} eV".format(dEcorr)) # dEcorr = 2.0 / 3.0 * q * q * dipoleinf["Ecorr(charge)"] / epstotav # print(" (2/3) * q^2 * aM / L / eps(tot,av): {:12.6g} eV".format(dEcorr)) print("") print(f"Save summary to {summaryfile}") ini = tkIniFile() ini.write_from_scratch(summaryfile, 'results', Car_dir1 = base_path1, Car_dir2 = base_path2, Car_dir3 = base_path3, Car_dir4 = base_path4, eps_av_el = eps_av_el, epsionav = epsionav, eps_tot = eps_tot, q = q, L = L, Error_charge = dipoleinf["Ecorr(charge)"], dEcorr = dEcorr ) app.terminate("", usage = usage, pause = True) def main(): global mode global cifpath, poscarpath global summaryfile updatevars() vasp = tkVASP() base_path = vasp.getdir(CAR_IDIPOL_dir) # logfile = app.replace_path(infile) logfile = os.path.join(base_path, 'charge_correction-out.txt') summaryfile = os.path.join(base_path, 'charge_correction-summary.prm') print("") print(f"Open logfile [{logfile}]") app.redirect(targets = ["stdout", logfile], mode = 'w') print("") print("=============== Estimate the correction of charged defects and dipole interaction by VASP ============") print("") print("mode: ", mode) if mode == 'Z': correct_charged_defects(mode, CAR_IDIPOL_dir, CAR_DIEL_el_dir, CAR_DIEL_ion_dir, CAR_charge_dir) else: app.terminate(f"Error: Invalide mode [{mode}]", usage = usage, pause = True) if __name__ == "__main__": main()