"""
非縮退半導体のEFの温度依存性を二分法で計算。
移動度-温度ホールデータ（キャリア移動度と温度の関係）を解析し、複数の散乱機構による移動度を評価するスクリプトです。
粒界散乱、光学フォノン散乱、および複数のべき乗則散乱機構を考慮した移動度モデルを用いて、測定された移動度-温度データにフィッティングを行います。
初期パラメータの推定（`init`モード）、詳細な非線形フィッティング（`fit`モード）、および既存のパラメータに基づくシミュレーション（`sim`モード）をサポートします。
:doc:`mu-T-Fit2_usage`
"""
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
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():
"""アプリケーションのグローバル変数と初期設定を初期化する。
tkApplication、tkMobility、tkParamsなどのオブジェクトを生成し、
デフォルトのパラメータ値、ファイルパス、グラフ設定などを設定します。
:returns: app, cparams, mu: 初期化されたアプリケーションオブジェクト、設定パラメータオブジェクト、移動度モデルオブジェクト。
:rtype: tuple[tkApplication, tkParams, tkMobility]
"""
#================================
# 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_str`に定義された使用法文字列を標準出力に表示します。
:param app: tkApplicationオブジェクト。
:type app: tklib.tkapplication.tkApplication
"""
# 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):
"""コマンドライン引数から設定パラメータを更新する。
sys.argvを解析し、モード、入力ファイル、データラベル、温度範囲、最適化パラメータなどを
設定オブジェクト`cparams`に適用します。入力ファイルが指定されていない場合はエラーで終了します。
:param app: tkApplicationオブジェクト。
:type app: tklib.tkapplication.tkApplication
:param cparams: 設定パラメータオブジェクト。
:type cparams: tklib.tkparams.tkParams
:param mu: tkMobilityオブジェクト。
:type mu: tklib.tktransport.tkmobility.tkMobility
"""
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):
"""移動度と各散乱機構による分解された移動度をプロットする。
観測された移動度、計算された総移動度、粒内移動度、光学フォノン移動度、
および複数のべき乗則散乱移動度を温度に対してプロットします。
`tkPlotEvent`を通じてデータを管理します。
:param app: tkApplicationオブジェクト。
:type app: tklib.tkapplication.tkApplication
:param cparams: 設定パラメータオブジェクト。
:type cparams: tklib.tkparams.tkParams
:param mu: tkMobilityオブジェクト。
:type mu: tklib.tktransport.tkmobility.tkMobility
:param ax: matplotlibのAxesオブジェクト。
:type ax: matplotlib.axes.Axes
:param xT: 温度データのリストまたはNumpy配列。
:type xT: list[float] or numpy.ndarray
:param ymu: 観測された移動度データのリストまたはNumpy配列。
:type ymu: list[float] or numpy.ndarray
:param ymucal: (オプション) 計算された総移動度データのリストまたはNumpy配列。
:type ymucal: list[float] or numpy.ndarray or None
:param ymu_ingrain: (オプション) 粒内移動度データのリストまたはNumpy配列。
:type ymu_ingrain: list[float] or numpy.ndarray or None
:param ymuop: (オプション) 光学フォノン散乱による移動度データのリストまたはNumpy配列。
:type ymuop: list[float] or numpy.ndarray or None
:param ymus: (オプション) 複数のべき乗則散乱による移動度データのリストまたはNumpy配列。
:type ymus: list[list[float]] or None
:param xlabel: (オプション) X軸のラベル文字列。
:type xlabel: str
:param ylabel: (オプション) Y軸のラベル文字列。
:type ylabel: str
:param colors: (オプション) プロットに使用する色のリスト。
:type colors: list[str] or None
:param markersize: (オプション) マーカーのサイズ。
:type markersize: float
:param plot_event: (オプション) tkPlotEventオブジェクト。
:type plot_event: tklib.tkgraphic.tkplotevent.tkPlotEvent or 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):
"""移動度を構成する各散乱機構の相対的な重みをプロットする。
粒界散乱、光学フォノン散乱、および複数のべき乗則散乱の各機構が全体移動度に与える重み（逆数和における割合）を
温度に対してプロットします。`tkPlotEvent`を通じてデータを管理します。
:param app: tkApplicationオブジェクト。
:type app: tklib.tkapplication.tkApplication
:param cparams: 設定パラメータオブジェクト。
:type cparams: tklib.tkparams.tkParams
:param mu: tkMobilityオブジェクト。
:type mu: tklib.tktransport.tkmobility.tkMobility
:param ax: matplotlibのAxesオブジェクト。
:type ax: matplotlib.axes.Axes
:param xT: 温度データのリストまたはNumpy配列。
:type xT: list[float] or numpy.ndarray
:param ywmugb: 粒界散乱の重みデータのリストまたはNumpy配列。
:type ywmugb: list[float] or numpy.ndarray
:param ywmuop: 光学フォノン散乱の重みデータのリストまたはNumpy配列。
:type ywmuop: list[float] or numpy.ndarray
:param ywmus: 複数のべき乗則散乱の重みデータのリストまたはNumpy配列。
:type ywmus: list[list[float]]
:param xlabel: (オプション) X軸のラベル文字列。
:type xlabel: str
:param ylabel: (オプション) Y軸のラベル文字列。
:type ylabel: str
:param colors: (オプション) プロットに使用する色のリスト。
:type colors: list[str] or None
:param markersize: (オプション) マーカーのサイズ。
:type markersize: float
:param plot_event: (オプション) tkPlotEventオブジェクト。
:type plot_event: tklib.tkgraphic.tkplotevent.tkPlotEvent or 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):
"""与えられたパラメータセットで移動度を計算する。
`tkMobility`オブジェクトの一時的なパラメータ更新を行い、指定された温度における
総移動度を計算します。計算後、パラメータは元の値に戻されます。
:param mu: tkMobilityオブジェクト。
:type mu: tklib.tktransport.tkmobility.tkMobility
:param T: 移動度を計算する温度 (K)。
:type T: float
:param pk: 可変長引数として渡される移動度モデルパラメータのタプル。
:type pk: tuple
:returns: 指定された温度Tにおける計算された移動度 (cm^2/Vs)。
:rtype: float
"""
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):
"""観測された移動度-温度データに非線形最小二乗フィッティングを行う。
指定された入力ファイルから移動度-温度データを読み込み、
`scipy.optimize.minimize`を用いて、粒界散乱、光学フォノン散乱、
複数のべき乗則散乱を考慮した移動度モデルのパラメータを最適化します。
フィッティング結果、各散乱機構の貢献度、およびプロットを出力します。
:param app: tkApplicationオブジェクト。
:type app: tklib.tkapplication.tkApplication
:param cparams: 設定パラメータオブジェクト。
:type cparams: tklib.tkparams.tkParams
:param mu: tkMobilityオブジェクト。
:type mu: tklib.tktransport.tkmobility.tkMobility
:param fit: tkFitオブジェクト。
:type fit: tklib.tksci.tkFit.tkFit
"""
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):
"""移動度モデルパラメータを線形最小二乗法で初期化する。
観測された移動度-温度データに対し、特定の形式に変換して線形最小二乗フィッティングを行い、
移動度モデルのべき乗則散乱項の係数を初期推定します。
活性化エネルギーの計算とプロットも行います。
:param app: tkApplicationオブジェクト。
:type app: tklib.tkapplication.tkApplication
:param cparams: 設定パラメータオブジェクト。
:type cparams: tklib.tkparams.tkParams
:param mu: tkMobilityオブジェクト。
:type mu: tklib.tktransport.tkmobility.tkMobility
:param fit: tkFitオブジェクト。
:type fit: tklib.tksci.tkFit.tkFit
"""
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):
"""既存のパラメータを用いて移動度-温度特性をシミュレーションする。
指定された入力ファイルから温度と観測移動度データを読み込み、
現在のパラメータ設定に基づいて移動度、各散乱機構の貢献度、
および活性化エネルギーを計算します。
結果をファイルに出力し、詳細なプロットを表示します。
:param app: tkApplicationオブジェクト。
:type app: tklib.tkapplication.tkApplication
:param cparams: 設定パラメータオブジェクト。
:type cparams: tklib.tkparams.tkParams
:param mu: tkMobilityオブジェクト。
:type mu: tklib.tktransport.tkmobility.tkMobility
:param fit: tkFitオブジェクト。
:type fit: tklib.tksci.tkFit.tkFit
"""
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():
"""プログラムのメインエントリポイント。
アプリケーションの初期化、コマンドライン引数の解析と設定の更新、
ログファイルの設定、パラメータファイルの読み込み、および
選択されたモード（`init`, `fit`, `sim`）に応じた処理の実行を制御します。
"""
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(cparams.mode), usage = usage, pause = True)
if __name__ == "__main__":
main()