import sys import numpy as np from numpy import sin, cos, tan, pi, exp, log, sqrt from scipy.optimize import minimize import pandas as pd from matplotlib import pyplot as plt from tklib.tksci.tksci import log10, e, kB from tklib.tkutils import getarg, getintarg, getfloatarg, pconv_by_type, print_line, print_data from tklib.tkapplication import tkApplication from tklib.tkparams import tkParams from tklib.tkvariousdata import tkVariousData from tklib.tkfilter import tkFilter from tklib.tktransport.tkTransport import meff2NC_FEA, meff2DC0_FEA from tklib.tktransport.tkWeightedMobility import weighted_mobility, weighted_mobility_new, weighted_mobility_exact from tklib.tktransport.tkDOS_FEA import tkDOS from tklib.tktransport.tkmobility_tau import tkMobility, split_optstr from tklib.tkgraphic.tkplotevent import tkPlotEvent from tklib.tksci.tkFit_mxy import tkFit_mxy """ 半導体の２キャリアモデル解析（熱電＋伝導度） """ def initialize(): #================================ # Global variables #================================ app = tkApplication() argv, narg = app.get_argv() app.cparams = tkParams() app.config = tkParams() cparams = app.cparams config = app.config config.debug = 0 config.print_level = 0 config.module_dir = './filter' config.module_name = 'S-sigma-T2xlsx' app.dos = tkDOS() app.dos.EV = 0.0 app.mu_e = tkMobility() app.mu_e.debug = 0 app.mu_e.use_simple = 0 app.mu_e.charge = -1.0 # in e, 1.0 for hole, -1.0 for electron app.mu_h = tkMobility() app.mu_h.debug = 0 app.mu_h.use_simple = 0 app.mu_h.charge = 1.0 # in e, 1.0 for hole, -1.0 for electron app.varname = [ "EC", "ED", "ND", "EA", "NA", "l0e", "l0h", "re", "rh"] # app.varname = [ "EC", "ED", "log(ND)", "EA", "log(NA)", "log(l0)", "l0h", "re", "rh"] app.unit = [ "eV", "eV", "cm^-3", "eV", "cm^-3", "m", "m", "re", "rh"] app.pk_convert = [ "", "", "", "", "", "", "", "", ""] # app.pk_convert = [ "", "", "log", "", "log", "log", "l0h", "re", "rh"] app.optid = [ 0, 0, 0, 1, 1, 1, 1, 0, 0] app.kmin = [ 0.0, 0.1, 1.0, 0.0, 1.0, 1.0e-15, 1.0e-15, 0.0, 0.0] # app.kmin = [ 0.0, 0.1, 1.0, 0.0, 1.0, log(1.0e-12), "l0h", 2.0, 2.0] app.kmax = [ 2.2, 2.2, 1.0e23, 1.5, 1.0e23, 1.0e-4, 1.0e-4, 2.0, 2.0] # app.kmax = [ 2.2, 2.2, log(1.0e23), 1.5, log(1.0e23), log(1.0e-4), "l0h", "re", "rh"] app.kpenalty = [ 1.0, 1.0, 1.0e30, 1.0, 1.0e30, 1.0e30, 1.0e30, 1.0, 1.0] # app.kpenalty = [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, "l0h", "re", "rh"] # app.optid_llsq = [ 0, 0, 1, 1, 0, 0, "l0h", "re", "rh"] app.add_argument(opt = None, type = "str", var_name = 'infile', opt_str = "input .xlsx file", desc = 'input file', defval = None, optional = False); app.add_argument(opt = "--mode", type = "str", var_name = 'mode', opt_str = "--mode=[EF|T|fit]", desc = 'task mode', defval = 'EF', optional = True); app.add_argument(opt = "--T_label", type = "str", var_name = 'T_label', opt_str = "--T_label=reg_exp", desc = 'regular expression for T data', defval = 'T.*\(.*K', optional = True); app.add_argument(opt = "--sigma_label", type = "str", var_name = 'sigma_label', opt_str = "--sigma_label=reg_exp", desc = 'regular expression for sigma data', defval = 'sigma', optional = True); app.add_argument(opt = "--S_label", type = "str", var_name = 'S_label', opt_str = "--S_label=reg_exp", desc = 'regular expression for S data', defval = 'S.*\(.*V.*K', optional = True); app.add_argument(opt = "--T0", type = "float", var_name = 'T0', opt_str = "--T0=val", desc = 'Temperature for mode=EF', defval = 300.0, optional = True); app.add_argument(opt = "--Eg", type = "float", var_name = 'Eg', opt_str = "--Eg=val", desc = 'band gap', defval = 1.1, optional = True); app.add_argument(opt = "--ED", type = "float", var_name = 'ED', opt_str = "--ED=val", desc = 'donar level from EV', defval = 1.0, optional = True); app.add_argument(opt = "--ND", type = "float", var_name = 'ND', opt_str = "--ND=val", desc = 'donar density in cm^-3', defval = 1.0e17, optional = True); app.add_argument(opt = "--EA", type = "float", var_name = 'EA', opt_str = "--EA=val", desc = 'acceptor level from EV', defval = 0.1, optional = True); app.add_argument(opt = "--NA", type = "float", var_name = 'NA', opt_str = "--NA=val", desc = 'acceptor density in cm^-3', defval = 0.0, optional = True); app.add_argument(opt = "--me", type = "float", var_name = 'me', opt_str = "--me=val", desc = 'electron effective mass in me0', defval = 1.0, optional = True); app.add_argument(opt = "--mh", type = "float", var_name = 'mh', opt_str = "--mh=val", desc = 'hole effective mass in me0', defval = 1.0, optional = True); app.add_argument(opt = "--l0e", type = "float", var_name = 'l0e', opt_str = "--l0e=val", desc = 'electron mean free path prefactor in m', defval = 1.0e-10, optional = True); app.add_argument(opt = "--l0h", type = "float", var_name = 'l0h', opt_str = "--l0h=val", desc = 'hole mean free path prefactor in m', defval = 1.0e-10, optional = True); app.add_argument(opt = "--rfac_e", type = "float", var_name = 'rfac_e', opt_str = "--rfac_e=val", desc = 'electron scattering factor', defval = 0.5, optional = True); app.add_argument(opt = "--rfac_h", type = "float", var_name = 'rfac_h', opt_str = "--rfac_h=val", desc = 'hole scattering factor', defval = 0.5, optional = True); # app.add_argument(opt = "--Tmin", type = "float", var_name = 'Tmin', opt_str = "--Tmin=val", desc = 'calculation T minimum', # defval = 0.0, optional = True); # app.add_argument(opt = "--Tmax", type = "float", var_name = 'Tmax', opt_str = "--Tmax=val", desc = 'calculation T maximum', # defval = 0.0, optional = True); app.add_argument(opt = "--method", type = "str", var_name = 'method', opt_str = "--method=[nelder-mead]", desc = 'optimization algorism', defval = 'nelder-mead', optional = True); app.add_argument(opt = "--wS", type = "float", var_name = 'wS', opt_str = "--wS=value", desc = 'weight of S w.r.t. that of ln(sigma)', defval = 1000, optional = True); app.add_argument(opt = "--nmaxiter", type = "int", var_name = 'nmaxiter', opt_str = "--nmaxiter=int(>1)", desc = 'maximum interation number for optimization', defval = 1000, optional = True); app.add_argument(opt = "--tol", type = "float", var_name = 'tol', opt_str = "--tol=val", desc = 'eps for optimization', defval = 1.0e-5, optional = True); # Temperature range cparams.Tmin = 100 # K cparams.Tmax = 1000 # Calclation range cparams.Tcalmin = 50 # K cparams.Tcalmax = 1000 cparams.Tstep = 10 cparams.nT = int((cparams.Tmax - cparams.Tmin) / cparams.Tstep + 1) cparams.h_diff = 1.0e-3 cparams.print_interval = 5 cparams.plot_interval = 1 #============================= # Graph configuration #============================= config.figsize = [12, 10] config.fontsize = 12 config.legend_fontsize = 8 config.graphupdateinterval = 10 return app, cparams, config #============================= # Treat argments #============================= def update_vars(app, cparams, config, apply_default = True): args_opt, args_idx, args_vars = app.read_args(vars = cparams, check_allowed_args = True, apply_default = apply_default) if args_opt is None: error_no = args_idx error_message = args_vars app.terminate("\n\n" + "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" + f"! {error_message}\n" + "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!", usage = app.usage) def read_fitting_parameter(path, cparams): fit = tkFit_mxy(tol = cparams.tol, nmaxiter = cparams.nmaxiter, print_interval = cparams.print_interval, plot_interval = cparams.plot_interval) # fit.read_parameters(cparams.parameterfile, heading = "\nRead parameters from [{cparams.parameterfile}]", # section = None, params = cparams, # keys = app.varname) cparams.read_parameters(path = cparams.parameterfile, section = None, AddSection = 0, terminator = ':', IsPrint = False, follow_vartype = True) return cparams, fit def print_simulation_parameters(cparams): print("Parameters:") print(f" Eg: {cparams.Eg} eV") print(f" ED: {cparams.ED} eV") print(f" ND: {cparams.ND} cm^-3") print(f" EA: {cparams.EA} eV") print(f" NA: {cparams.NA} cm^-3") print(f" me: {cparams.me} me") print(f" mh: {cparams.mh} me") print(f" electron:") print(f" rfac={cparams.rfac_e}") print(f" l0={cparams.l0e} m") print(f" hole:") print(f" rfac={cparams.rfac_h}") print(f" l0={cparams.l0h} m") def set_fitting_parameters(app, cparams): dos = app.dos mu_e = app.mu_e mu_h = app.mu_h dos.meeff = cparams.me dos.mheff = cparams.mh dos.EC = cparams.Eg dos.EA = cparams.EA dos.NA = cparams.NA dos.ED = cparams.ED dos.ND = cparams.ND dos.EF0 = 0.0 # initial EF to find EF mu_e.meff = dos.meeff # in me mu_e.rfac = cparams.rfac_e # tau = (meff/2)^0.5 * l0(T) * E^(r-0.5) mu_e.l0 = cparams.l0e # m mu_h.meff = dos.mheff # in me mu_h.rfac = cparams.rfac_h # tau = (meff/2)^0.5 * l0(T) * E^(r-0.5) mu_h.l0 = cparams.l0h # m # mu_h.l0 = mu_e.l0 * dos.mheff / dos.meeff klatt = 0.5 # W/m/K return dos, mu_e, mu_h, klatt def cal_EF(xEF, cparams, dos, mu_e, mu_h, klatt): yNe = [] yNh = [] ymue = [] ymuh = [] ysigmae = [] ysigmah = [] ySe = [] ySh = [] yStot = [] yRHe = [] yRHh = [] yRHtot = [] ysigma = [] yNHall = [] print() print(f"T={cparams.T0} K") labels = ["EF(eV)", "sigma,e", "sigma,h", "sigma,tot", "n,e", "n,h", "n,s", "mu,e", "mu,h", "S,e", "S,h", "S,tot", "RH,e", "RH,h", "RH,tot"] print_line(labels, format = '{:^12}') for EF in xEF: kBT = kB * cparams.T0 kT = e / kBT kBTe = kBT / e dos.NC = meff2NC_FEA(cparams.me, cparams.T0) dos.DC0 = meff2DC0_FEA(cparams.me, cparams.T0) dos.NV = meff2NC_FEA(cparams.mh, cparams.T0) dos.DV0 = meff2DC0_FEA(cparams.mh, cparams.T0) xe = (EF - dos.EC) * kT xh = (dos.EV - EF) * kT sigmae, ne, mue, tau_avge, Se, kappae, kappa_tote, Le, PFe, ZTe, infe = \ dos.cal_transport_S(xe, cparams.T0, dos.meeff, dos.NC, mu_e, klatt = klatt, validate_error_str = None, charge = None) sigmah, nh, muh, tau_avgh, Sh, kappah, kappa_toth, Lh, PFh, ZTh, infh = \ dos.cal_transport_S(xh, cparams.T0, dos.mheff, dos.NV, mu_h, klatt = klatt, validate_error_str = None, charge = None) sigmaHe, nHe, muHe, FHe, RH0e, RHe, nHalle, muHalle = \ dos.cal_transport_properties_Hall(xe, cparams.T0, dos.meeff, dos.NC, mu_e) sigmaHh, nHh, muHh, FHh, RH0h, RHh, nHallh, muHallh = \ dos.cal_transport_properties_Hall(xh, cparams.T0, dos.mheff, dos.NV, mu_h) sigma = sigmae + sigmah # Stot = (Sh + sigmae / sigmah * Se) / (1.0 + sigmae / sigmah) Stot = dos.cal_S_2(Sh, sigmah, Se, sigmae) # RHtot = (RHe * sigmae * sigmae + RHh * sigmah * sigmah) / sigma / sigma RHtot = dos.cal_RH(muHh, nHallh, muHe, nHalle) NeHall = 1.0 / e / RHtot yNe.append(ne) yNh.append(nh) ymue.append(mue) ymuh.append(muh) ysigmae.append(sigmae) ysigmah.append(sigmah) ySe.append(Se) ySh.append(Sh) yRHe.append(RHe) yRHh.append(RHh) ysigma.append(sigma) yStot.append(Stot) yRHtot.append(RHtot) yNHall.append(NeHall) print_line([EF, sigmae, sigmah, ne, nh, NeHall, mue, muh, Se, Sh, Stot, RHe, RHh, RHtot], format = '{:>12.4g}') # print(f" e(S) : s={sigmae:8.4g} n={ne:8.4g} mu={mue:8.4g} S={Se:8.4g}") # print(f" h(S) : s={sigmah:8.4g} n={nh:8.4g} mu={muh:8.4g} S={Sh:8.4g}") # print(f" e(Hall): s={sigmaHe:8.4g} n={nHe:8.4g} mu={muHe:8.4g} FH={FHe:8.4g} RH0={RH0e:8.4g} RH={RHe:8.4g} nHall={nHalle:8.4g} muHall={muHalle:8.4g}") # print(f" h(Hall): s={sigmaHh:8.4g} n={nHh:8.4g} mu={muHh:8.4g} FH={FHh:8.4g} RH0={RH0h:8.4g} RH={RHh:8.4g} nHall={nHallh:8.4g} muHall={muHallh:8.4g}") # print(f" RH: e={RHe:8.4g} h={RHh:8.4g} tot={RHtot:8.4g}") # print(f" S : e={Se:8.4g} h={Sh:8.4g} tot={Stot:8.4g}") yNHall = [abs(N) for N in yNHall] return labels, xEF, ysigmae, ysigmah, ysigma, yNe, yNh, yNHall, ymue, ymuh, ySe, ySh, yStot, yRHe, yRHh, yRHtot def sim_EF(app): config = app.config cparams = app.cparams cparams.outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-simulate.xlsx"]) cparams.EFmin = -0.2 cparams.EFmax = cparams.Eg + 0.2 cparams.nEF = 201 EFstep = (cparams.EFmax - cparams.EFmin) / (cparams.nEF - 1) print("") print("Simulate mu-EF dependence and compare with Hall data") cparams, fit = read_fitting_parameters(cparams.parameterfile, cparams) # cparams.print_parameters() print(f"module_dir : {config.module_dir}") print(f"module_name: {config.module_name}") print("mode : ", cparams.mode) print("infile : ", cparams.infile) print("outfile: ", cparams.outfile) print() print_simulation_parameters(cparams) print("") print(f"T: {cparams.T0} K") print(f"EF range: {cparams.EFmin:8.4g} - {cparams.EFmax:8.4g} at {EFstep:8.4g}, {cparams.nEF} points") dos, mu_e, mu_h, klatt = set_fitting_parameters(app, cparams) xEF = [cparams.EFmin + i * EFstep for i in range(cparams.nEF)] labels, xEF, ysigmae, ysigmah, ysigma, yNe, yNh, yNHall, ymue, ymuh, ySe, ySh, yStot, yRHe, yRHh, yRHtot = \ cal_EF(xEF, cparams, dos, mu_e, mu_h, klatt) outxlsxfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-EF.xlsx"]) print("") print(f"Save to [{outxlsxfile}]") tkVariousData().to_excel(outfile = outxlsxfile, labels = labels, data_list = [xEF, ysigmae, ysigmah, ysigma, yNe, yNh, yNHall, ymue, ymuh, ySe, ySh, yStot, yRHe, yRHh, yRHtot]) #======================================== # plot #======================================== fig, axes = plt.subplots(2, 2, figsize = config.figsize) axes = [ax for axes1 in axes for ax in axes1] plot_event = tkPlotEvent(plt) for ax in axes: ax.tick_params(labelsize = config.fontsize) axa = axes[0] axb = axa.twinx() p1, = axa.plot(xEF, yNe, label = 'Ne', linestyle = '-', color = 'red') p2, = axa.plot(xEF, yNh, label = 'Nh', linestyle = 'dashed', color = 'black') p1, = axa.plot(xEF, yNHall, label = 'N(Hall)', linestyle = '-', color = 'black') p3, = axb.plot(xEF, ymue, label = 'mue', linestyle = '-', color = 'blue') p4, = axb.plot(xEF, ymue, label = 'muh', linestyle = 'dashed', color = 'blue') axa.set_yscale('log') axa.set_xlabel('EF (eV)', fontsize = config.fontsize) axa.set_ylabel('N', fontsize = config.fontsize) axb.set_ylabel('mu', fontsize = config.fontsize) axa.legend(handles = [p1, p2, p3, p4], fontsize = config.legend_fontsize) axa = axes[1] # axb = axa.twinx() axa.plot(xEF, ysigmae, label = 'sigmae', linestyle = '-', color = 'red') axa.plot(xEF, ysigmah, label = 'sigmah', linestyle = 'dashed', color = 'blue') axa.plot(xEF, ysigma, label = 'sigma(tot)', linestyle = '-', color = 'black') axa.set_yscale('log') axa.set_xlabel('EF (eV)', fontsize = config.fontsize) axa.set_ylabel('sigma', fontsize = config.fontsize) axa.legend(fontsize = config.legend_fontsize) axa = axes[2] # axb = axa.twinx() axa.plot(xEF, ySe, label = 'Se', linestyle = '-', color = 'red') axa.plot(xEF, ySh, label = 'Sh', linestyle = 'dashed', color = 'blue') axa.plot(xEF, yStot, label = 'S(tot)', linestyle = '-', color = 'black') axa.set_xlabel('EF (eV)', fontsize = config.fontsize) axa.set_ylabel('S', fontsize = config.fontsize) axa.legend(fontsize = config.legend_fontsize) axa = axes[3] # axb = axa.twinx() axa.plot(xEF, yRHe, label = 'RHe', linestyle = '-', color = 'red') axa.plot(xEF, yRHh, label = 'RHh', linestyle = 'dashed', color = 'blue') axa.plot(xEF, yRHtot, label = 'RH(tot)', linestyle = '-', color = 'black') axa.set_xlabel('EF (eV)', fontsize = config.fontsize) axa.set_ylabel('RH', fontsize = config.fontsize) axa.legend(fontsize = config.legend_fontsize) plot_event.register_event(fig) plt.tight_layout() plt.pause(0.001) def read_data(infile, module_dir, module_name, app, cparams): fit = tkFit_mxy(tol = cparams.tol, nmaxiter = cparams.nmaxiter, print_interval = cparams.print_interval, plot_interval = cparams.plot_interval) print(f"Read data from [{infile}]") print(f" Module: {module_dir}/{module_name}.py") fit.read_data(infile = infile, module_dir = module_dir, module_name = module_name, imodule = 0, iregion = 0, x_indexes = [cparams.T_label], y_indexes = [cparams.sigma_label, cparams.S_label], app = app, cparams = cparams, is_print = False) labels_input = [fit.xlabels[0], fit.ylabels[0], fit.ylabels[1]] print("labels_input:", labels_input) print_data(labels = labels_input, data_list = [fit.x_list[0], fit.y_list[0], fit.y_list[1]], label_format = '{:^15}', data_format = '{:>15.4g}', header = "Input data:", print_level = 0) fit.w_list = [1.0, cparams.wS] fit.y_convert = ['', ''] # fit.y_convert = ['log', ''] # fit.y_list[0] = [log(y) for y in fit.y_list[0]] return fit def cal_T(xT, sigma_obs, S_obs, cparams, dos, mu_e, mu_h, klatt): EFmin = dos.EV - 0.2 EFmax = dos.EC + 0.2 Estep = 0.01 nrange = 6.0 yEF = [] ydQ = [] yNe = [] yNh = [] ymue = [] ymuh = [] ysigmae = [] ysigmah = [] ySe = [] ySh = [] yStot = [] yRHe = [] yRHh = [] yRHtot = [] ysigma = [] yNHall = [] print() labels = ["T(K)", "EF(eV)", "dQ(e)", "sigma,e", "sigma,h", "sigma,tot", "sigma,obs", "n,e", "n,h", "n,s", "mu,e", "mu,h", "S,e", "S,h", "S,tot", "S,obs", "RH,e", "RH,h", "RH,tot"] print_line(labels, format = '{:^12}') for i, _T in enumerate(xT): kBT = kB * _T kT = e / kBT kBTe = kB * _T / e dE = nrange * kBTe E0 = dos.EV - dE E1 = dos.EC + dE dos.NC = meff2NC_FEA(cparams.me, _T) dos.DC0 = meff2DC0_FEA(cparams.me, _T) dos.NV = meff2NC_FEA(cparams.mh, _T) dos.DV0 = meff2DC0_FEA(cparams.mh, _T) EF, dQ = dos.FindEF_semi(_T, EFmin, EFmax, E0, dos.EV, dos.EC, E1, Estep, print_level = 1) xe = (EF - dos.EC) * kT xh = (dos.EV - EF) * kT sigmae, ne, mue, tau_avge, Se, kappae, kappa_tote, Le, PFe, ZTe, infe = \ dos.cal_transport_S(xe, _T, dos.meeff, dos.NC, mu_e, klatt = klatt, validate_error_str = None, charge = None) sigmah, nh, muh, tau_avgh, Sh, kappah, kappa_toth, Lh, PFh, ZTh, infh = \ dos.cal_transport_S(xh, _T, dos.mheff, dos.NV, mu_h, klatt = klatt, validate_error_str = None, charge = None) sigmaHe, nHe, muHe, FHe, RH0e, RHe, nHalle, muHalle = \ dos.cal_transport_properties_Hall(xe, _T, dos.meeff, dos.NC, mu_e) sigmaHh, nHh, muHh, FHh, RH0h, RHh, nHallh, muHallh = \ dos.cal_transport_properties_Hall(xh, _T, dos.mheff, dos.NV, mu_h) sigma = sigmae + sigmah Stot = dos.cal_S_2(Sh, sigmah, Se, sigmae) RHtot = dos.cal_RH(muHh, nHallh, muHe, nHalle) NeHall = 1.0 / e / RHtot yEF.append(EF) ydQ.append(dQ) yNe.append(ne) yNh.append(nh) ymue.append(mue) ymuh.append(muh) ysigmae.append(sigmae) ysigmah.append(sigmah) ySe.append(Se) ySh.append(Sh) yRHe.append(RHe) yRHh.append(RHh) ysigma.append(sigma) yStot.append(Stot) yRHtot.append(RHtot) yNHall.append(NeHall) print_line([_T, EF, sigmae, sigmah, sigma, sigma_obs[i], ne, nh, NeHall, mue, muh, Se, Sh, Stot, S_obs[i], RHe, RHh, RHtot], format = '{:>12.4g}') yNHall = [abs(N) for N in yNHall] return labels, xT, yEF, ydQ, ysigmae, ysigmah, ysigma, sigma_obs, yNe, yNh, yNHall, ymue, ymuh, ySe, ySh, yStot, S_obs, yRHe, yRHh, yRHtot def plot_T(fig, axes, config, cparams, T_obs, yEF, ydQ, ysigmae, ysigmah, ysigma, sigma_obs, yNe, yNh, yNHall, ymue, ymuh, ySe, ySh, yStot, S_obs, yRHe, yRHh, yRHtot): plot_event = tkPlotEvent(plt) for ax in axes: ax.tick_params(labelsize = config.fontsize) T = T_obs axa = axes[0] axb = axa.twinx() p1, = axa.plot(T, yNe, label = 'Ne', linestyle = '-', color = 'red', marker = 'o') p2, = axa.plot(T, yNh, label = 'Nh', linestyle = 'dashed', color = 'black', marker = '^') p3, = axa.plot(T, yNHall, label = 'N,s', linestyle = '-', color = 'black', marker = 's', markersize = 3.0, markerfacecolor = 'w') p4, = axb.plot(T, ymue, label = 'mue', linestyle = '-', color = 'blue', marker = '+') p5, = axb.plot(T, ymue, label = 'muh', linestyle = 'dashed', color = 'blue', marker = 'x') axa.set_yscale('log') axa.set_xlabel('T (K)', fontsize = config.fontsize) axa.set_ylabel('N', fontsize = config.fontsize) axb.set_ylabel('mu', fontsize = config.fontsize) axa.legend(handles = [p1, p2, p3, p4, p5], fontsize = config.legend_fontsize) axa = axes[1] # axb = axa.twinx() axa.plot(T, ysigmae, label = 'sigmae', linestyle = 'dashed', color = 'red', marker = '<', markersize = 5.0) axa.plot(T, ysigmah, label = 'sigmah', linestyle = 'dashed', color = 'blue', marker = '^', markersize = 5.0) axa.plot(T, ysigma, label = 'sigma(tot)', linestyle = 'dashed', color = 'darkgreen', marker = 's', markersize = 3.0, markerfacecolor = 'w') axa.plot(T_obs, sigma_obs, label = 'sigma(meas)', linestyle = '-', color = 'black', marker = 'o', markersize = 8.0) axa.set_yscale('log') axa.set_xlabel('T (K)', fontsize = config.fontsize) axa.set_ylabel('sigma', fontsize = config.fontsize) axa.legend(fontsize = config.legend_fontsize) axa = axes[2] # axb = axa.twinx() axa.plot(T, ySe, label = 'Se', linestyle = 'dashed', color = 'red', marker = '<', markersize = 5.0) axa.plot(T, ySh, label = 'Sh', linestyle = 'dashed', color = 'blue', marker = '^', markersize = 5.0) axa.plot(T, yStot, label = 'S(tot)', linestyle = 'dashed', color = 'darkgreen', marker = 's', markersize = 3.0, markerfacecolor = 'w') axa.plot(T_obs, S_obs, label = 'S(meas)', linestyle = '-', color = 'black', marker = 'o', markersize = 8.0) axa.set_xlabel('T (K)', fontsize = config.fontsize) axa.set_ylabel('S', fontsize = config.fontsize) axa.legend(fontsize = config.legend_fontsize) axa = axes[3] # axb = axa.twinx() axa.plot(T, yRHe, label = 'RHe', linestyle = '-', color = 'red', marker = 'o') axa.plot(T, yRHh, label = 'RHh', linestyle = 'dashed', color = 'blue', marker = '^') axa.plot(T, yRHtot, label = 'RH(tot)', linestyle = '-', color = 'black', marker = 's', markersize = 3.0, markerfacecolor = 'w') axa.set_xlabel('T (K)', fontsize = config.fontsize) axa.set_ylabel('RH', fontsize = config.fontsize) axa.legend(fontsize = config.legend_fontsize) axa = axes[4] axb = axa.twinx() p4a, = axa.plot(T, yEF, label = 'EF', linestyle = '-', color = 'black', marker = 'o') p4b, = axb.plot(T, ydQ, label = 'dQ', linestyle = 'dashed', color = 'blue', linewidth = 0.5, marker = '^', markersize = 1.5) axa.axhline(0.0, linestyle = 'dashed', color = 'red', linewidth = 0.5) axa.axhline(cparams.Eg, linestyle = 'dashed', color = 'red', linewidth = 0.5) axa.set_xlabel('T (K)', fontsize = config.fontsize) axa.set_ylabel('EF (eV)', fontsize = config.fontsize) axb.set_ylabel(r'$\Delta$$Q$ (e)', fontsize = config.fontsize) axa.legend(handles = [p4a, p4b], fontsize = config.legend_fontsize) plot_event.register_event(fig) plt.tight_layout() plt.pause(0.001) def sim_T(app): config = app.config cparams = app.cparams cparams.outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-simulate.xlsx"]) print("") print("Simulate mu-T dependence and compare with Hall data") cparams, fit = read_fitting_parameter(cparams.parameterfile, cparams) # cparams.print_parameters() print(f"module_dir : {config.module_dir}") print(f"module_name: {config.module_name}") print("mode : ", cparams.mode) print("infile : ", cparams.infile) print(" T label : ", cparams.T_label) print(" sigma label: ", cparams.sigma_label) print(" S label : ", cparams.S_label) # print(f"T range for simulation: {cparams.Tcalmin} - {cparams.Tcalmax} K") print("outfile: ", cparams.outfile) print() print_simulation_parameter(cparams) print("") fit = read_data(cparams.infile, config.module_dir, config.module_name, app, cparams) T_label, sigma_label, S_label = [fit.xlabels[0], fit.ylabels[0], fit.ylabels[1]] T_obs, sigma_obs, S_obs = [fit.x_list[0], fit.y_list[0], fit.y_list[1]] muw = [] New = [] for _T, _sigma, _S in zip(T_obs, sigma_obs, S_obs): _mu, _ne, sign, carriertype = weighted_mobility(_sigma / 0.01, _S * 1.0e-6, _T, unit = 'cm') muw.append(_mu) New.append(_ne) print_data(labels = [T_label, sigma_label, S_label, "mu_w (cm2/Vs)", "Ne_w (cm-3)"], data_list = [T_obs, sigma_obs, S_obs, muw, New], label_format = '{:^15}', data_format = '{:>15.4g}', header = "Input data:", print_level = 0) dos, mu_e, mu_h, klatt = set_fitting_parameters(app, cparams) print() labels = ["T(K)", "NV", "NC", "DV0", "DC0"] print_line(labels, format = '{:^12}') for i, _T in enumerate(T_obs): dos.NC = meff2NC_FEA(cparams.me, _T) dos.DC0 = meff2DC0_FEA(cparams.me, _T) dos.NV = meff2NC_FEA(cparams.mh, _T) dos.DV0 = meff2DC0_FEA(cparams.mh, _T) print_line([_T, dos.NV, dos.NC, dos.DV0, dos.DC0], format = '{:^12.4g}') labels, xT, yEF, ydQ, ysigmae, ysigmah, ysigma, sigma_obs, yNe, yNh, yNHall, ymue, ymuh, ySe, ySh, yStot, S_obs, yRHe, yRHh, yRHtot = \ cal_T(T_obs, sigma_obs, S_obs, cparams, dos, mu_e, mu_h, klatt) outxlsxfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-T.xlsx"]) print("") print(f"Save to [{outxlsxfile}]") tkVariousData().to_excel(outfile = outxlsxfile, labels = labels, data_list = [T_obs, yEF, ydQ, ysigmae, ysigmah, ysigma, sigma_obs, yNe, yNh, yNHall, ymue, ymuh, ySe, ySh, yStot, S_obs, yRHe, yRHh, yRHtot]) #======================================== # plot #======================================== fig, axes = plt.subplots(2, 3, figsize = config.figsize) axes = [ax for axes1 in axes for ax in axes1] plot_T(fig, axes, config, cparams, T_obs, yEF, ydQ, ysigmae, ysigmah, ysigma, sigma_obs, yNe, yNh, yNHall, ymue, ymuh, ySe, ySh, yStot, S_obs, yRHe, yRHh, yRHtot) def cal_S_sigma(x_list, dos, mu_e, mu_h, dEF, Estep, nrange, pk): dos.EC, dos.ED, dos.ND, dos.EA, dos.NA, mu_e.l0, mu_h.l0, mu_e.rfac, mu_h.rfac = [*pk] # dos.ND = exp(dos.ND) # dos.NA = exp(dos.NA) # mu_e.l0 = exp(mu_e.l0) # mu_h.l0 = mu_e.l0 * dos.mheff / dos.meeff T = x_list[0] kBT = kB * T ekT = e / kBT EFmin = -dEF EFmax = dos.EC + dEF dE = nrange * kBT / e E0 = dos.EV - dE E1 = dos.EC + dE dos.NC = meff2NC_FEA(dos.meeff, T) dos.DC0 = meff2DC0_FEA(dos.meeff, T) dos.NV = meff2NC_FEA(dos.mheff, T) dos.DV0 = meff2DC0_FEA(dos.mheff, T) EF, dQ = dos.FindEF_semi(T, EFmin, EFmax, E0, dos.EV, dos.EC, E1, Estep, print_level = 1) xe = (EF - dos.EC) * ekT xh = (dos.EV - EF) * ekT sigmae, ne, mue, tau_avge, Se, kappae, kappa_tote, Le, PFe, ZTe, infe = \ dos.cal_transport_S(xe, T, dos.meeff, dos.NC, mu_e, klatt = 0.5, validate_error_str = None, charge = None) sigmah, nh, muh, tau_avgh, Sh, kappah, kappa_toth, Lh, PFh, ZTh, infh = \ dos.cal_transport_S(xh, T, dos.mheff, dos.NV, mu_h, klatt = 0.5, validate_error_str = None, charge = None) sigma = sigmae + sigmah # Stot = (Sh + sigmae / sigmah * Se) / (1.0 + sigmae / sigmah) Stot = dos.cal_S_2(Sh, sigmah, Se, sigmae) print(f"{ne=:8.3g} {nh=:8.3g} {Se=:8.3g} {Sh=:8.3g} {mue=:8.3g} {mue=:8.3g} {sigmae=:8.3g} {sigmah=:8.3g}") # sigma = log(sigma) return [sigma, Stot] def fit_T(app): config = app.config cparams = app.cparams cparams.outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-simulate.xlsx"]) print("") print("Fitting to mu-T Hall data by nonlinear least-squares method") cparams, fit = read_fitting_parameter(cparams.parameterfile, cparams) # cparams.print_parameters() print(f"module_dir : {config.module_dir}") print(f"module_name: {config.module_name}") print("mode : ", cparams.mode) print("infile : ", cparams.infile) print(" T label : ", cparams.T_label) print(" sigma label: ", cparams.sigma_label) print(" S label : ", cparams.S_label) print("outfile: ", cparams.outfile) print("parameter file : ", cparams.parameterfile) # Load modules print("") fit = read_data(cparams.infile, config.module_dir, config.module_name, app, cparams) dos, mu_e, mu_h, klatt = set_fitting_parameters(app, cparams) if dos.ND <= 100.0: dos.ND = 100.0 if dos.NA <= 100.0: dos.NA = 100.0 fit.varname = app.varname fit.unit = app.unit fit.pk = [dos.EC, dos.ED, dos.ND, dos.EA, dos.NA, mu_e.l0, mu_h.l0, mu_e.rfac, mu_h.rfac] # fit.pk = [dos.EC, dos.ED, log(dos.ND), dos.EA, log(dos.NA), log(mu_e.l0)] fit.pk_convert = app.pk_convert fit.optid = app.optid fit.kmin = app.kmin fit.kmax = app.kmax fit.kpenalty = app.kpenalty # fit.optid_llsq = [ 0, 0, 1, 1, 0, 0] nvars = len(fit.pk) method = cparams.method optid_original = fit.optid.copy() print(f"Read fitting parameters from {cparams.parameterfile}") fit.read_parameters(cparams.parameterfile, section = "S-T", keys = fit.varname, values = fit.pk, optid = fit.optid, kmin = fit.kmin, kmax = fit.kmax, kpenalty = fit.kpenalty) print() fit.print_variables(heading = "Fitting parameters:") print(f"Constants") print(f" me: {cparams.me} me") print(f" mh: {cparams.mh} me") print(f" electron:") print(f" rfac={cparams.rfac_e}") print(f" l0={cparams.l0e} m") print(f" hole:") print(f" rfac={cparams.rfac_h}") print(f" l0={cparams.l0h} m") print(f"Fitting configuration") print(f" T range for fitting : {cparams.Tmin} - {cparams.Tmax} K") print(f" T range for simulation: {cparams.Tcalmin} - {cparams.Tcalmax} K") print(f" method : {cparams.method}") print(f" wS : {cparams.wS}") print(f" tol : {cparams.tol}") print(f" nmaxiter: {cparams.nmaxiter}") print(f" h_diff : {cparams.h_diff}") Estep = 0.01 dEF = 0.2 nrange = 6.0 fit.func = lambda x_list, *pk: cal_S_sigma(x_list, dos, mu_e, mu_h, dEF, Estep, nrange, *pk) # print("") # fit.print_attributes() yini_list = fit.cal_ylist(fit.pk) optpk = fit.extract_parameters() fini = fit.minimize_func(optpk, x_list = fit.x_list, y_list = fit.y_list, w_list = fit.w_list) print("") print(f"Initial function: fmin={fini:10.4g}") print_data(labels = ["sigma", "sigma,cal", fit.ylabels[1], "S(cal)"], # print_data(labels = ["log(s)", "log(s),cal", fit.ylabels[1], "S(cal)"], data_list = [fit.x_list[0], fit.y_list[0], yini_list[0], fit.y_list[1], yini_list[1]], label_format = '{:^15}', data_format = '{:>15.4g}', header = "Initial data:", print_level = 0) #======================================== # plot #======================================== fig, axes = plt.subplots(len(fit.y_list), 3, figsize = config.figsize) # axes = [axes] # axes = axes.flatten() title = app.replace_path(cparams.infile, template = "{filename}") axes[0][0].set_title(title, fontsize = config.fontsize) fit.initial_plot(data_axes = [axes[0][0], axes[1][0]], error_axis = axes[0][2], yini_list = yini_list, label_ini = 'initial', fmin = fini, plt = plt, fig = fig, fontsize = config.fontsize) fit.plot_event.remove('error') #======================================== # Optimize #======================================== print("") print("Optimize:") fit.print_variables(heading = "Variables") pfin, ffin, success, res = fit.minimize(method) if success: print(f"Converged at iteration: {res.nit}") else: print(f"Function did not converge") #======================================== # Final result #======================================== dos.EC, dos.ED, dos.ND, dos.EA, dos.NA, mu_e.l0, mu_h.l0, mu_e.rfac, mu_h.rfac = [*pfin] fit.pk = pfin.copy() # dos.ND = exp(dos.ND) # dos.NA = exp(dos.NA) # mu_e.l0 = exp(mu_e.l0) print("Final parameters") for i in range(nvars): print(f" {fit.varname[i]:10}: {pfin[i]:10.4g} {fit.unit[i]}") print(f" f={ffin:12.6g}") print("") print(f"Save parameters to {cparams.parameterfile}") print("fit.pk=", fit.pk) fit.save_parameters(cparams.parameterfile, heading = "\nSave parameters:", section = "S-T", keys = fit.varname, values = fit.pk, optid = fit.optid, kmin = fit.kmin, kmax = fit.kmax, kpenalty = fit.kpenalty) cparams.save_parameters(cparams.parameterfile, section = 'Fitting', keys = ["Tmin", "Tmax", "print_interval", "plot_interval", "h_diff", "nmaxiter", "tol", "wS", "T_label", "sigma_label", "S_label", "me", "mh"]) yfin_list = fit.cal_ylist(pfin) print("") print_data(labels = [fit.xlabels[0], fit.ylabels[0], "sigma(cal)", fit.ylabels[1], "S(cal)"], data_list = [fit.x_list[0], fit.y_list[0], yfin_list[0], fit.y_list[1], yfin_list[1]], label_format = '{:^15}', data_format = '{:>15.4g}', header = "Final data:", print_level = 0) # fit.print_scores(heading = "\nScores between y(input) and y(fit)", y1 = fit.y, y2 = yfin) # Calculate T dep labels, xT, yEF, ydQ, ysigmae, ysigmah, ysigma, sigma_obs, yNe, yNh, yNHall, ymue, ymuh, ySe, ySh, yStot, S_obs, yRHe, yRHh, yRHtot = \ cal_T(fit.x_list[0], fit.y_list[0], fit.y_list[1], cparams, dos, mu_e, mu_h, klatt) outxlsxfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-T-final.xlsx"]) print("") print(f"Save to [{outxlsxfile}]") tkVariousData().to_excel(outfile = outxlsxfile, labels = labels, data_list = [xT, yEF, ydQ, ysigmae, ysigmah, ysigma, sigma_obs, yNe, yNh, yNHall, ymue, ymuh, ySe, ySh, yStot, S_obs, yRHe, yRHh, yRHtot]) #============================= # グラフの表示 #============================= print("") print("Plot optimized") fit.finalize_plot(yfin_list, iter = res.nit, fmin = ffin) fit.layout() # plt.tight_layout() plt.pause(0.01) def main(): print() app, cparams, config = initialize() update_vars(app, cparams, config, apply_default = True) cparams.logfile = app.replace_path(cparams.infile) print(f"Open logfile [{cparams.logfile}]") app.redirect(targets = ["stdout", cparams.logfile], mode = 'w') # app.redirect(targets = ["stdout"], mode = 'a', # redirect_traceback = True, output_traceback = 'stdout', display_type_traceback = 'colored') # app.redirect_exception() cparams.parameterfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}.in"]) cparams.parameterbkfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-back.in"]) print("") print("infile : {}".format(cparams.infile)) print("parameter file : {}".format(cparams.parameterfile)) print("parameter backup file: {}".format(cparams.parameterbkfile)) print("mode : {}".format(cparams.mode)) print("") print("Read [{}]".format(cparams.parameterfile)) cparams.read_parameters(cparams.parameterfile, section = "Preferences", ignore_keys = ["logfile", "outfile", "parameterfile", "parameterbkfile"]) cparams.read_parameters(cparams.parameterfile, section = "Fitting", ignore_keys = ["logfile", "outfile", "parameterfile", "parameterbkfile"]) cparams.read_parameters(cparams.parameterfile, section = "S-T", terminator = ':', ignore_keys = ["logfile", "outfile", "parameterfile", "parameterbkfile"]) # check args again to use given parameters update_vars(app, cparams, config, apply_default = False) config.print_parameters(heading = "\nconfig:") # cparams.print_parameters(heading = "\ncparams:") if cparams.mode == 'T': sim_T(app) elif cparams.mode == 'EF': sim_EF(app) elif cparams.mode == 'fit': fit_T(app) else: app.terminate("Error in main: Invalide mode [{}]".format(cparams.mode), usage = app.usage, pause = True) app.terminate("", usage = app.usage, pause = True) if __name__ == "__main__": main()