import os import sys from math import exp, sqrt, log import numpy as np #from numpy import arange from scipy import optimize # newton関数はscipy.optimizeモジュールに入っている from scipy.interpolate import interp1d from matplotlib import pyplot as plt from tklib.tkfile import tkFile from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg, save_csv from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me from tklib.tktransport.tkTransport import integrate_Simpson_list, fe, fh, NC2meff_FEA, meff2DC0_FEA from tklib.tktransport.tkDOS_VASP import integrate_Simpson, integrate_Simpson_list, tkDOS """ Calculate carrier densities and Fermi level based on DOSCAR of VASP. """ #================================ # Global variables #================================ Debug = 0 #mode: 'me', 'EF', 'T' mode = 'me' # For mode = 'me'. E range for fitting to me*(DOS) dE_me_lsq = 0.5 # measured from EC/EV, eV # EF range for mode = 'EF' T0 = 300 # K dEFmin = -0.5 # measured from EV, eV dEFmax = 0.5 # measured from EC, eV nEF = 50 EFmin = None EFmax = None EFstep = None # Temperature range for mode = 'T' Tmin = 300 # K Tmax = 600 nT = 11 Tstep = None # read from DOSCAR dos = tkDOS() file = "DOSCAR" outcsvfile = None # Flag to specify the data to determine band edges #dos.data_for_bandedges = 'EIGENVAL' dos.data_for_bandedges = 'DOSCAR' dos.dEA = 0.05 # Acceptor level measured from EV, eV dos.NA = 0.8e17 dos.dED = -0.28 # Donor level measured from EC, eV dos.ND = 0.0e17 dos.EA = None # Acceptor states are not considered if EA = None dos.ED = None # Donor states are not considered if ED = None # other parameters dos.EF0 = 0.0 # EF at charge neutral dos.Egth = 1.0e18 # Critera DOS value to search band edges, cm^-3 # Interpolation type dos.interp_type = 'linear' #dos.interp_type = 'cubic' # integration parameters # Integration range w.r.t. kBT dos.Einteg0 = -3.0 dos.Einteg1 = 3.0 dos.nrange = 6.0 # Relative accuracy of the quad functin dos.eps_quad = 1.0e-5 dos.nmaxdiv_quad = 150 # Bisection method parameters # Initial search range dos.dE_bisec = 0.4 # EPS dos.eps_bisec = 1.0e-7 # max iteration number dos.nmaxiter_bisec = 200 # Output cycle dos.iprintiterval_bisec = 1 E_raw = None dos_raw = None nDOS = None Emin = None Emax = None Estep = None EV = None EC = None #============================= # Treat argments #============================= def usage(): global mode, dos global Vcell, T0 global file global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep global dE_me_lsq print("") print("Usage: Variables in () are optional") print(" python {} mode (file args)".format(sys.argv[0], )) print(" data_for_bandedges: Data to determine band edges [EIGENVAL|DOSCAR]") print(" (i) python {} me file dE data_for_bandedges".format(sys.argv[0])) print(" Plot DOS and fit for me*(DOS)") print(" dE: Energy range from EC/EV for fitting to m*(DOS)") print(" ex: python {} me {} {} {}".format(sys.argv[0], file, dE_me_lsq, dos.data_for_bandedges)) print(" (ii) python {} T file Tmin Tmax nT data_for_bandedges)".format(sys.argv[0])) print(" Calculate T dependences of ne, nh, ND+, NA-") print(" Tmin, Tmax, nT: T range and mesh") print(" ex: python {} T {} {} {} {} {}".format(sys.argv[0], file, Tmin, Tmax, nT, dos.data_for_bandedges)) print(" (iii) python {} EF file nEF T data_for_bandedges)".format(sys.argv[0])) print(" Calculate EF dependences of ne, nh, ND+, NA-") print(" nEF: # of EF mesh") print(" T: Temperature") print(" ex: python {} EF {} {} {} {}".format(sys.argv[0], file, nEF, T0, dos.data_for_bandedges)) def updatevars(): global mode, dos global Vcell, T0 global file, outcsvfile global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep global Tmin, Tmax, nT, Tstep global dE_me_lsq argv = sys.argv if len(argv) == 1: usage() exit() mode = getarg( 1, mode) file = getarg( 2, file) if mode == 'me': dE_me_lsq = getfloatarg(3, dE_me_lsq) dos.data_for_bandedges = getarg(4, dos.data_for_bandedges) elif mode == 'EF': nEF = getintarg(3, nEF) T0 = getfloatarg(4, T0) dos.data_for_bandedges = getarg(5, dos.data_for_bandedges) elif mode == 'T': Tmin = getintarg(3, Tmin) Tmax = getintarg(4, Tmax) nT = getintarg(5, nT) Tstep = (Tmax - Tmin) / (nT - 1) dos.data_for_bandedges = getarg(6, dos.data_for_bandedges) header, ext = os.path.splitext(file) filebody = os.path.basename(header) outcsvfile = filebody + f'-EFT-{mode}.csv' def exec_T(): global mode, dos, T0 print("") print("T dependence") print(" Temperature range: {} - {}, {} K step".format(Tmin, Tmax, Tstep)) xT = [Tmin + iT * Tstep for iT in range(nT)] xTinv, yEF, yNe, yNh, yNAm, yNDp, yRH, yNs, yNsabs = dos.cal_T_lists(xT, print_level = 0) print("") yEaNe = [] yEaNeTh = [] yEaNh = [] yEaNhTh = [] print ("%8s %10s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s" % ("T(K)", "EF(eV)", "Ea(Ne,meV)", "Ea(Ne,T1/2)", "Ea(Nh,meV)", "Ea(Nh,T1/2)", "Ne", "Nh", "ND+", "NA-", "RH", "Ns")) for i in range(nT): if i == 0: i1 = i i2 = i+1 elif i == nT - 1: i1 = i-1 i2 = i else: i1 = i-1 i2 = i+1 Th1 = sqrt(xT[i1]) Th2 = sqrt(xT[i2]) dT = 1.0 / (xTinv[i2] - xTinv[i1]) eps = 1.0e-300 EaNe = (log(yNe[i2]+eps) - log(yNe[i1])+eps) * dT # meV EaNeTh = (log(yNe[i2]/Th2+eps) - log(yNe[i1]/Th1)+eps) * dT # meV EaNh = (log(yNh[i2]+eps) - log(yNh[i1])+eps) * dT # meV EaNhTh = (log(yNh[i2]/Th2+eps) - log(yNh[i1]/Th1)+eps) * dT # meV yEaNe.append (-1000.0 * 1000.0 * kB / e * EaNe) yEaNeTh.append(-1000.0 * 1000.0 * kB / e * EaNeTh) yEaNh.append (-1000.0 * 1000.0 * kB / e * EaNh) yEaNhTh.append(-1000.0 * 1000.0 * kB / e * EaNhTh) print("%8.3f %10.4f %12.4g %12.4g %12.4g %12.4g " % (xT[i], yEF[i], yEaNe[i], yEaNeTh[i], yEaNh[i], yEaNhTh[i]) + "%12g %12g %12g %12g %12g %12g" % (yNe[i], yNh[i], yNDp[i], yNAm[i], yRH[i], yNs[i])) print("") print("Save fitting data to [{}]".format(outcsvfile)) save_csv(outcsvfile, ["T(K)", "EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH", "Ns", "Ea(Ne,meV)", "Ea(Ne,T1/2)", "Ea(Nh,meV)", "Ea(Nh,T1/2)"], [xT, yEF, yNe, yNh, yNDp, yNAm, yRH, yNs, yEaNe, yEaNeTh, yEaNh, yEaNhTh]) #============================= # Plot graphs #============================= fig = plt.figure(figsize = (12, 8)) ax1 = fig.add_subplot(2, 3, 1) ax2 = fig.add_subplot(2, 3, 2) ax3 = fig.add_subplot(2, 3, 4) ax4 = fig.add_subplot(2, 3, 3) ax5 = fig.add_subplot(2, 3, 6) ax1.plot(xT, yEF, label = '$E_F$', color = 'black') ax1.set_xlabel("T (K)") ax1.set_ylabel("$E_F$ (eV)") ax1.legend() ax2.plot(dos.Elist, dos.DOSlist, label = 'DOS', color = 'black') ax2.plot([dos.EF0, dos.EF0], [min(dos.DOSlist), max(dos.DOSlist)], color = 'red', linestyle = 'dashed') ax2.set_xlabel("E (eV)") ax2.set_ylabel("DOS (states/cm$^3$)") ax2.legend() ax3.plot(xT, yRH, label = 'RH', color = 'black', marker = 'o') ax3.set_xlabel("T (K)") ax3.set_ylabel("R$_H$ (C$^{-1}$cm$^{-3}$)") ax3.legend() ax4.plot(xTinv, yNe, label = '$N_e$', color = 'black', marker = 'o', linewidth = 0.5, markersize = 5.0) ax4.plot(xTinv, yNh, label = '$N_h$', color = 'red', marker = 'o', linewidth = 0.5, markersize = 5.0) ax4.plot(xTinv, yNsabs, label = '|$N_s$|', color = 'purple', linewidth = 0.5, linestyle = 'dashed', marker = '^', markersize = 2.0) ax4.plot(xTinv, yNDp, label = '$N_D^+$', color = 'blue', marker = 'x', linewidth = 0.5, markersize = 5.0) ax4.plot(xTinv, yNAm, label = '$N_A^-$', color = 'green', marker = 'x', linewidth = 0.5, markersize = 5.0) ax4.set_yscale('log') minNs = [min(yNe), min(yNh), min(yNsabs)] maxNs = [max(yNe), max(yNh), max(yNsabs)] ylim = [ max([min(minNs)*0.5, 1.0e8]), min([max(maxNs)*2.0, 1.0e23]) ] ax4.set_ylim(ylim) # ax4.plot([0, 0], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5) ax4.set_xlabel("1000/T (K$^{-1}$)") ax4.set_ylabel("N (cm$^{-3}$)") ax4.legend() ax5.plot(xTinv, yEaNe, label = '$E_a(N_e)$', linewidth = 0.5, color = 'black', marker = 'o', markersize = 5.0) ax5.plot(xTinv, yEaNeTh, label = '$E_a(N_e*T^{-0.5}$)', linewidth = 0.5, color = 'black', marker = 'x', markersize = 5.0) ax5.plot(xTinv, yEaNh, label = '$E_a(N_h)$', linewidth = 0.5, color = 'red', marker = 'o', markersize = 5.0) ax5.plot(xTinv, yEaNhTh, label = '$E_a(N_h*T^{-0.5}$)', linewidth = 0.5, color = 'red', marker = 'x', markersize = 5.0) ax5.set_xlabel("1000/T (K$^{-1}$)") ax5.set_ylabel("Ea (meV)") ax5.legend() # Rearange the graph axes so that they are not overlapped plt.tight_layout() plt.pause(0.1) print("") usage() print("") print("Press ENTER to exit>>", end = '') input() def exec_EF(): global mode, dos, T0 global EFmin, EFmax, EFstep print("") print("EF dependence at {} K".format(T0)) print(" EF range: {} - {}, {} eV step".format(EFmin, EFmax, EFstep)) xEF = [EFmin + iEF * EFstep for iEF in range(nEF)] yNh, yNe, yNAm, yNDp, yRH, yNs, yNsabs = dos.cal_EF_lists(T0, xEF, print_level = 0) yT0Ne = [] yT0Nh = [] print("") print ("%5s\t%8s\t%12s\t%12s\t%12s\t%12s\t%12s\t%12s\t%12s" % ("EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH", "T0(Ne)", "T0(Nh)", "Ns")) for i in range(nEF): EF = xEF[i] if i == 0: i1 = i i2 = i+1 elif i == nEF - 1: i1 = i-1 i2 = i else: i1 = i-1 i2 = i+1 dEF = xEF[i2] - xEF[i1] dlnNedEF = (log(yNe[i2]+1.0e-300) - log(yNe[i1]+1.0e-300)) / dEF dlnNhdEF = (log(yNh[i2]+1.0e-300) - log(yNh[i1]+1.0e-300)) / dEF if dlnNedEF == 0.0: dlnNedEF = 1.0-4 T0Ne = 1.0 / kB / dlnNedEF * e if dlnNhdEF == 0.0: dlnNhdEF = 1.0-4 T0Nh = -1.0 / kB / dlnNhdEF * e yT0Ne.append(T0Ne) yT0Nh.append(T0Nh) print ("%10.4f\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g" % (EF, yNe[i], yNh[i], yNDp[i], yNAm[i], yRH[i], T0Ne, T0Nh, yNs[i])) Emid = (dos.EV + dos.EC) / 2.0 imid = int((Emid - EFmin) / EFstep + 1.0e-5) Eimid = EFmin + imid * EFstep Ne = yNe[imid] Nh = yNh[imid] T0Ne = yT0Ne[imid] T0Nh = yT0Ne[imid] NC = Ne * exp((dos.EC - Eimid) * e / kB / T0Ne) NV = Nh * exp((Eimid - dos.EV) * e / kB / T0Nh) meeff = NC2meff_FEA(NC, T0) mheff = NC2meff_FEA(NV, T0) DC0 = meff2DC0_FEA(meeff, T0) DV0 = meff2DC0_FEA(mheff, T0) DOS_FEA = [] print("") print ("{:10}\t{:10}".format("E(eV)", "DOS,FEA(1/cm^3/eV)")) for i in range(nEF): EF = xEF[i] if EF < dos.EV: DOS_FEA.append(DV0 * sqrt(dos.EV - EF)) elif EF > dos.EC: DOS_FEA.append(DC0 * sqrt(EF - dos.EC)) else: DOS_FEA.append(0.0) print ("{:10.4g}\t{:10.4g}".format(xEF[i], DOS_FEA[i])) print("") print("Effective density of states at the mid gap {} eV (i={}, E={:8.4g}):".format(Emid, imid, Eimid)) print(" NC = {:12.6g} cm-3 (T0 = {:8.4g} K)".format(NC, T0Ne)) print(" NV = {:12.6g} cm-3 (T0 = {:8.4g} K)".format(NV, T0Nh)) print(" DC0 = {:12.6g} cm^-3/eV^1.5 (T0 = {:8.4g} K)".format(DC0, T0Ne)) print(" DV0 = {:12.6g} cm^-3/eV^1.5 (T0 = {:8.4g} K)".format(DV0, T0Nh)) print(" me* = {:8.4g} me".format(meeff)) print(" mh* = {:8.4g} me".format(mheff)) xEFapprox = [EFmin + i * EFstep for i in range(nEF)] yNeapprox = [NC * exp(-(dos.EC - xEFapprox[i]) * e / kB / T0) for i in range(nEF)] yNhapprox = [NV * exp(-(xEFapprox[i] - dos.EV) * e / kB / T0) for i in range(nEF)] print("") print("Save fitting data to [{}]".format(outcsvfile)) save_csv(outcsvfile, ["EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH", "T0(Ne)", "T0(Nh)", "Ns"], [xEF, yNe, yNh, yNDp, yNAm, yRH, yT0Ne, yT0Nh, yNs]) #============================= # Plot graphs #============================= fig = plt.figure(figsize = (8, 8)) ax2 = fig.add_subplot(2, 2, 1) ax3 = fig.add_subplot(2, 2, 2) ax4 = fig.add_subplot(2, 2, 3) ax1 = fig.add_subplot(2, 2, 4) ax2.plot(dos.Elist, dos.DOSlist, label = 'DOS', color = 'black') ax2.plot(xEF, DOS_FEA, label = 'DOS$_{,FEA}$', color = 'green', linewidth = 0.5, linestyle = 'dashed') ax2.plot([dos.EF0, dos.EF0], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5) ax2.plot([dos.EV, dos.EV], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5) ax2.plot([dos.EC, dos.EC], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5) ax2.set_xlabel("E (eV)") ax2.set_ylabel("DOS (states/cm$^3$)") ax2.legend() ax3.plot(xEF, yNe, label = '$N_e$', color = 'black', linewidth = 0.5) #, marker = 'o') ax3.plot(xEF, yNh, label = '$N_h$', color = 'red', linewidth = 0.5) #, marker = 'x') ax3.plot(xEF, yNsabs, label = '|$N_s$|', color = 'purple', linewidth = 0.5, linestyle = 'dashed') #, marker = '^') ax3.plot(xEF, yNDp, label = '$N_D^+$', color = 'blue') #, marker = 'x') ax3.plot(xEF, yNAm, label = '$N_A^-$', color = 'green') #, marker = 'x') ax3.plot(xEFapprox, yNeapprox, label = 'Ne(Boltz)', color = 'gray', linewidth = 0.5) #, marker = 'x') ax3.plot(xEFapprox, yNhapprox, label = 'Nh(Boltz)', color = 'purple', linewidth = 0.5) #, marker = 'x') ax3.plot([dos.EV, dos.EV], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5) ax3.plot([dos.EC, dos.EC], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5) ax3.set_yscale('log') ax3.set_xlabel("EF (eV)") ax3.set_ylabel("N (cm$^{-3}$)") ax3.legend() ax4.plot(xEF, yRH, label = 'R$_H$', color = 'black', linewidth = 0.5, marker = 'o') ax4.set_xlabel("EF (eV)") ax4.set_ylabel("R$_H$ (C$^{-1}$cm$^{-3}$)") ax4.legend() ax1.plot(xEF, yT0Ne, label = 'T$_0$(Ne)', color = 'black') #, marker = 'o') ax1.plot(xEF, yT0Nh, label = 'T$_0$(Nh)', color = 'red') #, marker = 'x') ax1.plot([dos.EV, dos.EV], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5) ax1.plot([dos.EC, dos.EC], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5) ax1.set_ylim([0.0, 2.0*T0]) ax1.set_xlabel("EF (eV)") ax1.set_ylabel("T$_0$ (K)") ax1.legend() # Rearange the graph axes so that they are not overlapped plt.tight_layout() plt.pause(0.1) print("") usage() print("") print("Press ENTER to exit>>", end = '') input() def exec_me(): global dos global file print("") print("Plot DOS and calculate m*") print(" dE_me_lsq: {:8.4g} eV".format(dE_me_lsq)) print(" EV: {:8.4g} eV dEFmin: {:8.4g} eV".format(dos.EV, dEFmin)) print(" EC: {:8.4g} eV dEFmax: {:8.4g} eV".format(dos.EC, dEFmax)) plotE0 = min([dos.EV + dEFmin, dos.EV - dE_me_lsq]) plotE1 = max([dos.EC + dEFmax, dos.EC + dE_me_lsq]) print("Plot E range: {:8.4g} - {:8.4g} eV".format(plotE0, plotE1)) ne_from_dos = integrate_Simpson_list(dos.Elist, dos.DOSlist) E_me = [] dos2 = [] Emid = (dos.EV + dos.EC) / 2.0 c = 0 for i in range(dos.nE): E = dos.Elist[i] # print("ne=", dos.Emin, E, dos.DOS(E), ne) if E < plotE0 or plotE1 < E: continue E_me.append(E) dos2.append(dos.DOSlist[i]**2) me_me = [] nme = len(E_me) for i in range(nme): if i == 0: i0 = 0 i1 = 1 elif i == nme - 1: i0 = i - 1 i1 = i else: i0 = i - 1 i1 = i + 1 DC0 = (dos2[i1] - dos2[i0]) / (E_me[i1] - E_me[i0]) if DC0 < 0.0: sign = -1.0 else: sign = 1.0 DC0 = sqrt(abs(DC0)) # DC0 = sqrt(2.0) / pi / pi * pow(me * meff, 1.5) / hbar / hbar / hbar * 1.0e-6 * pow(e, 1.5) meff = pow(DC0 / sqrt(2.0) * pi * pi * hbar * hbar * hbar / 1.0e-6 / pow(e, 1.5), 2.0/3.0) / me meff *= sign # print("m=", DC0, meff) me_me.append(meff) #============================= # Plot graphs #============================= fig = plt.figure(figsize = (8, 8)) ax1 = fig.add_subplot(2, 2, 1) ax2 = fig.add_subplot(2, 2, 2) ax3 = fig.add_subplot(2, 2, 3) ax4 = fig.add_subplot(2, 2, 4) ax1.plot(dos.Elist, dos.DOSlist, label = 'DOS', color = 'black', linewidth = 0.5) ax1.plot([dos.EF0, dos.EF0], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5) ax1.plot([dos.EV, dos.EV], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5) ax1.plot([dos.EC, dos.EC], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5) ax1.set_xlabel("E (eV)") ax1.set_ylabel("DOS (states/cm$^3$/eV)") ax1.legend() # print("len=", len(dos.Elist), len(ne_from_dos)) # print("ne=", ne_from_dos) ax2.plot(dos.Elist, ne_from_dos, label = 'N', color = 'black', linewidth = 0.5) ylim = [min(ne_from_dos), max(ne_from_dos)] ax2.plot([dos.EF0, dos.EF0], ylim, label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5) ax2.plot([dos.EV, dos.EV], ylim, label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5) ax2.plot([dos.EC, dos.EC], ylim, label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5) ax2.set_xlabel("E (eV)") ax2.set_ylabel("N (states/cm$^3$)") ax2.legend() ax3.plot(E_me, dos2, label = 'DOS$^2$', color = 'black', linewidth = 0.5) ylim = [min(dos2), max(dos2)] ax3.plot([dos.EF0, dos.EF0], ylim, label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5) ax3.plot([dos.EV, dos.EV], ylim, label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5) ax3.plot([dos.EC, dos.EC], ylim, label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5) # ax3.plot([dos.EV - dE_me_lsq, dos.EV - dE_me_lsq], ylim, label = '$E_{fit}$', color = 'green', linestyle = 'dashed', linewidth = 0.5) # ax3.plot([dos.EC + dE_me_lsq, dos.EC + dE_me_lsq], ylim, label = '$E_{fit}$', color = 'green', linestyle = 'dashed', linewidth = 0.5) ax3.set_xlabel("E (eV)") ax3.set_ylabel("DOS$^2$ (states$^2$/cm$^6$)") ax3.legend() ax4.plot(E_me, me_me, label = '$m^*$', color = 'black', linewidth = 0.5) ylim = [min(me_me), max(me_me)] ax4.plot([dos.EF0, dos.EF0], ylim, label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5) ax4.plot([dos.EV, dos.EV], ylim, label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5) ax4.plot([dos.EC, dos.EC], ylim, label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5) # ax4.plot([dos.EV - dE_me_lsq, dos.EV - dE_me_lsq], ylim, label = '$E_{fit}$', color = 'green', linestyle = 'dashed', linewidth = 0.5) # ax4.plot([dos.EC + dE_me_lsq, dos.EC + dE_me_lsq], ylim, label = '$E_{fit}$', color = 'green', linestyle = 'dashed', linewidth = 0.5) ax4.set_xlabel("E (eV)") ax4.set_ylabel("$m^*$") ax4.legend() # Rearange the graph axes so that they are not overlapped plt.tight_layout() plt.pause(0.1) print("") usage() print("") print("Press ENTER to exit>>", end = '') input() def main(): global mode, dos, T0 global Vcell, file, E_raw, dos_raw global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep updatevars() print("") print("===============================================================================") print(" Calculate T dependences of semiconductor statistics using VASP results") print("===============================================================================") print("mode: ", mode) dos.read_vaspfiles(file) EFmin = dos.EV + dEFmin EFmax = dos.EC + dEFmax EFstep = (EFmax - EFmin) / (nEF - 1) print("") print("Integration configuration") print(" Integration E range: {} - {} eV, or EF+-{}*kBT eV".format(dos.Einteg0, dos.Einteg1, dos.nrange)) print(" quad() parameters") print(" releps=", dos.eps_quad) print(" nmaxdiv=", dos.nmaxdiv_quad) print("") print("EF at 0 K = ", dos.EF0, " eV") Ne0K = dos.Ne(0.0, dos.EF0, -5.0, dos.EF0, dos.Estep) # Ne0K = Ne(0.0, EF0, Einteg0, EF0) print(" Ne at 0 K = ", Ne0K, " /u.c.") print("") if mode == 'me': exec_me() elif mode == 'T': exec_T() elif mode == 'EF': exec_EF() else: terminate("Error: Invalid mode [{}]".format(mode), usage) if __name__ == '__main__': main()