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, getarg, getintarg, getfloatarg, save_csv #from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg from tklib.tkapplication import tkApplication from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me from tklib.tktransport.tkTransport import integrate_Simpson_list, fe, fh, NC2meff_FEA from tklib.tktransport.tkDOS_FEA import integrate_Simpson, integrate_Simpson_list, tkDOS """ Simulate N-T Hall data with Free Electron Approximation """ def pint(str): return int(str) def pfloat(str): return float(str) #================================ # Global variables #================================ app = tkApplication() argv, narg = app.get_args() params = app.get_param_dic() #params.debug = 0 #mode: 'me', 'EF', 'T' mode = 'me' outcsvfile = None # read from DOSCAR dos = tkDOS() dos.meeff = 0.3 dos.mheff = 1.0 dos.EV = 0.0 dos.EC = 1.1 dos.EA = 0.05 # Acceptor level, eV dos.NA = 0.0e17 dos.ED = 1.05 # Donor level dos.ND = 0.8e17 dos.EF0 = 0.0 # initial EF to find EF T0 = 300.0 # 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 dos.nEF = 51 dos.Estep = 0.01 # Integration step, eV # Temperature range for mode = 'T' Tmin = 300 # K Tmax = 600 nT = 11 Tstep = None # Integration range w.r.t. kBT dos.Einteg0 = -3.0 dos.Einteg1 = 3.0 dos.nrange = 6.0 # 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 figsize = (6, 6) #============================= # Treat argments #============================= def usage(): global mode, dos global T0 global dEFmin, dEFmax global dE_me_lsq print("") print("Usage: Variables in () are optional") print(" python {} mode (file args)".format(sys.argv[0], )) print(" (i) python {} me dE".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], dE_me_lsq)) print(" (ii) python {} T EC EA NA ED ND Tmin Tmax nT)".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], dos.EC, dos.EA, dos.NA, dos.ED, dos.ND, Tmin, Tmax, nT)) print(" (iii) python {} EF T EC EA NA ED ND nEF)".format(sys.argv[0])) print(" Calculate EF dependences of ne, nh, ND+, NA-") print(" nEF: # of EF mesh") print(" ex: python {} EF {} {} {} {} {} {} {}".format(sys.argv[0], T0, dos.EC, dos.EA, dos.NA, dos.ED, dos.ND, dos.nEF)) def updatevars(): global mode, dos, outcsvfile global T0 global dEFmin, dEFmax global Tmin, Tmax, nT, Tstep global dE_me_lsq argv = sys.argv if len(argv) == 1: usage() exit() mode = getarg( 1, mode) if mode == 'me': dE_me_lsq = getfloatarg(2, dE_me_lsq) elif mode == 'EF': T0 = getfloatarg(2, T0) dos.EC = getfloatarg(3, dos.EC) dos.EA = getfloatarg(4, dos.EA) dos.NA = getfloatarg(5, dos.NA) dos.ED = getfloatarg(6, dos.ED) dos.ND = getfloatarg(7, dos.ND) dos.nEF = getintarg(8, dos.nEF) dos.Eg = dos.EC - dos.EV elif mode == 'T': dos.EC = getfloatarg(2, dos.EC) dos.EA = getfloatarg(3, dos.EA) dos.NA = getfloatarg(4, dos.NA) dos.ED = getfloatarg(5, dos.ED) dos.ND = getfloatarg(6, dos.ND) dos.Eg = dos.EC - dos.EV Tmin = getfloatarg(7, Tmin) Tmax = getfloatarg(8, Tmax) nT = getintarg(9, nT) Tstep = (Tmax - Tmin) / (nT - 1) outcsvfile = f'EF-T-semi_FEA-{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, ydQ = dos.cal_T_lists(xT, print_level = 1) yEaNe = [] yEaNeTh = [] yEaNh = [] yEaNhTh = [] print("") print ("%8s %10s %12s %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", "dQ")) 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 %12g" % (yNe[i], yNh[i], yNDp[i], yNAm[i], yRH[i], yNs[i], ydQ[i])) print("") print("Save simulation data to [{}]".format(outcsvfile)) save_csv(outcsvfile, ["T(K)", "1000/T(K^-1)", "EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH", "Ns", "Ea(Ne,meV)", "Ea(Ne,T1/2)", "Ea(Nh,meV)", "Ea(Nh,T1/2)", "dQ"], [xT, xTinv, yEF, yNe, yNh, yNDp, yNAm, yRH, yNs, yEaNe, yEaNeTh, yEaNh, yEaNhTh, ydQ]) #============================= # Plot graphs #============================= fig = plt.figure(figsize = figsize) 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 print("") print("EF dependence at {} K".format(T0)) print(" EF range: {} - {}, {:12.6g} eV step, nE={}".format(dos.EFmin, dos.EFmax, dos.EFstep, dos.nEF)) xEF = [dos.EFmin + i * dos.EFstep for i in range(dos.nEF)] yNh, yNe, yNAm, yNDp, yRH, yNs, yNsabs = dos.cal_EF_lists(T0, xEF, print_level = 1) yT0Ne = [] yT0Nh = [] eps = 1.0e-300 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(dos.nEF): EF = xEF[i] if i == 0: i1 = i i2 = i+1 elif i == dos.nEF - 1: i1 = i-1 i2 = i else: i1 = i-1 i2 = i+1 dEF = xEF[i2] - xEF[i1] dlnNedEF = (log(yNe[i2]+eps) - log(yNe[i1]+eps)) / dEF dlnNhdEF = (log(yNh[i2]+eps) - log(yNh[i1]+eps)) / 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])) print("") Emid = (dos.EV + dos.EC) / 2.0 imid = int((Emid - dos.EFmin) / dos.EFstep + 1.0e-5) Eimid = dos.EFmin + imid * dos.EFstep print("Effective density of states at the mid gap {} eV (i={}, E={:8.4g}):".format(Emid, imid, Eimid)) 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) 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(" me* = {:8.4g} me".format(meeff)) print(" mh* = {:8.4g} me".format(mheff)) yNeapprox = [] yNhapprox = [] for i in range(dos.nEF): EF = xEF[i] kexpc = (dos.EC - EF) * e / kB / T0 kexpv = (EF - dos.EV) * e / kB / T0 yNeapprox.append(dos.NC * exp(-kexpc)) yNhapprox.append(dos.NV * exp(-kexpv)) print("") print("Save simulation 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 = figsize) 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([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(xEF, yNeapprox, label = 'Ne(Boltz)', color = 'gray', linewidth = 0.5) #, marker = 'x') ax3.plot(xEF, 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 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 = figsize) 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.Eg = dos.EC - dos.EV dos.NC = dos.meff2NC_FEA(dos.meeff, T0) dos.NV = dos.meff2NC_FEA(dos.mheff, T0) dos.DC0 = dos.meff2DC0_FEA(dos.meeff, T0) dos.DV0 = dos.meff2DC0_FEA(dos.mheff, T0) print("") print("Eg: {} eV".format(dos.Eg)) print("Valence band:") print(" EV : {} eV".format(dos.EV)) print(" mh*: {} me".format(dos.mheff)) print(" NV : {} me at {} K".format(dos.NV, T0)) print(" DV0: {} me at {} K".format(dos.DV0, T0)) print("Conduction band:") print(" EC : {} eV".format(dos.EC)) print(" me*: {} me".format(dos.meeff)) print(" NC : {} me at {} K".format(dos.NC, T0)) print(" DC0: {} me at {} K".format(dos.DC0, T0)) print("") print("Acceptor states:") print(" {} cm^-3 at {} eV".format(dos.NA, dos.EA)) print("Donoar states:") print(" {} cm^-3 at {} eV".format(dos.ND, dos.ED)) dos.EFmin = dos.EV + dEFmin dos.EFmax = dos.EC + dEFmax dos.EFstep = (dos.EFmax - dos.EFmin) / (dos.nEF - 1) dos.Emin = dos.EV + dEFmin dos.Emax = dos.EC + dEFmax dos.nE = int((dos.Emax - dos.Emin) / dos.Estep + 1.00001) dos.Elist = [dos.Emin + i * dos.Estep for i in range(dos.nE)] dos.DOSlist = [dos.DOS(dos.Elist[i]) for i in range(dos.nE)] print("") print("DOS range: {:8.4g} - {:8.4g}, {:8.4g} steps, nE={}".format(dos.Emin, dos.Emax, dos.Estep, dos.nE)) print("") print("Integration configuration") print(" Integration E range: {} - {} eV, or EF+-{}*kBT eV".format(dos.Einteg0, dos.Einteg1, dos.nrange)) 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()