import sys import numpy as np from numpy import sin, cos, tan, pi, exp, log, sqrt from scipy.optimize import minimize from matplotlib import pyplot as plt from tklib.tksci.tksci import log10, e, kB from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me, Ea_Arrhenius from tklib.tkutils import getarg, getintarg, getfloatarg, pconv_by_type from tklib.tkapplication import tkApplication from tklib.tkparams import tkParams from tklib.tkvariousdata import tkVariousData from tklib.tksci.tkFit import tkFit from tklib.tkgraphic.tkplotevent import tkPlotEvent from tklib.tktransport.tkTransport import integrate_Simpson_list, fe, fh, NC2meff_FEA from tklib.tktransport.tkmobility import tkMobility from tklib.tksci.tkmatrix import make_matrix1 from tklib.tksci.tkoptimize import mlsq_general """ 非縮退半導体のEFの温度依存性を二分法で計算。 """ usage_str = ''' "" "(i) usage: python {} mode infile_path(.xlsx|.csv|.in) method mu0 Eb Aop Eop a0 a1 a2 a3 a4".format(sys.argv[0]) " mode : init, sim, fit" " Eb : Grain boundary photential barrier height [eV]" " Eop: Optical phonon energy [eV]" '''[1:-1] def initialize(): #================================ # Global variables #================================ app = tkApplication(usage_str = usage_str) argv, narg = app.get_argv() app.mu = tkMobility() cparams = app.get_params() cparams.debug = 0 cparams.use_simple = 0 cparams.print_level = 0 # mode: 'init', 'fit', 'sim' cparams.mode = 'init' # cparams.mode = 'fit' # data file to fit (.xlsx) cparams.infile = None cparams.parameterfile = None cparams.parameterbkfile = None cparams.outcsvfile = None cparams.outcomponentcsvfile = None cparams.T_label = None cparams.n_label = None cparams.mu_label = None cparams.sigma_label = None mu = app.mu # For debug purpose. 1 will use 3-terms polynomial for the inverse of mobility mu.debug = cparams.debug mu.use_simple = cparams.use_simple # 移動度パラメータ # 係数は、移動度（緩和時間）の逆数に関する係数であることに注意 # mu.Tnorm = 300.0 # 移動度の基準温度 mu.mu0 = 10.0 # Tnormにおける電子移動度, cm2/Vs mu.Eb = 0.0 # ポテンシャル障壁高さ, eV mu.Eop = 0.0001 # 光学フォノン振動エネルギー, eV。光学フォノン散乱を無視するときは小さい値を入れる mu.Aop = 1.0e-10 # 光学フォノン散乱前指数因子。光学フォノン散乱を無視するときは非常に小さい値を入れる mu.ppolyi = [0.0, -1.5, 1.5, -1.0, 1.0] mu.apolyi = [1.0, 0.0, 0.0, 0.0, 0.0] mu.varkeys = ["Eb", "Aop", "Eop", "apolyi[0]", "ppolyi[0]", "apolyi[1]", "ppolyi[1]", "apolyi[2]", "ppolyi[2]", "apolyi[3]", "ppolyi[3]", "apolyi[4]", "ppolyi[4]"] mu.header = None mu.xT = None mu.yN = None mu.ymu = None mu.ys = None app.ymuini = None # Temperature range cparams.Tmin = 50 # K cparams.Tmax = 1000 cparams.Tstep = 10 cparams.nT = int((cparams.Tmax - cparams.Tmin) / cparams.Tstep + 1) # 最適化パラメータの初期値 # Eb以外のパラメータは、係数、固定変数の順番で組み合わせになっていないといけない mu.varname = ("Eb", "Aop", "Eop", "a0", "p0", "a1", "p1", "a2", "p2", "a3", "p3", "a4", "p4") mu.varunit = ["eV", "", "eV", "", " ", "", "", "", "", "", "", "", ""] mu.ai0 = mu.parameter_list() mu.optid = [ 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0] mu.id_llsq = [ 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0] # scipy.optimize.minimizeで使うアルゴリズム cparams.method = "nelder-mead" #method = 'cg' #method = 'bfgs' cparams.maxiter = 1000 cparams.tol = 1.0e-4 cparams.h_diff = 1.0e-5 cparams.print_interval = 1 cparams.plot_interval = 1 #============================= # Graph configuration #============================= cparams.fig = None cparams.colors = ['black', 'red', 'blue', 'green', 'orange', 'darkgreen', 'darkorange', 'navy', 'blue', 'darkgreen', 'darkorange', 'navy', 'slategray', 'hotpink', 'olive', 'chocolate', 'magenta', 'green', 'yellow', 'cyan'] cparams.linewidth = 2.0 cparams.figsize = [12, 8] cparams.fontsize = 14 cparams.legend_fontsize = 12 cparams.graphupdateinterval = 10 return app, cparams, mu #============================= # 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, mu): cparams.mode = getarg(1, cparams.mode) cparams.method = getarg(2, cparams.method) cparams.infile = getarg(3, cparams.infile) cparams.T_label = getarg(4, cparams.T_label) cparams.n_label = getarg(5, cparams.n_label) cparams.mu_label = getarg(6, cparams.mu_label) cparams.sigma_label = getarg(7, cparams.sigma_label) cparams.Tmin = getfloatarg(8, cparams.Tmin) cparams.Tmax = getfloatarg(9, cparams.Tmax) cparams.maxiter = getintarg(10, cparams.maxiter) cparams.tol = getfloatarg(11, cparams.tol) cparams.h_diff = getfloatarg(12, cparams.h_diff) if cparams.infile is None: app.terminate("Error in updatevars: infile must be specified", usage = usage) def plot_muT_decomposed(app, cparams, mu, ax, xT, ymu, ymucal = None, ymu_ingrain = None, ymuop = None, ymus = None, xlabel = "T (K)", ylabel = "$\mu$ (cm$^2$/Vs)", colors = None, markersize = 1.0, plot_event = None): ax.clear() ax.tick_params(labelsize = cparams.fontsize) color = colors[0] data = ax.plot(xT, ymu, label = '$\mu(obs)$', linestyle = '', color = color, marker = 'o', markersize = markersize, markerfacecolor = color, markeredgecolor = color) plot_event.add_data({"label": "mu_obs", "plot_type": "plot", "axis": ax, "data": data}) if ymucal: color = colors[1] data = ax.plot(xT, ymucal, label = '$\mu(cal)$', color = color, linewidth = 1.0, linestyle = '-') plot_event.add_data({"label": "mu_cal", "plot_type": "plot", "axis": ax, "data": data}) if ymu_ingrain: color = colors[0] data = ax.plot(xT, ymu_ingrain, label = '$\mu_{in-grain}$', linewidth = cparams.linewidth, linestyle = 'dashed', color = color) plot_event.add_data({"label": "mu_ingrain", "plot_type": "plot", "axis": ax, "data": data}) if ymuop: color = colors[1] data = ax.plot(xT, ymuop, label = '$\mu_{op}$', linewidth = cparams.linewidth, linestyle = 'dashed', color = color) plot_event.add_data({"label": "mu_op", "plot_type": "plot", "axis": ax, "data": data}) if ymus: for i in range(len(ymus)): color = colors[2 + i] label = "$\mu$($T^{" + f"{mu.ppolyi[i]}" + "})$" data = ax.plot(xT, ymus[i], label = label, linewidth = cparams.linewidth, linestyle = 'dashed', color = color, marker = 'o', markersize = 2.0) plot_event.add_data({"label": label, "plot_type": "plot", "axis": ax, "data": data}) ax.set_xlabel(xlabel, fontsize = cparams.fontsize) ax.set_ylabel(ylabel, fontsize = cparams.fontsize) ax.legend(fontsize = cparams.legend_fontsize) def plot_muT_weight(app, cparams, mu, ax, xT, ywmugb, ywmuop, ywmus, xlabel = "T (K)", ylabel = "Linear weight", colors = None, markersize = 1.0, plot_event = None): color = colors[0] wgb_data = ax.plot(xT, ywmugb, label = '$w_{GB}$', linewidth = cparams.linewidth, linestyle = 'dashed', color = color) color = colors[1] wop_data = ax.plot(xT, ywmuop, label = '$w_{op}$', linewidth = cparams.linewidth, linestyle = 'dashed', color = color) plot_event.add_data({"label": "w_GB", "plot_type": "plot", "axis": ax, "data": wgb_data}) plot_event.add_data({"label": "w_op", "plot_type": "plot", "axis": ax, "data": wop_data}) for i in range(len(ywmus)): color = colors[2 + i] label = "w($T^{" + f"{mu.ppolyi[i]}" + "})$" data = ax.plot(xT, ywmus[i], label = label, linewidth = cparams.linewidth, linestyle = 'dashed', color = color, marker = 'o', markersize = 2.0) plot_event.add_data({"label": label, "plot_type": "plot", "axis": ax, "data": data}) ax.set_xlabel(xlabel, fontsize = cparams.fontsize) ax.set_ylabel(ylabel, fontsize = cparams.fontsize) ax.legend(fontsize = cparams.legend_fontsize) def cal_mu(mu, T, *pk): aiorg = mu.parameter_list() # ainew = mu.recover_parameters(pk) mu.set_parameters(pk) _mu = mu.cal_mu(T) mu.set_parameters(aiorg) return _mu def fitting(app, cparams, mu, fit): ekB = e / kB print("") print("Fitting to mu-T Hall data by nonlinear least-squares method") print("infile : ", cparams.infile) print("T_label : ", cparams.T_label) print("mu_label : ", cparams.mu_label) print(f"T range for fitting : {cparams.Tmin} - {cparams.Tmax} K") print("method : ", cparams.method) print("maxiter : ", cparams.maxiter) print("tol : ", cparams.tol) print("h : ", cparams.h_diff) if '***' in cparams.method: app.terminate("\nError: Choose method", pause = True) if '***' in cparams.T_label: app.terminate("\nError: Choose T_label", pause = True) if '***' in cparams.mu_label: app.terminate("\nError: Choose mu_label", pause = True) print("") print(f"Read data from [{cparams.infile}]") fit.read_data(cparams.infile, xlabel = cparams.T_label, ylabel = cparams.mu_label, xmin = cparams.Tmin, xmax = cparams.Tmax, usage = usage) xT = fit.x ymu = fit.y # fit.read_parameters(cparams.parameterfile, section = "mu", # keys = fit.varname, values = fit.pk, optid = fit.optid) print("ndata=", fit.ndata) print(f"{fit.xlabel:10} {fit.ylabel:10}") for i in range(fit.ndata): print(f"{fit.x[i]:10.4g} {fit.y[i]:10.4g}") fit.varname = mu.varname fit.unit = mu.varunit fit.pk = mu.parameter_list() fit.optid = mu.optid fit.id_llsq = mu.id_llsq # fit.kmin = [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # fit.kmax = [ 1.0e100, 1.0e100, 1.0e100, 1.0e100, 1.0e100, 1.0e100] # fit.kpenalty = [ 1.0e4, 1.0e6, 1.0e6, 1.0, 1.0e2, 1.0e6] nvars = len(fit.pk) fit.func = lambda T, *pk: cal_mu(mu, T, *pk) optpk = fit.extract_parameters() fini = fit.minimize_func(optpk, x = fit.x, y = fit.y) print("") print("Initial parameters:") fit.print_variables() print(f" f={fini:12.4g}") # calculate initial values print("") print("Calculate mu-T from the initial parameters") mu.print_parameters() print(f"{'T(K)':8}\t{'mu(cm2/Vs)':12}\t{'mu(cm2/Vs)(ini)':12}") # app.ymuini = fit.cal_ylist(fit.pk) app.ymuini = [] for iT in range(len(xT)): T = xT[iT] muc = cal_mu(mu, T, *fit.pk) # muc = mu.cal_mu(T) app.ymuini.append(muc) print(f"{T:8.4g}\t{ymu[iT]:12.6g}\t{muc:12.6g}") xTinv = [] print("") print("{:8}\t{:8}\t{:14}\t{:14}".format('T', '1000/T', 'mu (cm2/Vs)', 'mu (cal)')) for i in range(len(xT)): xTinv.append(1000.0 / xT[i]) print(f"{xT[i]:8.4f}\t{xTinv[i]:8.4f}\t{ymu[i]:14.6g}\t{app.ymuini[i]:14.6g}") #============================= # グラフの表示 #============================= print("") print("plot") app.plt = plt # fig, axes = plt.subplots(1, 2, figsize = cparams.figsize) fig = plt.figure(figsize = cparams.figsize) app.fig = fig ax1 = fig.add_subplot(1, 3, 1) ax2 = fig.add_subplot(1, 3, 2) ax3 = fig.add_subplot(1, 3, 3) axes = [ax1, ax2, ax3] title = app.replace_path(cparams.infile, template = "{filename}") plt.title(title, fontsize = cparams.fontsize) # axes = [axes] # axes = axes.flatten() fit.initial_plot(data_axis = axes[0], error_axis = axes[1], yini = app.ymuini, label_ini = 'initial', fmin = fini, fig = fig, fontsize = cparams.fontsize) fit.plot_event.remove('error') #======================================== # Optimize #======================================== print("") print("Nonlinear least-squares fitting by ", cparams.method) print(" tol=", cparams.tol) fit.print_variables(heading = "Variables") pfin, ffin, success, res = fit.minimize(cparams.method) if success: print("") print(f"Converged at iteration: {res.nit}") else: print(f"Function did not converge") mu.ai0 = pfin mu.Eb, mu.Aop, mu.Eop, mu.apolyi[0], mu.ppolyi[0], mu.apolyi[1], mu.ppolyi[1], \ mu.apolyi[2], mu.ppolyi[2], mu.apolyi[3], mu.ppolyi[3], \ mu.apolyi[4], mu.ppolyi[4] = [*pfin] print("") print("Final parameters") mu.print_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("Save parameters to [{}]".format(cparams.parameterfile)) cparams.save_parameters(cparams.parameterfile, "Preferences") mu.save_parameters(cparams.parameterfile, optid = None, parameters_section = 'Parameters', otherparams = {"S2mu": ffin}, other_section = 'OtherParameters') #======================================== # Final result #======================================== print("") print("Calculate final result") ymufin = fit.cal_ylist(pfin) # print(f"{'i':4} {'T':10} {'mu':12} {'mu(ini)':12} {'mu(fin)':12}") # for i in range(len(fit.x)): # print(f"{i:4} {fit.x[i]:10.4g} {fit.y[i]:12.4g} {app.ymuini[i]:12.4g} {ymufin[i]:12.4g}") outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-mu-fit.xlsx"]) print("") print(f"Save results to [{outfile}]") fit.to_excel(outfile, ['T(K)', 'mu', 'mu(ini)', 'mu(fin)'], [fit.x, fit.y, app.ymuini, ymufin]) fit.print_scores(heading = "\nScores between y(input) and y(fit)", y1 = fit.y, y2 = ymufin) #============================= # 貢献度の計算 #============================= # ymufin = [] yKgb = [] ymu_ingrain = [] ymuop = [] ymus = make_matrix1(len(mu.ppolyi), type = 'list') print("") print("Components of mobility") print("{:6} {:10} {:10} {:10} {:10} {:10} {:10} {:10} {:10} {:10} {:10}" .format("T(K)", "mu,obs(cm2/Vs)", "mu,tot", "Kgb", "mu,in-grain", "mu,op", f"mu(T^{mu.ppolyi[0]})", f"mu(T^{mu.ppolyi[1]})", f"mu(T^{mu.ppolyi[2]})", f"mu(T^{mu.ppolyi[3]})", f"mu(T^{mu.ppolyi[4]})")) for iT in range(len(xT)): T = xT[iT] mu_tot, Kgb, mu_ingrain, muop, mus = mu.cal_mu_components(T) # print("mu=", mu_tot, mus[0], mus[1], mus[2]) # ymufin.append(mu_tot) yKgb.append(Kgb) ymu_ingrain.append(mu_ingrain) ymuop.append(muop) for i in range(len(mus)): ymus[i].append(mus[i]) print("{:6.3g} {:10.4g} {:10.4g} {:10.4g} {:10.4g} {:10.4g} {:10.4g} {:10.4g} {:10.4g} {:10.4g} {:10.4g}" .format(T, ymu[iT], mu_tot, Kgb, mu_ingrain, muop, mus[0], mus[1], mus[2], mus[3], mus[4])) ywmugb = [] ywmuop = [] ywmus = make_matrix1(len(mu.ppolyi), type = 'list') print("") print("{:6} {:10} {:10} {:10} {:10} {:10} {:10} {:10}" .format("T(K)", "w,gb", "w,op", f"w(T^{mu.ppolyi[0]})", f"w(T^{mu.ppolyi[1]})", f"w(T^{mu.ppolyi[2]})", f"w(T^{mu.ppolyi[3]})", f"w(T^{mu.ppolyi[4]})")) for iT in range(len(xT)): T = xT[iT] totinv = 1.0 / abs(ymuop[iT]) for i in range(len(mus)): totinv += 1.0 / abs(ymus[i][iT]) ywmugb.append(1.0 - yKgb[iT]) ywmuop.append(yKgb[iT] / ymuop[iT] / totinv) for i in range(len(mus)): ywmus[i].append(yKgb[iT] / ymus[i][iT] / totinv) print("{:6.3g} {:10.4g} {:10.4g} {:10.4g} {:10.4g} {:10.4g} {:10.4g} {:10.4g}" .format(T, ywmugb[iT], ywmuop[iT], ywmus[0][iT], ywmus[1][iT], ywmus[2][iT], ywmus[3][iT], ywmus[4][iT])) print("") print(f"Save fitting data to [{cparams.outxlsfile}]") fit.to_excel(cparams.outxlsfile, ["T(K)", "mu,obs(cm2/Vs)", "mu,ini(cm2/Vs)", "mu,fin(cm2/Vs)"], [xT, ymu, app.ymuini, ymufin]) print(f"Save component data to [{cparams.outcomponentxlsfile}]") fit.to_excel(cparams.outcomponentxlsfile, ["T(K)", "mu,obs(cm2/Vs)", "mu,tot", "Kgb", "mu,in-grain", "mu,op", f"mu(T^{mu.ppolyi[0]})", f"mu(T^{mu.ppolyi[1]})", f"mu(T^{mu.ppolyi[2]})", f"mu(T^{mu.ppolyi[3]})", f"mu(T^{mu.ppolyi[4]})", "w,GB", "w,op", f"w(T^{mu.ppolyi[0]})", f"w(T^{mu.ppolyi[1]})", f"w(T^{mu.ppolyi[2]})", f"w(T^{mu.ppolyi[3]})", f"w(T^{mu.ppolyi[4]})"], [xT, ymu, ymufin, yKgb, ymu_ingrain, ymuop, ymus[0], ymus[1], ymus[2], ymus[3], ymus[4], ywmugb, ywmuop, ywmus[0], ywmus[1], ywmus[2], ywmus[3], ywmus[4], ]) #============================= # グラフの表示 #============================= print("") print("Plot optimized") ymin = min(ymu) yinmax = max(ymu) ymax = yinmax for m in [ymu_ingrain, ymuop, *ymus]: mumax = max(m) mumin = min(m) if mumax < 1e9 and ymax < mumax: ymax = mumax if mumin > -abs(yinmax) and ymin > mumin: ymin = mumin if 5.0 * yinmax < ymax: ymax = 5.0 * yinmax if ymin < 0.0: ymin *= 1.1 elif ymin == 0.0: ymin = -yinmax * 0.1 ax1.tick_params(labelsize = cparams.fontsize) ax3.tick_params(labelsize = cparams.fontsize) fit.finalize_plot(ymufin, iter = res.nit, fmin = ffin) plot_muT_decomposed(app, cparams, mu, ax2, xT, ymu, ymufin, ymu_ingrain, ymuop, ymus, colors = cparams.colors, markersize = 3.0, plot_event = fit.plot_event) ax2.set_ylim([ymin * 0.8, ymax * 1.1]) plot_muT_weight(app, cparams, mu, ax3, xT, ywmugb, ywmuop, ywmus, colors = cparams.colors, markersize = 3.0, plot_event = fit.plot_event) fit.layout() # plt.tight_layout() plt.pause(0.01) app.terminate(usage = usage, pause = True) def init(app, cparams, mu, fit): ekB = e / kB print("") print("Initialize mobility parameters for mu-T Hall data by linear least-squares fitting") print("infile : ", cparams.infile) print("T_label : ", cparams.T_label) print("mu_label : ", cparams.mu_label) if '***' in cparams.T_label: app.terminate("\nError: Choose T_label", pause = True) if '***' in cparams.mu_label: app.terminate("\nError: Choose mu_label", pause = True) # read data print("") print(f"Read data from [{cparams.infile}]") fit.read_data(cparams.infile, xlabel = cparams.T_label, ylabel = cparams.mu_label, xmin = cparams.Tmin, xmax = cparams.Tmax, usage = usage) xT = fit.x ymu = fit.y print("") print("Paramters") mu.print_parameters() # calculate initial values ys = [] #fit.cal_ylist(fit.pk) ymuT05 = [] print("Calculate mu-T from the initial parameters") print(f"{'T':8}\t{'mu (cm2/Vs)':14}\t{'mu(cal)':14}") for iT in range(len(xT)): T = xT[iT] # muc = cal_mu(mu, T, *fit.pk) muc = mu.cal_mu(T) ys.append(muc) ymuT05.append(ymu[iT] * sqrt(T) * exp(mu.Eb * ekB / T)) print(f"{xT[iT]:8.4f}\t{ymu[iT]:14.6g}\t{ys[iT]:14.6g}") print("") print("Activation energy vs T") xTinv, yEa = Ea_Arrhenius(xT, ymu, eps = 1.0e-300, kTinv = 1000.0) xTinv, yEaT05 = Ea_Arrhenius(xT, ymuT05, eps = 1.0e-300, kTinv = 1000.0) print("f{:8}\t{:8}\t{:14}\t{:14}".format('T', '1000/T', 'Ea,mu (eV)', 'Ea,mu*T^0.5 (eV)')) for iT in range(len(xT)): T = xT[iT] print("f{:8.4f}\t{:8.4f}\t{:14.6g}\t{:14.6g}".format(xT[iT], xTinv[iT], yEa[iT], yEaT05[iT])) print("") print("Linear least-squares fitting") if mu.use_simple == 1: ai, aipoly, aiall, Si, Sij = mu.linearfit3(xT, ymu) else: ai, aipoly, aiall, Si, Sij = mu.linearfit(xT, ymu) mu.apolyi = aipoly mu.ai = aiall mu.set_parameters() print("") print("Fitted parameters:") mu.print_parameters() #varname, varunit, mu.parameter_list()) #, optid) print("") print("Calculate final result") ymufin = [] # print(f"{'T(K)':8}\t{'mu,obs(cm2/Vs)':12}\t{'mu,ini(cm2/Vs)':12}\t{'mu,fin(cm2/Vs)':12}") for iT in range(len(xT)): T = xT[iT] # muc = cal_mu(mu, T, *fit.pk) muc = mu.cal_mu(T) ymufin.append(muc) # print(f"{T:8.4g}\t{ymu[iT]:12.6g}\t{ys[iT]:12.6g}\t{ymufin[iT]:12.6g}") fit.print_scores(heading = "\nScores between y(input) and y(fit)", y1 = fit.y, y2 = ymufin) print("") print("Save parameters to [{}]".format(cparams.parameterfile)) cparams.save_parameters(cparams.parameterfile, "Preferences") mu.save_parameters(cparams.parameterfile) #, optid) print("Save fitting data to [{}]".format(cparams.outxlsfile)) fit.to_excel(cparams.outxlsfile, ["T(K)", "mu,obs(cm2/Vs)", "mu,ini(cm2/Vs)", "mu,fin(cm2/Vs)"], [xT, ymu, ys, ymufin]) #============================= # グラフの表示 #============================= print("") app.plt = plt app.fig = app.plt.figure(figsize = cparams.figsize) ax1 = app.fig.add_subplot(1, 2, 1) ax1.plot(xT, ymu, label = '$\mu(obs)$', linestyle = '', marker = 'o') # ax1.plot(xT, app.ymuini, label = '$\mu(ini)$', color = 'blue', linewidth = cparams.linewidth, linestyle = 'dashed') ax1.plot(xT, ymufin, label = '$\mu(fin)$', color = 'red', linewidth = cparams.linewidth, linestyle = '-') ax1.set_xlabel("T (K)", fontsize = cparams.fontsize) ax1.set_ylabel("$\mu$ (cm$^2$/Vs)", fontsize = cparams.fontsize) ax1.legend(fontsize = cparams.legend_fontsize) ax1.tick_params(labelsize = cparams.fontsize) app.plt.tight_layout() app.plt.pause(0.1) app.terminate("", usage = usage, pause = True) def sim(app, cparams, mu, fit): print("") print("Initialize mobility parameters for mu-T Hall data by linear least-squares fitting") print("infile : ", cparams.infile) print("T_label : ", cparams.T_label) # print("n_label : ", cparams.n_label) print("mu_label : ", cparams.mu_label) # print("sigma_label: ", cparams.sigma_label) if '***' in cparams.T_label: app.terminate("\nError: Choose T_label", pause = True) if '***' in cparams.mu_label: app.terminate("\nError: Choose mu_label", pause = True) # read data print("") print(f"Read data from [{cparams.infile}]") fit.read_data(cparams.infile, xlabel = cparams.T_label, ylabel = cparams.mu_label, xmin = cparams.Tmin, xmax = cparams.Tmax, usage = usage) xT = fit.x ymu = fit.y # calculate initial values print("") print("Calculate mu-T from the initial parameters") print("{:8}\t{:12}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}" .format("T(K)", "mu,obs(cm2/Vs)", "mu,tot", "Kgb", "mu,in-grain", "mu,op", f"mu(T^{mu.ppolyi[0]})", f"mu(T^{mu.ppolyi[1]})", f"mu(T^{mu.ppolyi[2]})", f"mu(T^{mu.ppolyi[3]})", f"mu(T^{mu.ppolyi[4]})")) app.ymuini = [] yKgb = [] ymu_ingrain = [] ymuop = [] ymus = make_matrix1(len(mu.ppolyi), type = 'list') # nskip = int(len(xT) / 10) nskip = 1 for iT in range(0, len(xT), nskip): T = xT[iT] mu_tot, Kgb, mu_ingrain, muop, mus = mu.cal_mu_components(T) # mu_tot = mu.cal_mu(T) # muc_tot = cal_mu(mu, T, *fit.pk) app.ymuini.append(mu_tot) yKgb.append(Kgb) ymu_ingrain.append(mu_ingrain) ymuop.append(muop) for i in range(len(mus)): ymus[i].append(mus[i]) print("{:8.4g}\t{:12.6g}\t{:10.6g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}" .format(T, ymu[iT], mu_tot, Kgb, mu_ingrain, muop, mus[0], mus[1], mus[2], mus[3], mus[4])) xTinv = [] ymuT05 = [] ymuiniT05 = [] print("") print(f"{'T':8}\t{'1000/T':8}\t{'mu (cm2/Vs)':14}\t{'mu*T^0.5 (cm2K0.5/Vs)':14}\t{'mu(cal)':14}\t{'mu(cal)*T^0.5 (cm2K0.5/Vs)':14}") for i in range(len(xT)): xTinv.append(1000.0 / xT[i]) ymuT05.append(ymu[i] * sqrt(xT[i])) ymuiniT05.append(app.ymuini[i] * sqrt(xT[i])) print(f"{xT[i]:8.4f}\t{xTinv[i]:8.4f}\t{ymu[i]:14.6g}\t{ymuT05[i]:14.6g}\t{app.ymuini[i]:14.6g}\t{ymuiniT05[i]:14.6g}") print("") print("Activation energy vs T") Tinv, yEa = Ea_Arrhenius(xT, ymu, eps = 1.0e-300, kTinv = 1000.0) Tinv, yEaT05 = Ea_Arrhenius(xT, ymuT05, eps = 1.0e-300, kTinv = 1000.0) print("{:8}\t{:8}\t{:14}\t{:14}".format('T', '1000/T', 'Ea,mu (eV)', 'Ea,mu*T^0.5 (eV)')) for iT in range(len(xT)): T = xT[iT] print("{:8.4f}\t{:8.4f}\t{:14.6g}\t{:14.6g}".format(xT[iT], xTinv[iT], yEa[iT], yEaT05[iT])) print("") print("Save parameters to [{}]".format(cparams.parameterfile)) cparams.save_parameters(cparams.parameterfile, "Preferences") mu.save_parameters(cparams.parameterfile) print("Save simulation data to [{}]".format(cparams.outxlsfile)) fit.to_excel(cparams.outxlsfile, ["T(K)", "mu,obs(cm2/Vs)", "mu,ini(cm2/Vs)", "Ea,mu(eV)", "Ea,mu*T^0.5(eV)"], [xT, ymu, app.ymuini, yEa, yEaT05]) print("Save component data to [{}]".format(cparams.outcomponentxlsfile)) fit.to_excel(cparams.outcomponentxlsfile, ["T(K)", "mu,obs(cm2/Vs)", "mu,tot", "Kgb", "mu,in-grain", "mu,op", f"mu(T^{mu.ppolyi[0]})", f"mu(T^{mu.ppolyi[1]})", f"mu(T^{mu.ppolyi[2]})", f"mu(T^{mu.ppolyi[3]})", f"mu(T^{mu.ppolyi[4]})"], [xT, ymu, app.ymuini, yKgb, ymu_ingrain, ymuop, ymus[0], ymus[1], ymus[2], ymus[3], ymus[4]]) #============================= # グラフの表示 #============================= print("") fig = plt.figure(figsize = cparams.figsize) app.plt = plt app.fig = fig fit.initialize_plot(fig = fig) ax1 = fig.add_subplot(2, 3, 1) ax2 = fig.add_subplot(2, 3, 4) ax3 = fig.add_subplot(2, 3, 2) ax4 = fig.add_subplot(2, 3, 5) ax5 = fig.add_subplot(2, 3, 3) ax6 = fig.add_subplot(2, 3, 6) ymin = min([min(app.ymuini), min(ymu)]) ymax = max([max(app.ymuini), max(ymu)]) plot_muT_decomposed(app, cparams, mu, ax1, xT, ymu, app.ymuini, ymu_ingrain, ymuop, ymus, colors = cparams.colors, plot_event = fit.plot_event) ax1.set_ylim([ymin * 0.5, ymax * 2.0]) ax1.legend(fontsize = 8) #legend_fontsize) ax1.tick_params(labelsize = cparams.fontsize) plot_muT_decomposed(app, cparams, mu, ax2, xT, ymu, app.ymuini, ymu_ingrain, ymuop, ymus, colors = cparams.colors, plot_event = fit.plot_event) ax2.set_ylim([ymin * 0.1, ymax * 10.0]) ax2.set_xscale('log') ax2.set_yscale('log') ax2.legend(fontsize = 8) #legend_fontsize) ax2.tick_params(labelsize = cparams.fontsize) plot_muT_decomposed(app, cparams, mu, ax3, xTinv, ymu, app.ymuini, ymu_ingrain, ymuop, ymus, xlabel = "1000/T (K$^{-1}$)", colors = cparams.colors, plot_event = fit.plot_event) ax3.set_ylim([ymin * 0.5, ymax * 2.0]) ax3.set_yscale('log') ax3.legend(fontsize = 8) #legend_fontsize) ax3.tick_params(labelsize = cparams.fontsize) ax3.set_ylim([ymin * 0.5, ymax * 2.0]) ax4.plot(xTinv, ymuT05, label = '$\mu(obs)\cdot$$T^{1/2}$', linestyle = '', marker = 'o', markersize = 1.0) ax4.plot(xTinv, ymuiniT05, label = '$\mu(ini)\cdot$$T^{1/2}$', color = 'blue', linewidth = 1.0, linestyle = 'dashed') ax4.set_xlabel("1000/T (K$^{-1}$)", fontsize = cparams.fontsize) ax4.set_ylabel("$\mu\cdot$$T^{1/2}$ (cm$^2$/Vs)", fontsize = cparams.fontsize) ax4.set_yscale('log') ax4.legend(fontsize = cparams.legend_fontsize) ax4.tick_params(labelsize = cparams.fontsize) ax5.plot(xT, yEa, label = '$E_a$ (eV)', linestyle = '-', linewidth = 1.0, marker = 'o', markersize = 1.0) ax5.set_xlabel("$T$ (K)", fontsize = cparams.fontsize) ax5.set_ylabel("$E_a$ (eV)", fontsize = cparams.fontsize) ax5.legend(fontsize = cparams.legend_fontsize) ax5.tick_params(labelsize = cparams.fontsize) ax6.plot(xT, yEaT05, label = '$E_a$ (eV)', linestyle = '-', linewidth = 1.0, marker = 'o', markersize = 1.0) ax6.set_xlabel("$T$ (K)", fontsize = cparams.fontsize) ax6.set_ylabel("$E_a$($\mu$$T^{1/2}$) (eV)", fontsize = cparams.fontsize) ax6.legend(fontsize = cparams.legend_fontsize) ax6.tick_params(labelsize = cparams.fontsize) plt.tight_layout() plt.pause(0.1) app.terminate("", usage = usage, pause = True) def main(): app, cparams, mu = initialize() updatevars(app, cparams, mu) cparams.logfile = app.replace_path(cparams.infile) print(f"Open logfile [{cparams.logfile}]") app.redirect(targets = ["stdout", cparams.logfile], mode = 'w') cparams.parameterfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}.in"]) cparams.parameterbkfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}.in.back"]) 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, section = "Preferences") # check args again to use given parameters updatevars(app, cparams, mu) fit = tkFit(tol = cparams.tol, nmaxiter = cparams.maxiter, print_interval = cparams.print_interval, plot_interval = cparams.plot_interval) cparams.outxlsfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-muT-{mode}.xlsx"], ext_dict = {"mode": cparams.mode}) cparams.outcomponentxlsfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-mu-components.xlsx"]) cparams.outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-out.xlsx"]) mu.read_parameters(cparams.parameterfile) #, optid) mu.set_misc_parameters(cparams.get_dict()) mu.ai0 = mu.parameter_list() # mu.ai0 = mu.parameter_tuple() # mu.optid = tuple(mu.optid) print("") print("Initial parameters") if mu.Eop <= 0.0: print("#####################################################################################################") print(" Warning: Eop={:8.3g} <= 0.0 is not allowed.".format(mu.Eop)) print(" Eop is set to 1.0e-6.") print(" Aop is set to 1.0e-10 so that optical phonon scattering will not affect the total mobility.") print("#####################################################################################################") mu.Eop = 1.0e-6 mu.Aop = 1.0e-10 mu.print_parameters() #varname, varunit, mu.parameter_list()) #, optid) if cparams.mode == 'init': init(app, cparams, mu, fit) elif cparams.mode == 'fit': fitting(app, cparams, mu, fit) elif cparams.mode == 'sim': sim(app, cparams, mu, fit) else: app.terminate("Error in main: Invalide mode [{}]".format(mode), usage = usage, pause = True) if __name__ == "__main__": main()