import os import sys import signal import numpy as np from numpy import sin, cos, tan, pi, exp, log from matplotlib import pyplot as plt from scipy import optimize # newton関数はscipy.optimizeモジュールに入っている from scipy.optimize import minimize sys.path.append("c:/tkProg/tklib/python") sys.path.append("d:/tkProg/tklib/python") from tklib.tkfile import tkFile from tklib.tkutils import terminate, getarg, getintarg, getfloatarg, save_csv from tklib.tkutils import is_exist, is_file, is_dir #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, Ea_Arrhenius from tklib.tktransport.tkTransport import read_Hall_excel from tklib.tktransport.tkTransport import integrate_Simpson_list, fe, fh, NC2meff_FEA from tklib.tktransport.tkDOS_FEA import integrate_Simpson, integrate_Simpson_list, tkDOS from tklib.tktransport.tkmobility_tau import tkMobility, split_optstr """ Fit N-T Hall data with Free Electron Approximation """ usage_str = ''' "" "Usage: Variables in () are optional" " python {} mode (file args)".format(sys.argv[0]) " (i) python {} sim infile(.xlsx) EC EA NA ED ND".format(sys.argv[0]) " Plot N-T with simulated data" " ex: python {} sim {} {} {} {} {} {}".format(sys.argv[0], cparams.infile, dos.EC, dos.EA, dos.NA, dos.ED, dos.ND) " (ii) python {} fit|fit_h|fit_e infile(.xlsx) EC EA NA ED ND Eb a1 a2)".format(sys.argv[0]) " Fit to N-mu-T data" " mode = fit : Fit to Ns = 1/eRH by two-carrier (hole-electron) model. Scattering factor is fixed to r = 0.5" " fit_h : Fit to Nh by single-carrier model. Scattering factor is fixed to r = 0.5" " fit_e : Fit to Ne by single-carrier model. Scattering factor is fixed to r = 0.5" " fit_RH: Fit to RH by two-carrier (hole-electron) model with scattering factor" " ex: python {} fit {} {} {} {} {} {}".format(sys.argv[0], cparams.infile, dos.EC, dos.EA, dos.NA, dos.ED, dos.ND) '''[1:-1] def pint(str): return int(str) def pfloat(str): return float(str) def sigint_handler(signal, frame): print ('KeyboardInterrupt is caught') terminate(usage = usage) signal.signal(signal.SIGINT, sigint_handler) # EPS to avoid log(0.0) => log(x + cparams.eps) eps = 1.0e-700 def initialize(): #================================ # Global variables #================================ app = tkApplication(usage_str = usage_str) argv, narg = app.get_argv() app.dos = tkDOS() app.mu = tkMobility() cparams = app.get_params() cparams.debug = 0 cparams.use_simple = 0 cparams.print_level = 0 # mode: 'sim', 'fit', 'fit_e', 'fit_h', "fit_RH" cparams.mode = 'sim' # cparams.mode = 'fit' # cparams.mode = 'fit_e' # cparams.mode = 'fit_h' # data file to fit (.xlsx) cparams.infile = None cparams.parameterfile = None cparams.parameterbkfile = None cparams.outcsvfile = None dos = app.dos # Assume me = mh to calculate RH for two-carrier model dos.cal_mobility_e = lambda T: 1.0 dos.cal_mobility_h = lambda T: 1.0 dos.cal_mobility_e_description = "Same mobility for hole and electron" dos.cal_mobility_h_description = "Same mobility for hole and electron" 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.0e17 cparams.EF0 = 0.0 # initial EF to find EF #移動度パラメータ mobility = app.mu # For debug purpose. 1 will use 3-terms polynomial for the inverse of mobility mobility.debug = cparams.debug mobility.use_simple = cparams.use_simple #mobility.use_simple = 1 # charge, effective mass, scattering factor mobility.charge = 1.0 # in e, 1.0 for hole, -1.0 for electron mobility.meff = dos.meeff # in me mobility.rfac = 0.5 # tau = (meff/2)^0.5 * l0(T) * E^(r-0.5) mobility.l0 = 1.0e-8 # m dos.varkeys = ["meeff", "mheff", "EV", "EC", "EA", "NA", "ED", "ND"] # For mode = 'me'. E range for fitting to me*(DOS) cparams.dE_me_lsq = 0.5 # measured from EC/EV, eV # EF range for mode = 'EF' cparams.T0 = 300 # K cparams.dEFmin = -1.0 # measured from EV, eV cparams.dEFmax = 1.0 # measured from EC, eV cparams.nEF = 50 cparams.Estep = 0.01 # Integration step, eV # Temperature range for mode = 'T' cparams.Tmin = 300 # K cparams.Tmax = 600 cparams.nT = 11 cparams.Tstep = None # Integration range w.r.t. kBT cparams.Einteg0 = -3.0 cparams.Einteg1 = 3.0 cparams.nrange = 6.0 # Bisection method parameters # Initial search range cparams.dE_bisec = 0.4 # EPS, max iteration number, Output cycle cparams.eps_bisec = 1.0e-7 cparams.nmaxiter_bisec = 200 cparams.iprintiterval_bisec = 1 # 最適化パラメータの初期値 dos.varname = ["meeff", "mheff", "EV", "EC", "EA", "NA", "ED", "ND"] dos.varunit = [ "me", "me", "eV", "eV", "eV", "cm-3", "eV", "cm-3"] dos.ai0 = dos.parameter_list() dos.optid = [ 1, 1, 0, 0, 1, 1, 0, 0] dos.xT = [] dos.ys = [] dos.ymu = [] dos.yNhini = [] dos.yNeini = [] dos.yNsini = [] dos.yNsabsini = [] #============================================= # scipy.optimize.minimizeで使うアルゴリズム #============================================= #nelder-mead Downhill simplex #powell Modified Powell #cg conjugate gradient (Polak-Ribiere method) #bfgs BFGS法 #newton-cg Newton-CG #trust-ncg 信頼領域 Newton-CG 法 #dogleg 信頼領域 dog-leg 法 #L-BFGS-B’ (see here) #TNC’ (see here) #COBYLA’ (see here) #SLSQP’ (see here) #trust-constr’(see here) #dogleg’ (see here) #trust-exact’ (see here) #trust-krylov’ (see here) cparams.method = "nelder-mead" #method = 'cg' #method = 'powell' #method = 'bfgs' cparams.maxiter = 300 cparams.tol = 1.0e-4 cparams.h_diff = 1.0e-3 cparams.outputinterval = 1 #============================= # Graph configuration #============================= cparams.fig = None cparams.figsize = [12, 8] cparams.fontsize = 14 cparams.legend_fontsize = 12 cparams.graphupdateinterval = 10 return app, cparams, app.dos, app.mu def print_parameters(varname, varunit, ai, optid): for i in range(len(varname)): print(" {:10s}: {:14.8g} {:10}: optid={}".format(varname[i], ai[i], varunit[i], optid[i])) #============================= # Treat argments #============================= def usage(app): # app.usage(infile = cparams.infile) for s in app.usage_str.split('\n'): cmd = 'print({})'.format(s.rstrip()) eval(cmd) def updatevars(app, cparams, dos, mu): argv = sys.argv if len(argv) == 1: usage(app) exit() cparams.mode = getarg( 1, cparams.mode) cparams.infile = getarg( 2, cparams.infile) 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.Eg = dos.EC - dos.EV ai0 = dos.parameter_list() if cparams.infile is None: terminate("Error in updatevars: infile must be specified", usage = usage) header, ext = os.path.splitext(cparams.infile) filebody = os.path.basename(header) cparams.infile = filebody + '.xlsx' if not is_file(cparams.infile): cparams.infile = filebody + '.csv' cparams.parameterfile = filebody + '.in' cparams.parameterbkfile = filebody + '.in.bak' cparams.outcsvfile = filebody + f'-EFT-{cparams.mode}.csv' #============================= # other functions #============================= def recover_parameters(ais, optid, aidef0): aidef = aidef0.copy() ai = [] c = 0 for i in range(len(optid)): if optid[i] == 1: ai.append(ais[c]) c += 1 else: ai.append(aidef[i]) return ai def choose_parameters(ais0, optid): ais = ais0.copy() #ai0 = [mu0, Eb, a1, a2] ai = [] c = 0 for i in range(len(optid)): if optid[i] == 1: ai.append(ais[i]) return ai def calS2(app, cparams, dos, mobility, ai): aiorg = dos.parameter_list() ai0 = dos.parameter_list() ainew = recover_parameters(ai, dos.optid, ai0) dos.set_parameters(ainew) xTinv, yEFfin, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, yNsabsfin, ydQ = dos.cal_T_lists(dos.xT, print_level = cparams.print_level) # print("") if cparams.print_level >= 1: print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}" .format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', '|Ns,fin|')) for i in range(len(dos.xT)): print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}" .format(dos.xT[i], xTinv[i], yEFfin[i], dos.yN[i], yNefin[i], yNhfin[i], yNAmfin[i], yNDpfin[i], yRHfin[i], yNsabsfin[i])) if cparams.mode == 'fit_h': yNc = yNhfin elif cparams.mode == 'fit_e': yNc = yNefin else: yNc = yNsabsfin S2 = 0.0 for iT in range(len(dos.xT)): T = dos.xT[iT] S2 += (log(dos.yN[iT]+eps) - log(yNc[iT]+eps))**2 S2 /= (len(dos.xT) - 1) dos.set_parameters(aiorg) return S2 def calS2RH(app, cparams, dos, mobility, ai): aiorg = dos.parameter_list() ai0 = dos.parameter_list() ainew = recover_parameters(ai, dos.optid, ai0) dos.set_parameters(ainew) # xTinv, yEFfin, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, yNsabsfin, ydQ = dos.cal_T_lists(dos.xT, print_level = cparams.print_level) inf = dos.cal_T_lists2(dos.xT, mobility, print_level = cparams.print_level) xTinv = inf['Tinv'] yEFfin = inf['EF'] yNefin = inf['n_e'] yNhfin = inf['n_h'] # yNeHallfin = inf['nHall_e'] # yNhHallfin = inf['nHall_h'] yNsfin = inf['n_s'] yNsabsfin = inf['n_s(abs)'] yNAmfin = inf['NAm'] yNDpfin = inf['NDp'] yRHfin = inf['RH'] # yRHefin = inf['RH_e'] # yRHhfin = inf['RH_h'] # yFHefin = inf['FHall_e'] # yFHhfin = inf['FHall_h'] # print("") if cparams.print_level >= 1: print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}" .format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', '|Ns,fin|')) for i in range(len(dos.xT)): print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}" .format(dos.xT[i], xTinv[i], yEFfin[i], dos.yN[i], yNefin[i], yNhfin[i], yNAmfin[i], yNDpfin[i], yRHfin[i], yNsabsfin[i])) if cparams.mode == 'fit_RH': yNc = yRHfin if cparams.mode == 'fit_h': yNc = yNhfin elif cparams.mode == 'fit_e': yNc = yNefin else: yNc = yNsabsfin S2 = 0.0 for iT in range(len(dos.xT)): T = dos.xT[iT] if cparams.mode == 'fit_RH': if dos.yRH[iT] < 0.0: lnRH = -log(-dos.yRH[iT]) else: lnRH = log(dos.yRH[iT]+eps) if yRHfin[iT] < 0.0: lnRHfin = -log(-yRHfin[iT]) else: lnRHfin = log(yRHfin[iT]+eps) # S2 += (log(dos.yRH[iT]+eps) - log(yRHfin[iT]+eps))**2 S2 += (lnRH - lnRHfin)**2 else: S2 += (log(dos.yN[iT]+eps) - log(yNc[iT]+eps))**2 S2 /= (len(dos.xT) - 1) dos.set_parameters(aiorg) return S2 # First derivatives to be used e.g. for cg method # Approximate by forward difference method with the delta h = h_diff def diff1(app, cparams, dos, mu, ai): n = len(ai) f0 = calS2(ai) df = np.empty(n) for i in range(n): aii = ai.copy() aii[i] = ai[i] + h_diff df[i] = (calS2(app, cparams, mu, aii) - f0) / h_diff return df def plot(app, cparams, dos, mu, xTinv, yN, yNhini, yNeini, yNsabsini, yNhfin = None, yNefin = None, yNsabsfin = None, ymu = None, ymu_fin = None, yNhHallfin = None, yNeHallfin = None): kplot = [0.01, 2.0] plt = app.plt fig = app.fig fontsize = cparams.fontsize legend_fontsize = cparams.legend_fontsize plt.clf() ax1 = fig.add_subplot(2, 2, 1) ax2 = fig.add_subplot(2, 2, 2) axmu = fig.add_subplot(2, 2, 3) plt.title(cparams.infile, fontsize = fontsize) ax1.plot(xTinv, yN, label = 'N(obs)', linestyle = 'none', marker = 'o') ax1.plot(xTinv, yNhini, label = '$N_h$(ini)', linestyle = 'dashed', color = 'red', linewidth = 0.5) ax1.plot(xTinv, yNeini, label = '$N_e$(ini)', linestyle = 'dashed', color = 'blue', linewidth = 0.5) ax1.plot(xTinv, yNsabsini, label = '|$N_s$(ini)|', linestyle = '-', color = 'green', linewidth = 0.5) ax1.set_yscale('log') ax1.set_ylim([min(yN) * kplot[0], max(yN) * kplot[1]]) ax1.set_xlabel('1000/T (K$^{-1}$)', fontsize = fontsize) ax1.set_ylabel('N (cm$^{-3}$)', fontsize = fontsize) ax1.legend(fontsize = legend_fontsize) ax1.tick_params(labelsize = fontsize) if yNhfin: print_parameters(dos.varname, dos.varunit, dos.parameter_list(), dos.optid) ax2.plot(xTinv, yN, label = '$N$(obs,Hall)', linestyle = 'none', marker = 'o') ax2.plot(xTinv, yNhfin, label = '$N_h$(fin)', linestyle = 'dashed', color = 'red', linewidth = 0.5) ax2.plot(xTinv, yNefin, label = '$N_e$(fin)', linestyle = 'dashed', color = 'blue', linewidth = 0.5) ax2.plot(xTinv, yNsabsfin, label = '|$N_s$(fin)|', linestyle = '-', color = 'green', linewidth = 0.5) if yNhHallfin: ax2.plot(xTinv, yNhHallfin, label = '$N_{h,Hall}$(fin)', linestyle = 'dotted', color = 'red', linewidth = 0.3) ax2.plot(xTinv, yNeHallfin, label = '$N_{e,Hall}$(fin)', linestyle = 'dotted', color = 'blue', linewidth = 0.3) if yNhfin: ax2.set_yscale('log') ax2.set_ylim([min(yN) * kplot[0], max(yN) * kplot[1]]) ax2.set_xlabel('1000/T (K$^{-1}$)', fontsize = fontsize) ax2.set_ylabel('N (cm$^{-3}$)', fontsize = fontsize) ax2.legend(fontsize = legend_fontsize) ax2.tick_params(labelsize = fontsize) if ymu is not None: axmu.plot(dos.xT, dos.ymu, label = '$\mu$(obs,Hall)', color = 'black', linewidth = 0.5, linestyle = '-', marker = 'o') if ymu_fin is not None: axmu.plot(dos.xT, ymu_fin, label = '$\mu$(fin)', color = 'blue', linewidth = 0.5, linestyle = 'dashed', marker = '^') axmu.set_xlabel("T (K)", fontsize = fontsize) axmu.set_ylabel("$\mu$ (cm$^2/(Vs)$)", fontsize = fontsize) axmu.legend(fontsize = legend_fontsize) axmu.tick_params(labelsize = fontsize) plt.tight_layout() plt.pause(0.01) # Callback function for scipy.optimize.minimize() # Print variables every iteration, and update graph for every graphupdateinterval iterationsおｔ def callback(app, cparams, dos, mobility, xk): S2 = calS2(app, cparams, dos, mu, xk) ai0 = dos.parameter_list() ainew = recover_parameters(xk, dos.optid, ai0) dos.set_parameters(ainew) if app.iter % cparams.outputinterval == 0: print("callback {:04d}: ".format(app.iter), end = '') for i in range(len(ainew)): print("{:10.4g} ".format(ainew[i]), end = '') print(" S2={:10.6g}".format(S2)) if app.iter % cparams.graphupdateinterval == 0: ymufin = [] xTinv, yEFfin, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, yNsabsfin, ydQ \ = dos.cal_T_lists(dos.xT, print_level = cparams.print_level) ymu_fit = [] for i in range(len(dos.xT)): ymu_fit.append(dos.ys[i] / e / (yNefin[i] + yNhfin[i])) plot(app, cparams, dos, mobility, xTinv, dos.yN, dos.yNhini, dos.yNeini, dos.yNsabsini, yNhfin, yNefin, yNsabsfin, ymu = dos.ymu, ymu_fin = ymu_fit) app.iter += 1 try: pass except KeyboardInterrupt: print ('KeyboardInterrupt exception is caught') terminate(usage = usage) def callbackRH(app, cparams, dos, mobility, xk): S2 = calS2RH(app, cparams, dos, mu, xk) ai0 = dos.parameter_list() ainew = recover_parameters(xk, dos.optid, ai0) dos.set_parameters(ainew) if app.iter % cparams.outputinterval == 0: print("callback {:04d}: ".format(app.iter), end = '') for i in range(len(ainew)): print("{:10.4g} ".format(ainew[i]), end = '') print(" S2={:10.6g}".format(S2)) if app.iter % cparams.graphupdateinterval == 0: ymufin = [] # xTinv, yEFfin, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, yNsabsfin, ydQ = dos.cal_T_lists(dos.xT, print_level = cparams.print_level) inf = dos.cal_T_lists2(dos.xT, mobility, print_level = cparams.print_level) xTinv = inf['Tinv'] # yEFfin = inf['EF'] yNefin = inf['n_e'] yNhfin = inf['n_h'] yNeHallfin = inf['nHall_e'] yNhHallfin = inf['nHall_h'] yNsfin = inf['n_s'] yNsabsfin = inf['n_s(abs)'] # yNAmfin = inf['NAm'] # yNDpfin = inf['NDp'] # yRHfin = inf['RH'] # yRHefin = inf['RH_e'] # yRHhfin = inf['RH_h'] # yFHefin = inf['FHall_e'] # yFHhfin = inf['FHall_h'] ymu_fit = [] for i in range(len(dos.xT)): ymu_fit.append(dos.ys[i] / e / (yNefin[i] + yNhfin[i])) plot(app, cparams, dos, mobility, dos.xTinv, dos.yN, dos.yNhini, dos.yNeini, dos.yNsabsini, yNhfin, yNefin, yNsabsfin, ymu = dos.ymu, ymu_fin = ymu_fit, yNhHallfin = yNhHallfin, yNeHallfin = yNeHallfin) app.iter += 1 try: pass except KeyboardInterrupt: print ('KeyboardInterrupt exception is caught') terminate(usage = usage) #======================================================== #======================================================== #======================================================== # Fitting to Hall - T data by ne,Hall, nh,Hall and ns # Scattering factor is fixed to r = 0.5 #======================================================== #======================================================== #======================================================== def fit(app, cparams, dos, mobility): print("Fitting to mu-T Hall data by nonlinear least-squares method") print(" r: {}".format(mobility.rfac)) # read data print("") print("Read Hall data from {}".format(cparams.infile)) header, xT, yN, ymu, ys = dos.read_excel_Hall(cparams.infile, usage = usage) # header, xT, yN, ymu, ys, yRH = dos.read_excel_Hall2(cparams.infile, usage = usage) print(" Headers:", header) if xT is None: terminate(f"Temperature data is not found in [{cparams.infile}]", usage = usage) if yN is None: terminate(f"Carrier density data is not found in [{cparams.infile}]", usage = usage) if ymu is None: terminate(f"Mobility data is not found in [{cparams.infile}]", usage = usage) if ys is None: terminate(f"Conductivity data is not found in [{cparams.infile}]", usage = usage) # calculate initial values xTinv = [] ymuini = [] print("") print("Calculate mu-T from the initial parameters") # print("{:8}\t{:12}".format("T(K)", "mu(cm2/Vs)")) dos.xTinv, dos.yEFini, dos.yNeini, dos.yNhini, dos.yNAmini, dos.yNDpini, dos.yRHini, dos.yNsini, dos.yNsabsini, dos.ydQ \ = dos.cal_T_lists(xT, print_level = cparams.print_level) print("") print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}" .format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', 'Ns,fin')) for i in range(len(xT)): print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}" .format(dos.xT[i], dos.xTinv[i], dos.yEFini[i], dos.yN[i], dos.yNeini[i], dos.yNhini[i], dos.yNAmini[i], dos.yNDpini[i], dos.yRHini[i], dos.yNsini[i])) # exit() #============================= # グラフの表示 #============================= print("") print("plot") app.plt = plt app.fig = plt.figure(figsize = cparams.figsize) plot(app, cparams, dos, mobility, dos.xTinv, dos.yN, dos.yNhini, dos.yNeini, dos.yNsabsini) #============================= # Optimization #============================= print("") print("") print("Nonlinear least-squares fitting by ", cparams.method) print(" tol=", cparams.tol) app.iter = 0 ai0 = dos.parameter_list() ai = choose_parameters(ai0, dos.optid) ret = minimize(lambda xi: calS2(app, cparams, dos, mobility, xi), ai, method = cparams.method, jac = lambda xi: diff1(app, cparams, dos, mobility, xi), tol = cparams.tol, callback = lambda xk: callback(app, cparams, dos, mobility, xk), options = {'maxiter': cparams.maxiter, "disp": True}) # print("ret=", ret) if cparams.method == 'nelder-mead': simplex = ret['final_simplex'] ai = simplex[0][0] fmin = ret['fun'] elif cparams.method == 'cg': ai = ret['x'] fmin = ret['fun'] elif cparams.method == 'powell': ai = ret['x'] fmin = ret['fun'] elif cparams.method == 'bfgs': ai = ret['x'] fmin = CalS2(ai) aifin = recover_parameters(ai, dos.optid, ai0) dos.set_parameters(aifin) print("") print("Optimized at S2={:12.6g}:".format(fmin)) print_parameters(dos.varname, dos.varunit, dos.parameter_list(), dos.optid) dos.NC = dos.meff2NC_FEA(dos.meeff, cparams.T0) dos.NV = dos.meff2NC_FEA(dos.mheff, cparams.T0) dos.DC0 = dos.meff2DC0_FEA(dos.meeff, cparams.T0) dos.DV0 = dos.meff2DC0_FEA(dos.mheff, cparams.T0) print(" NV : {} me at {} K".format(dos.NV, cparams.T0)) print(" DV0: {} me at {} K".format(dos.DV0, cparams.T0)) print(" NC : {} me at {} K".format(dos.NC, cparams.T0)) print(" DC0: {} me at {} K".format(dos.DC0, cparams.T0)) print("") print("Save to [{}]".format(cparams.parameterfile)) cparams.save_parameters(cparams.parameterfile, "Preferences") dos.save_parameters(cparams.parameterfile, dos.optid, {"last_fitmode": cparams.mode, "rfac": mobility.rfac, "l0": mobility.l0, 'S2N': fmin}) print("") print("Calculate final data") xTinv, yEFfin, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, yNsabsfin, ydQ \ = dos.cal_T_lists(xT, print_level = cparams.print_level) print("") print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}" .format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', 'Ns,fin', 'mu(Hall)', 'mu,fin')) ymu_fit = [] for i in range(len(xT)): ymu_fit.append(ys[i] / e / (yNefin[i] + yNhfin[i])) print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}" .format(xT[i], xTinv[i], yEFfin[i], yN[i], yNefin[i], yNhfin[i], yNAmfin[i], yNDpfin[i], yRHfin[i], yNsfin[i], ymu[i], ymu_fit[i])) print("Save fitting data to [{}]".format(cparams.outcsvfile)) save_csv(cparams.outcsvfile, ['T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', 'Ns,fin', 'mu(Hall)(cm2/Vs)', 'mu,fin'], [xT, xTinv, yEFfin, yN, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, ymu, ymu_fit]) #============================= # グラフの表示 #============================= print("") print("Plot optimized") plot(app, cparams, dos, mobility, dos.xTinv, dos.yN, dos.yNhini, dos.yNeini, dos.yNsabsini, yNhfin, yNefin, yNsabsfin, ymu = ymu, ymu_fin = ymu_fit) print("") print("Press ENTER to exit>>", end = '') input() print("") app.terminate(usage = usage) #======================================================== #======================================================== #======================================================== # Fitting to Hall - T data by RH # Scattering factor is taken from mobility.rfac #======================================================== #======================================================== #======================================================== def fit_RH(app, cparams, dos, mobility): print("Fitting to mu-T Hall data by nonlinear least-squares method") print(" r: {}".format(mobility.rfac)) print(" Calculate two-carrier model RH:") print(f" mu_h: {dos.cal_mobility_h_description}") print(f" mu_e: {dos.cal_mobility_e_description}") # read data print("") print("Read Hall data from {}".format(cparams.infile)) header, xT, yN, ymu, ys, yRH = dos.read_excel_Hall2(cparams.infile, usage = usage) print(" Headers:", header) # print(" xT :", xT) # print(" yN :", yN) # print(" ymu:", ymu) # print(" ys :", ys) # print(" yRH:", yRH) if xT is None: app.terminate(f"Temperature data is not found in [{cparams.infile}]", usage = usage) if yN is None: app.terminate(f"Carrier density data is not found in [{cparams.infile}]", usage = usage) if ymu is None: app.terminate(f"Mobility data is not found in [{cparams.infile}]", usage = usage) if ys is None: app.terminate(f"Conductivity data is not found in [{cparams.infile}]", usage = usage) if yRH is None: app.terminate(f"RH data is not found in [{cparams.infile}]", usage = usage) # calculate initial values xTinv = [] ymuini = [] print("") print("Calculate mu-T from the initial parameters") # print("{:8}\t{:12}".format("T(K)", "mu(cm2/Vs)")) # xTinv, yEFini, yNeini, yNhini, yNAmini, yNDpini, yRHini, yNsini, yNsabsini, ydQ = dos.cal_T_lists(xT, print_level = cparams.print_level) inf = dos.cal_T_lists2(xT, mobility, print_level = cparams.print_level) dos.xTinv = inf['Tinv'] dos.yEFini = inf['EF'] dos.yNeini = inf['n_e'] dos.yNhini = inf['n_h'] dos.yNeHallini = inf['nHall_e'] dos.yNhHallini = inf['nHall_h'] dos.yNsini = inf['n_s'] dos.yNsabsini = inf['n_s(abs)'] dos.yNAmini = inf['NAm'] dos.yNDpini = inf['NDp'] dos.yRHini = inf['RH'] dos.yRHeini = inf['RH_e'] dos.yRHhini = inf['RH_h'] dos.yFHeini = inf['FHall_e'] dos.yFHhini = inf['FHall_h'] print("") print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}" .format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,ini', 'Nh,ini', 'NA-,ini', 'ND+,ini', 'RH,ini', 'Ns,ini')) for i in range(len(xT)): print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}" .format(dos.xT[i], dos.xTinv[i], dos.yEFini[i], dos.yN[i], dos.yNeini[i], dos.yNhini[i], dos.yNAmini[i], dos.yNDpini[i], dos.yRHini[i], dos.yNsini[i])) # exit() #============================= # グラフの表示 #============================= print("") print("plot") app.plt = plt app.fig = plt.figure(figsize = cparams.figsize) plot(app, cparams, dos, mobility, dos.xTinv, dos.yN, dos.yNhini, dos.yNeini, dos.yNsabsini) #============================= # Optimization #============================= print("") print("") print("Nonlinear least-squares fitting by ", cparams.method) print(" tol=", cparams.tol) app.iter = 0 ai0 = dos.parameter_list() ai = choose_parameters(ai0, dos.optid) ret = minimize(lambda xk: calS2RH(app, cparams, dos, mobility, xk), ai, method = cparams.method, jac = lambda xk: diff1(app, cparams, dos, mobility, xk), tol = cparams.tol, callback = lambda xk: callbackRH(app, cparams, dos, mobility, xk), options = {'maxiter': cparams.maxiter, "disp":True}) # print("ret=", ret) if cparams.method == 'nelder-mead': simplex = ret['final_simplex'] ai = simplex[0][0] fmin = ret['fun'] elif cparams.method == 'cg': ai = ret['x'] fmin = ret['fun'] elif cparams.method == 'powell': ai = ret['x'] fmin = ret['fun'] elif cparams.method == 'bfgs': ai = ret['x'] fmin = CalS2(ai) aifin = recover_parameters(ai, dos.optid, ai0) dos.set_parameters(aifin) print("") print("Optimized at S2={:12.6g}:".format(fmin)) print_parameters(dos.varname, dos.varunit, dos.parameter_list(), dos.optid) dos.NC = dos.meff2NC_FEA(dos.meeff, cparams.T0) dos.NV = dos.meff2NC_FEA(dos.mheff, cparams.T0) dos.DC0 = dos.meff2DC0_FEA(dos.meeff, cparams.T0) dos.DV0 = dos.meff2DC0_FEA(dos.mheff, cparams.T0) print(" NV : {} me at {} K".format(dos.NV, cparams.T0)) print(" DV0: {} me at {} K".format(dos.DV0, cparams.T0)) print(" NC : {} me at {} K".format(dos.NC, cparams.T0)) print(" DC0: {} me at {} K".format(dos.DC0, cparams.T0)) print("") print("Save to [{}]".format(cparams.parameterfile)) cparams.save_parameters(cparams.parameterfile, "Preferences") dos.save_parameters(cparams.parameterfile, dos.optid, {"last_fitmode": cparams.mode, "rfac": mobility.rfac, "l0": mobility.l0, 'S2N': fmin}) print("") print("Calculate final data") # xTinv, yEFfin, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, yNsabsfin, ydQ = dos.cal_T_lists(xT, print_level = cparams.print_level) inf = dos.cal_T_lists2(xT, mobility, print_level = cparams.print_level) xTinv = inf['Tinv'] yEFfin = inf['EF'] yNefin = inf['n_e'] yNhfin = inf['n_h'] yNeHallfin = inf['nHall_e'] yNhHallfin = inf['nHall_h'] yNsfin = inf['n_s'] yNsabsfin = inf['n_s(abs)'] yNAmfin = inf['NAm'] yNDpfin = inf['NDp'] yRHfin = inf['RH'] yRHefin = inf['RH_e'] yRHhfin = inf['RH_h'] yFHefin = inf['FHall_e'] yFHhfin = inf['FHall_h'] print("") print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}" .format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', 'Ns,fin', 'mu(Hall)', 'mu,fin')) ymu_fit = [] for i in range(len(xT)): ymu_fit.append(ys[i] / e / (yNefin[i] + yNhfin[i])) print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}" .format(xT[i], xTinv[i], yEFfin[i], yN[i], yNefin[i], yNhfin[i], yNAmfin[i], yNDpfin[i], yRHfin[i], yNsfin[i], ymu[i], ymu_fit[i])) print("Save fitting data to [{}]".format(cparams.outcsvfile)) save_csv(cparams.outcsvfile, ['T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', 'Ns,fin', 'mu(Hall)(cm2/Vs)', 'mu,fin', 'FH,e', 'FH,h', 'Ne,Hall', 'Nh,Hall'], [xT, xTinv, yEFfin, yN, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, ymu, ymu_fit, yFHefin, yFHhfin, yNeHallfin, yNhHallfin]) #============================= # グラフの表示 #============================= print("") print("Plot optimized") plot(app, cparams, dos, mobility, dos.xTinv, dos.yN, dos.yNhini, dos.yNeini, dos.yNsabsini, yNhfin, yNefin, yNsabsfin, ymu = dos.ymu, ymu_fin = ymu_fit, yNhHallfin = yNhHallfin, yNeHallfin = yNeHallfin) print("") print("Press ENTER to exit>>", end = '') input() print("") app.terminate(usage = usage) #======================================================== #======================================================== #======================================================== # Simulate Hall-related data #======================================================== #======================================================== #======================================================== def sim(app, cparams, dos, mobility): print("") print("Simulate T dependences") print(" r: {}".format(mobility.rfac)) # read data print("") print("Read Hall data from {}".format(cparams.infile)) # header, xT, yN, ymu, ys = dos.read_excel_Hall(cparams.infile, usage = usage) header, xT, yN, ymu, ys, yRH = dos.read_excel_Hall2(cparams.infile, usage = usage) print(" Headers:", header) print(" xT :", xT) print(" yN :", yN) print(" ymu:", ymu) print(" ys :", ys) print(" yRH:", yRH) # calculate initial values print("") # xTinv, yEFini, yNeini, yNhini, yNAmini, yNDpini, yRHini, yNsini, yNsabsini, ydQ = dos.cal_T_lists(xT, print_level = cparams.print_level) inf = dos.cal_T_lists2(xT, mobility, print_level = cparams.print_level) xTinv = inf['Tinv'] yEFini = inf['EF'] yNeini = inf['n_e'] yNhini = inf['n_h'] yNeHallini = inf['nHall_e'] yNhHallini = inf['nHall_h'] yNsini = inf['n_s'] yNsabsini = inf['n_s(abs)'] yNAmini = inf['NAm'] yNDpini = inf['NDp'] yRHini = inf['RH'] yRHeini = inf['RH_e'] yRHhini = inf['RH_h'] yFHeini = inf['FHall_e'] yFHhini = inf['FHall_h'] yNT23 = [] yNT43 = [] print("Calculate n-T from the initial parameters") print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}" .format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,ini', 'Nh,ini', 'NA-,ini', 'ND+,ini', 'RH,ini', 'Ns,ini', 'mu(Hall)', 'mu,ini')) ymu_ini = [] for i in range(len(xT)): print("i=", i, ys[i], yNeini[i], yNhini[i]) ymu_ini.append(ys[i] / e / (yNeini[i] + yNhini[i])) yNT23.append(yN[i] * pow(xT[i], -2.0 / 3.0)) yNT43.append(yN[i] * pow(xT[i], -4.0 / 3.0)) print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}" .format(xT[i], xTinv[i], yEFini[i], yN[i], yNeini[i], yNhini[i], yNAmini[i], yNDpini[i], yRHini[i], yNsini[i], ymu[i], ymu_ini[i])) print("") # xTinv, yEFini, yNeini, yNhini, yNAmini, yNDpini, yRHini, yNsini, yNsabsini, ydQ = dos.cal_T_lists(xT, print_level = dos.print_level) Tinv, yEa = Ea_Arrhenius(xT, yN, eps = 1.0e-300, kTinv = 1000.0) Tinv, yEaT23 = Ea_Arrhenius(xT, yNT23, eps = 1.0e-300, kTinv = 1000.0) Tinv, yEaT43 = Ea_Arrhenius(xT, yNT43, eps = 1.0e-300, kTinv = 1000.0) nT = len(xT) print("Activation energy") print("{:8}\t{:14}\t{:14}\t{:14}".format("T(K)", "Ea(eV)(log(N))", "Ea(eV)(log(N*T^-2/3))", "Ea(eV)(log(N*T^-4/3))")) for iT in range(nT): print("{:8.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}".format(xT[iT], yEa[iT], yEaT23[iT], yEaT43[iT])) print("") print("Save to [{}]".format(cparams.parameterfile)) cparams.save_parameters(cparams.parameterfile, "Preferences") dos.save_parameters(cparams.parameterfile, dos.optid, {"rfac": mobility.rfac, "l0": mobility.l0}) print("Save simulation data to [{}]".format(cparams.outcsvfile)) save_csv(cparams.outcsvfile, ['T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,ini', 'Nh,ini', 'NA-,ini', 'ND+,ini', 'RH,ini', 'Ns,ini', 'mu(Hall)(cm2/Vs)', 'mu,ini', 'FH,e', 'FH,h', 'Ne,Hall', 'Nh,Hall'], [xT, xTinv, yEFini, yN, yNeini, yNhini, yNAmini, yNDpini, yRHini, yNsini, ymu, ymu_ini, yFHeini, yFHhini, yNeHallini, yNhHallini]) #============================= # グラフの表示 #============================= app.fig = plt.figure(figsize = cparams.figsize) fig = app.fig fontsize = cparams.fontsize legend_fontsize = cparams.legend_fontsize ax1 = fig.add_subplot(2, 3, 1) ax3 = fig.add_subplot(2, 3, 2) axmu = fig.add_subplot(2, 3, 3) ax2 = fig.add_subplot(2, 3, 4) ax4 = fig.add_subplot(2, 3, 5) axFH = fig.add_subplot(2, 3, 6) axEF = axFH.twinx() ax1.set_title(cparams.infile, fontsize = fontsize) ins1 = axFH.plot(xTinv, yFHeini, label = '$F_{Hall,e}$', color = 'blue', linestyle = '-', linewidth = 0.5) ins2 = axFH.plot(xTinv, yFHhini, label = '$F_{Hall,h}$', color = 'red', linestyle = '-', linewidth = 0.5) ins3 = axEF.plot(xTinv, yEFini, label = '$E_F$', color = 'black', linestyle = '-', linewidth = 0.5) ins4 = axEF.plot(axFH.get_xlim(), [dos.EV, dos.EV], label = '$E_V$', color = 'green', linestyle = 'dashed', linewidth = 0.3) ins5 = axEF.plot(axFH.get_xlim(), [dos.EC, dos.EC], label = '$E_C$', color = 'green', linestyle = 'dashed', linewidth = 0.3) axFH.set_xlabel("1000/T (K$^{-1}$)", fontsize = fontsize) axFH.set_ylabel("$F_{Hall}$ (cm$^{-3}$/C)", fontsize = fontsize) axEF.set_ylabel("$E_F$ (eV)", fontsize = fontsize) ins = ins1 + ins2 + ins3 + ins4 + ins5 axFH.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'upper center') #loc = 'best') # axFH.legend(fontsize = legend_fontsize) axFH.tick_params(labelsize = fontsize) ax1.plot(xTinv, yN, label = '$N(obs)$', color = 'black', linestyle = '', marker = 'o') ax1.plot(xTinv, yNeini, label = '$N_e(ini)$', color = 'blue', linewidth = 0.5, linestyle = '-') ax1.plot(xTinv, yNhini, label = '$N_h(ini)$', color = 'red', linewidth = 0.5, linestyle = '-') ax1.plot(xTinv, yNeHallini, label = '$N_{e,Hall}(ini)$', color = 'blue', linewidth = 0.5, linestyle = 'dashed') ax1.plot(xTinv, yNhHallini, label = '$N_{h,Hall}(ini)$', color = 'red', linewidth = 0.5, linestyle = 'dashed') ax1.plot(xTinv, yNsabsini, label = '$N_s(ini)$', color = 'purple', linewidth = 0.5, linestyle = 'dashed') ax1.set_xlabel("1000/T (K$^{-1}$)", fontsize = fontsize) ax1.set_ylabel("N (cm$^{-3}$)", fontsize = fontsize) ax1.set_yscale('log') ax1.legend(fontsize = legend_fontsize) ax1.tick_params(labelsize = fontsize) ax3.plot(xT, yN, label = '$N(obs)$', color = 'black', linestyle = '', marker = 'o') ax3.plot(xT, yNeini, label = '$N_e(ini)$', color = 'blue', linewidth = 0.5, linestyle = '-') ax3.plot(xT, yNhini, label = '$N_h(ini)$', color = 'red', linewidth = 0.5, linestyle = '-') ax3.plot(xT, yNsabsini, label = '$N_s(ini)$', color = 'purple', linewidth = 0.5, linestyle = 'dashed') ax3.set_xlabel("T (K)", fontsize = fontsize) ax3.set_ylabel("N (cm$^{-3}$)", fontsize = fontsize) ax3.set_yscale('log') ax3.legend(fontsize = legend_fontsize) ax3.tick_params(labelsize = fontsize) axmu.plot(xT, ymu, label = '$\mu$(Hall)', color = 'black', linewidth = 0.5, linestyle = '-', marker = 'o') axmu.plot(xT, ymu_ini, label = '$\mu$(ini)', color = 'blue', linewidth = 0.5, linestyle = 'dashed', marker = '^') axmu.set_xlabel("T (K)", fontsize = fontsize) axmu.set_ylabel("$\mu$ (cm$^2/(Vs)$)", fontsize = fontsize) axmu.legend(fontsize = legend_fontsize) axmu.tick_params(labelsize = fontsize) ax2.plot(xTinv, yEa, label = '$E_a(obs)$', color = 'black', linestyle = '-', linewidth = 0.5, marker = 'o') ax2.plot(xTinv, yEaT23, label = '$E_a$(obs)(log($N$$T^{-2/3}$))', color = 'blue', linestyle = '-', linewidth = 0.5, marker = '^') ax2.plot(xTinv, yEaT43, label = '$E_a$(obs)(log($N$$T^{-4/3}$))', color = 'red', linestyle = '-', linewidth = 0.5, marker = '+') ax2.set_xlabel("1000/T (K$^{-1}$)", fontsize = fontsize) ax2.set_ylabel("$E_a$ (eV)", fontsize = fontsize) ax2.legend(fontsize = legend_fontsize) ax2.tick_params(labelsize = fontsize) ax4.plot(xT, yEa, label = '$E_a(obs)$', color = 'black', linestyle = '-', linewidth = 0.5, marker = 'o') ax4.plot(xT, yEaT23, label = '$E_a$(obs)(log($N$$T^{-2/3}$))', color = 'blue', linestyle = '-', linewidth = 0.5, marker = '^') ax4.plot(xT, yEaT43, label = '$E_a$(obs)(log($N$$T^{-4/3}$))', color = 'red', linestyle = '-', linewidth = 0.5, marker = '+') ax4.set_xlabel("T (K)", fontsize = fontsize) ax4.set_ylabel("$E_a$ (eV)", fontsize = fontsize) ax4.legend(fontsize = legend_fontsize) ax4.tick_params(labelsize = fontsize) plt.tight_layout() plt.pause(0.1) print("Press ENTER to exit>>", end = '') input() app.terminate("", usage = usage) def main(app, cparams, dos, mobility): updatevars(app, cparams, dos, mu) print("") print("===============================================================================") print(" Calculate T dependences of semiconductor statistics using VASP results") print("===============================================================================") print("") print("infile : {}".format(cparams.infile)) print("parameter file : {}".format(cparams.parameterfile)) print("parameter backup file: {}".format(cparams.parameterbkfile)) print("mode: ", cparams.mode) print("") print("") print("Read [{}]".format(cparams.parameterfile)) cparams.read_parameters(cparams.parameterfile) dos.read_parameters(cparams.parameterfile, dos.optid, mobility) dos.set_misc_parameters(cparams.get_dict()) # check args again to use given parameters updatevars(app, cparams, dos, mobility) print("") print("Initial parameters") print_parameters(dos.varname, dos.varunit, dos.parameter_list(), dos.optid) dos.Eg = dos.EC - dos.EV dos.NC = dos.meff2NC_FEA(dos.meeff, cparams.T0) dos.NV = dos.meff2NC_FEA(dos.mheff, cparams.T0) dos.DC0 = dos.meff2DC0_FEA(dos.meeff, cparams.T0) dos.DV0 = dos.meff2DC0_FEA(dos.mheff, cparams.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, cparams.T0)) print(" DV0: {} me at {} K".format(dos.DV0, cparams.T0)) print("Conduction band:") print(" EC : {} eV".format(dos.EC)) print(" me*: {} me".format(dos.meeff)) print(" NC : {} me at {} K".format(dos.NC, cparams.T0)) print(" DC0: {} me at {} K".format(dos.DC0, cparams.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)) print("") print("Scattering factor: r = {}".format(mobility.rfac)) dos.EFmin = dos.EV + cparams.dEFmin dos.EFmax = dos.EC + cparams.dEFmax print("d=", dos.EFmin, dos.EFmax,dos.nEF) dos.EFstep = (dos.EFmax - dos.EFmin) / (dos.nEF - 1) dos.Emin = dos.EV + cparams.dEFmin dos.Emax = dos.EC + cparams.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 cparams.mode == 'sim': sim(app, cparams, dos, mobility) elif cparams.mode == 'fit_RH': fit_RH(app, cparams, dos, mobility) elif 'fit' in cparams.mode: fit(app, cparams, dos, mobility) else: app.terminate("Error in main: Invalide cparams.mode [{}]".format(cparams.mode), usage = usage) if __name__ == '__main__': app, cparams, dos, mu = initialize() main(app, cparams, dos, mu)