"""
半導体の散乱モデルに基づいたホール移動度(μ-T)データの最小二乗法フィッティングツール。
概要:
このスクリプトは、半導体のホール移動度（μ-T）データの非線形最小二乗法フィッティングを実行します。
複数の散乱機構（光学的フォノン散乱、中性不純物散乱、格子振動散乱）と粒界散乱モデル（KGBモデル）
を考慮したモデルを用いて、移動度の温度依存性を解析します。
詳細説明:
与えられたホール移動度（μ）と温度（T）の関係データに対し、以下の散乱機構を組み合わせたモデルを適用します。
- 光学的フォノン散乱 (Optical Phonon Scattering): μ_op^(-1) = A_op / (exp(E_op / kBT) - 1)
- 中性不純物散乱 (Neutral Impurity Scattering): μ_T0^(-1) = A_T0
- 格子振動散乱 (Acoustic Phonon Scattering): μ_T3/2^(-1) = A_T3/2 * T^(-1.5)
- 粒界散乱 (Grain Boundary Scattering): μ = K_GB * μ_in-grain
ここで、K_GB = exp(-V_B / kBT + (s_phi / kBT)^2 / 2)
μ_in-grain^(-1) = μ_op^(-1) + μ_T0^(-1) + μ_T3/2^(-1)
スクリプトは、フィッティングモード ('fit')、初期パラメータ設定モード ('init')、シミュレーションモード ('sim')
の3つの主要なモードで動作します。
- 'fit'モードでは、非線形最小二乗法を用いて最適化を行い、モデルパラメータを決定します。
- 'init'モードでは、線形最小二乗法に近い形で初期パラメータを探索します。
- 'sim'モードでは、既存のパラメータセットに基づき移動度を計算し、観測データと比較します。
結果として、最適化されたパラメータ、各散乱機構による移動度成分、およびそれらの全体移動度に対する寄与度を計算し、
Excelファイルへの保存とグラフ表示を行います。
関連リンク:
:doc:`mu-T-Fit_tkFit_usage`
"""
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
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
usage_str = '''
""
f"(i) usage: python {sys.argv[0]} fit mode infile xlabel ylabel Tmin Tmax Tcalmin Tcalmax"
'''[1:-1]
[ドキュメント]
def initialize():
"""
アプリケーションとパラメータを初期化する。
概要:
`tkApplication`と`tkParams`オブジェクトを作成し、スクリプト全体の制御パラメータを設定する。
詳細説明:
フィッティングモード、入力ファイルパス、データ列指定（xlabel, ylabel）、
フィッティングおよび計算対象の温度範囲、最適化アルゴリズム、イテレーション数、許容誤差、
グラフ設定（色、サイズ、フォントサイズなど）、そして移動度モデルの初期パラメータ
（Eop, Aop, AT0, AT32, VB, s_phi）がここで定義されます。
戻り値:
:returns app: tkApplicationのインスタンス。アプリケーションの全体的な制御とユーティリティを提供。
:returns cparams: tkParamsのインスタンス。アプリケーション全体の制御パラメータを保持。
"""
#================================
# Global variables
#================================
app = tkApplication(usage_str = usage_str)
argv, narg = app.get_argv()
app.cparams = tkParams()
cparams = app.cparams
cparams.debug = 0
cparams.print_level = 0
# mode: 'init', 'fit', 'sim'
# cparams.mode = 'init'
cparams.mode = 'fit'
cparams.fix_N = False
# data file to fit (.xlsx)
cparams.infile = 'IO-Hall-02A01.xlsx'
cparams.parameterfile = None
cparams.parameterbkfile = None
cparams.xlabel = 2
cparams.ylabel = 6
cparams.T_label = 6
cparams.mu_label = 7
# 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)
#scipy.optimize.minimizeで使うアルゴリズム
#nelder-mead Downhill simplex
#cg conjugate gradient (Polak-Ribiere method)
#bfgs BFGS法
#newton-cg Newton-CG
cparams.method = "nelder-mead"
cparams.nmaxiter = 3000
cparams.tol = 1.0e-3
cparams.h_diff = 1.0e-3
cparams.print_interval = 5
cparams.plot_interval = 1
#=============================
# Graph configuration
#=============================
cparams.fig = None
cparams.linewidth = 2.0
cparams.colors = ['black', 'red', 'blue', 'green', 'orange', 'darkgreen', 'darkorange', 'navy',
'slategray', 'hotpink', 'olive', 'chocolate', 'magenta',
'green', 'yellow', 'cyan']
cparams.figsize = [12, 6]
cparams.fontsize = 14
cparams.legend_fontsize = 12
cparams.graphupdateinterval = 10
app.mu_params = tkParams()
app.mu_params.Eop = 0.046
app.mu_params.Aop = 0.0297
app.mu_params.AT0 = 0.0138
app.mu_params.AT32 = 0.7186
app.mu_params.VB = 0.0
app.mu_params.s_phi = 0.0
return app, cparams
#=============================
# Treat argments
#=============================
[ドキュメント]
def usage(app):
"""
スクリプトの利用方法メッセージを表示する。
概要:
`app.usage_str`に定義されたスクリプトの使用方法を標準出力に表示する。
詳細説明:
コマンドライン引数の指定方法や、各引数の意味についてユーザーに情報を提供します。
通常、エラー発生時やヘルプが要求された際に呼び出されます。
引数:
:param app: tkApplicationのインスタンス。使用方法の文字列を保持。
"""
cparams = app.cparams
for s in app.usage_str.split('\n'):
cmd = 'print({})'.format(s.rstrip())
eval(cmd)
[ドキュメント]
def update_vars(app, cparams):
"""
コマンドライン引数に基づいて制御パラメータを更新する。
概要:
スクリプト実行時に与えられたコマンドライン引数を解析し、`cparams`の対応する設定を上書きする。
詳細説明:
`tklib.tkutils`モジュールから提供される`getarg`, `getintarg`, `getfloatarg`関数を使用して、
モード、最適化手法、入力ファイル、プロットの軸、温度範囲、最大反復回数、許容誤差、差分計算のステップサイズなど、
様々な制御パラメータをコマンドライン引数の値で更新します。
引数:
:param app: tkApplicationのインスタンス。引数取得に使用。
:param cparams: tkParamsのインスタンス。更新対象の制御パラメータを保持。
"""
cparams.mode = getarg(1, cparams.mode)
cparams.method = getarg(2, cparams.method)
cparams.infile = getarg(3, cparams.infile)
cparams.xlabel = getarg(4, cparams.xlabel)
cparams.ylabel = getarg(5, cparams.ylabel)
cparams.Tmin = getfloatarg(6, cparams.Tmin)
cparams.Tmax = getfloatarg(7, cparams.Tmax)
cparams.Tcalmin = getfloatarg(8, cparams.Tcalmin)
cparams.Tcalmax = getfloatarg(9, cparams.Tcalmax)
cparams.nmaxiter = getintarg (10, cparams.nmaxiter)
cparams.tol = getfloatarg(11, cparams.tol)
cparams.h_diff = getfloatarg(11, cparams.h_diff)
[ドキュメント]
def muinv_op(T, Eop, Aop):
"""
光学的フォノン散乱による移動度の逆数を計算する。
概要:
温度Tにおける光学的フォノン散乱による移動度逆数 (1/μ_op) を計算する。
詳細説明:
光学的フォノン散乱は、高温領域で主要な散乱機構の一つです。
`ekBT = e / kB / T` は熱エネルギーの逆数に比例する量で、Boltzmann定数`kB`と素電荷`e`が使用されます。
引数:
:param T: float 温度 (K)。
:param Eop: float 光学的フォノンのエネルギー (eV)。
:param Aop: float 光学的フォノン散乱係数。
戻り値:
:returns: float 光学的フォノン散乱による移動度の逆数。
"""
ekBT = e / kB / T
return Aop / (exp(Eop * ekBT) - 1.0)
[ドキュメント]
def muinv_T0(T, AT0):
"""
中性不純物散乱（またはT^0依存性散乱）による移動度の逆数を計算する。
概要:
温度Tにおける中性不純物散乱による移動度逆数 (1/μ_T0) を計算する。
詳細説明:
この散乱機構は温度にほとんど依存しない成分として扱われ、主に低温領域で影響を与えることがあります。
引数:
:param T: float 温度 (K)。
:param AT0: float 中性不純物散乱係数。
戻り値:
:returns: float 中性不純物散乱による移動度の逆数。
"""
return AT0
[ドキュメント]
def muinv_T32(T, AT32):
"""
格子振動散乱（またはT^-1.5依存性散乱）による移動度の逆数を計算する。
概要:
温度Tにおける格子振動散乱（音響フォノン散乱など）による移動度逆数 (1/μ_T3/2) を計算する。
詳細説明:
音響フォノン散乱などの格子振動による散乱は、通常、移動度がT^(-1.5)に比例する形で記述されます。
これは中〜高温領域で重要な散乱機構です。
引数:
:param T: float 温度 (K)。
:param AT32: float 格子振動散乱係数。
戻り値:
:returns: float 格子振動散乱による移動度の逆数。
"""
return AT32 * pow(T, -1.5)
[ドキュメント]
def KGB(T, VB, s_phi):
"""
粒界散乱モデル (KGBモデル) の伝導率因子を計算する。
概要:
温度Tにおける粒界散乱モデル（KGBモデル）に基づく伝導率因子を計算する。
詳細説明:
KGB (Kainthla-Gandhi-Bhalla) モデルは、粒界ポテンシャル障壁による電子の散乱を記述します。
`ekBT = e / kB / T` は熱エネルギーの逆数に比例する量で、Boltzmann定数`kB`と素電荷`e`が使用されます。
`VB`は粒界ポテンシャルの平均値、`s_phi`は粒界ポテンシャルのばらつきを表します。
引数:
:param T: float 温度 (K)。
:param VB: float 粒界ポテンシャルエネルギー (eV)。
:param s_phi: float 粒界ポテンシャルのばらつき (eV)。
戻り値:
:returns: float KGBモデルによる伝導率因子。
"""
ekBT = e / kB / T
return exp(-VB * ekBT + (s_phi * ekBT)**2 / 2.0)
[ドキュメント]
def cal_mu_components(T, _Eop, _Aop, _AT0, _AT32, _VB, _s_phi):
"""
移動度の各散乱機構からの成分と合計移動度を計算する。
概要:
光学的フォノン、中性不純物、格子振動散乱による移動度成分と、粒界散乱の影響を考慮した合計移動度を計算する。
詳細説明:
各散乱機構（光学的フォノン、中性不純物、格子振動）による移動度逆数を合計して粒内移動度を求め、
それにKGBモデルの伝導率因子を乗じることで合計移動度を算出します。
計算中にゼロ除算が発生する可能性がある場合、警告を出力し、安定した結果を得るために
対応する移動度成分を大きな値（1e10）に補正します。
引数:
:param T: float 温度 (K)。
:param _Eop: float 光学的フォノンのエネルギー (eV)。
:param _Aop: float 光学的フォノン散乱係数。
:param _AT0: float 中性不純物散乱係数。
:param _AT32: float 格子振動散乱係数。
:param _VB: float 粒界ポテンシャルエネルギー (eV)。
:param _s_phi: float 粒界ポテンシャルのばらつき (eV)。
戻り値:
:returns mu_tot: float 合計移動度 (cm^2/Vs)。
:returns mu_KGB: float KGBモデルによる伝導率因子。
:returns mu_ingrain: float 粒内移動度 (cm^2/Vs)。
:returns mu_op: float 光学的フォノン散乱による移動度 (cm^2/Vs)。
:returns mu_T0: float 中性不純物散乱による移動度 (cm^2/Vs)。
:returns mu_T32: float 格子振動散乱による移動度 (cm^2/Vs)。
"""
mui_op = muinv_op(T, _Eop, _Aop)
mui_T0 = muinv_T0(T, _AT0)
mui_T32 = muinv_T32(T, _AT32)
mu_ingrain = 1.0 / (mui_op + mui_T0 + mui_T32)
mu_KGB = KGB(T, _VB, _s_phi)
mu_tot = mu_KGB * mu_ingrain
if mu_ingrain == 0.0:
print("\n***Warning: mu_ingrains is 0.0. Correct to 1e10")
mu_ingrain = 1.0e10
if mui_op == 0.0:
print("\n***Warning: mui_op is 0.0. Correct to 1e10")
mui_op = 1.0e10
if mui_T0 == 0.0:
print("\n***Warning: mui_T0 is 0.0. Correct to 1e10")
mui_T0 = 1.0e10
if mui_T32 == 0.0:
print("\n***Warning: mui_T32 is 0.0. Correct to 1e10")
mui_T32 = 1.0e10
return mu_tot, mu_KGB, mu_ingrain, 1.0 / mui_op, 1.0 / mui_T0, 1.0 / mui_T32
[ドキュメント]
def cal_mu_from_params(T, _Eop, _Aop, _AT0, _AT32, _VB, _s_phi):
"""
指定されたパラメータと温度に基づいて合計移動度を計算する。
概要:
`cal_mu_components`関数を呼び出し、その結果から合計移動度 (`mu_tot`) のみを返す。
詳細説明:
この関数は主に`scipy.optimize.minimize`などの最適化ルーチンに渡されるコールバック関数として設計されています。
フィッティングの目的関数内で、与えられたパラメータセットから計算される移動度を提供します。
引数:
:param T: float 温度 (K)。
:param _Eop: float 光学的フォノンのエネルギー (eV)。
:param _Aop: float 光学的フォノン散乱係数。
:param _AT0: float 中性不純物散乱係数。
:param _AT32: float 格子振動散乱係数。
:param _VB: float 粒界ポテンシャルエネルギー (eV)。
:param _s_phi: float 粒界ポテンシャルのばらつき (eV)。
戻り値:
:returns: float 合計移動度 (cm^2/Vs)。
"""
mu_tot, mu_KGB, mu_ingrain, mu_op, mu_T0, mu_T32 = cal_mu_components(T, _Eop, _Aop, _AT0, _AT32, _VB, _s_phi)
return mu_tot
# muinv = muinv_op(T, _Eop, _Aop)
# muinv += muinv_T0(T, _AT0)
# muinv += muinv_T32(T, _AT32)
# mu = KGB(T, _VB, _s_phi) / muinv
#
# return mu
[ドキュメント]
def plot_muT_decomposed(app, cparams, fit, ax, xT, ymu, ymucal = None, ymu_ingrain = None, ymuop = None, ymuT0 = None, ymuT32 = None,
xlabel = "T (K)", ylabel = "$\mu$ (cm$^2$/Vs)", colors = None, markersize = 1.0, plot_event = None):
"""
ホール移動度の観測値と計算値、および分解された移動度成分をプロットする。
概要:
指定されたMatplotlib Axesオブジェクトに、観測されたホール移動度、フィッティングされた移動度、
そして各散乱機構（光学的フォノン、中性不純物、格子振動）による移動度成分をプロットする。
詳細説明:
観測データはマーカーで、計算結果は実線でプロットされます。
粒内移動度、各散乱成分移動度は破線で表示され、各成分の寄与を視覚的に理解するのに役立ちます。
軸ラベル、凡例、マーカーサイズ、色などのプロット設定は`cparams`から取得されます。
引数:
:param app: tkApplicationのインスタンス。
:param cparams: tkParamsのインスタンス。グラフの描画設定を保持。
:param fit: tkFitのインスタンス。データおよびフィッティングに関する情報を保持。
:param ax: matplotlib.axes.Axes プロット対象のMatplotlib Axesオブジェクト。
:param xT: numpy.ndarray 温度データの配列 (K)。
:param ymu: numpy.ndarray 観測された移動度データの配列 (cm$^2$/Vs)。
:param ymucal: numpy.ndarray, optional 計算された（フィッティングされた）移動度データの配列。デフォルトはNone。
:param ymu_ingrain: numpy.ndarray, optional 粒内移動度データの配列。デフォルトはNone。
:param ymuop: numpy.ndarray, optional 光学的フォノン散乱による移動度データの配列。デフォルトはNone。
:param ymuT0: numpy.ndarray, optional 中性不純物散乱による移動度データの配列。デフォルトはNone。
:param ymuT32: numpy.ndarray, optional 格子振動散乱による移動度データの配列。デフォルトはNone。
:param xlabel: str, optional X軸のラベル。デフォルトは"T (K)"。
:param ylabel: str, optional Y軸のラベル。デフォルトは"$\mu$ (cm$^2$/Vs)"。
:param colors: list, optional プロットに使用する色のリスト。デフォルトはNone。
:param markersize: float, optional データポイントのマーカーサイズ。デフォルトは1.0。
:param plot_event: tkPlotEvent, optional プロットイベントハンドラ。プロットデータ登録に使用。デフォルトはNone。
戻り値:
None
"""
ax.clear()
ax.tick_params(labelsize = cparams.fontsize)
data = ax.plot(xT, ymu, label = '$\mu(obs)$', linestyle = '', marker = 'o', markersize = markersize)
plot_event.add_data({"label": "mu_obs", "plot_type": "plot", "axis": ax, "data": data})
if ymucal:
data = ax.plot(xT, ymucal, label = '$\mu(cal)$', color = 'red', 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 = 1.0, 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 = 1.0, linestyle = 'dashed', color = color,
marker = 'o', markersize = 2.0)
plot_event.add_data({"label": "mu_op", "plot_type": "plot", "axis": ax, "data": data})
if ymuT0:
color = colors[2]
data = ax.plot(xT, ymuT0, label = '$\mu_{T0}$', linewidth = 1.0, linestyle = 'dashed', color = color,
marker = 'o', markersize = 2.0)
plot_event.add_data({"label": "mu_op", "plot_type": "plot", "axis": ax, "data": data})
if ymuT32:
color = colors[3]
data = ax.plot(xT, ymuT32, label = '$\mu_{T3/2}$', linewidth = 1.0, linestyle = 'dashed', color = color,
marker = 'o', markersize = 2.0)
plot_event.add_data({"label": "mu_op", "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, fit, ax, xT, ywmugb, ywmuop, ywmuT0, ywmuT32,
xlabel = "T (K)", ylabel = "Linear weight", colors = None, markersize = 1.0, plot_event = None):
"""
各散乱機構の線形重み（寄与度）をプロットする。
概要:
指定されたMatplotlib Axesオブジェクトに、粒界散乱、光学的フォノン散乱、
中性不純物散乱、格子振動散乱の各機構の線形重みをプロットする。
詳細説明:
各散乱機構が全体移動度にどれだけ寄与しているかを線形重みとしてプロットします。
これにより、特定の温度範囲でどの散乱機構が支配的であるかを視覚的に分析できます。
軸ラベル、凡例、マーカーサイズ、色などのプロット設定は`cparams`から取得されます。
引数:
:param app: tkApplicationのインスタンス。
:param cparams: tkParamsのインスタンス。グラフの描画設定を保持。
:param fit: tkFitのインスタンス。データおよびフィッティングに関する情報を保持。
:param ax: matplotlib.axes.Axes プロット対象のMatplotlib Axesオブジェクト。
:param xT: numpy.ndarray 温度データの配列 (K)。
:param ywmugb: numpy.ndarray 粒界散乱の重みデータの配列。
:param ywmuop: numpy.ndarray 光学的フォノン散乱の重みデータの配列。
:param ywmuT0: numpy.ndarray 中性不純物散乱の重みデータの配列。
:param ywmuT32: numpy.ndarray 格子振動散乱の重みデータの配列。
:param xlabel: str, optional X軸のラベル。デフォルトは"T (K)"。
:param ylabel: str, optional Y軸のラベル。デフォルトは"Linear weight"。
:param colors: list, optional プロットに使用する色のリスト。デフォルトはNone。
:param markersize: float, optional データポイントのマーカーサイズ。デフォルトは1.0。
:param plot_event: tkPlotEvent, optional プロットイベントハンドラ。プロットデータ登録に使用。デフォルトはNone。
戻り値:
None
"""
color = colors[0]
wgb_data = ax.plot(xT, ywmugb, label = '$w_{GB}$', linewidth = 1.0, linestyle = 'dashed', color = color)
color = colors[1]
wop_data = ax.plot(xT, ywmuop, label = '$w_{op}$', linewidth = 1.0, linestyle = 'dashed', color = color,
marker = 'o', markersize = 2.0)
color = colors[2]
wT0_data = ax.plot(xT, ywmuT0, label = '$w_{T0}$', linewidth = 1.0, linestyle = 'dashed', color = color,
marker = 'o', markersize = 2.0)
color = colors[3]
wT32_data = ax.plot(xT, ywmuT32, label = '$w_{T3/2}$', linewidth = 1.0, linestyle = 'dashed', color = color,
marker = 'o', markersize = 2.0)
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})
plot_event.add_data({"label": "w_T0", "plot_type": "plot", "axis": ax, "data": wT0_data})
plot_event.add_data({"label": "w_T32", "plot_type": "plot", "axis": ax, "data": wT32_data})
ax.set_xlabel(xlabel, fontsize = cparams.fontsize)
ax.set_ylabel(ylabel, fontsize = cparams.fontsize)
ax.legend(fontsize = cparams.legend_fontsize)
[ドキュメント]
def fit_mu_T(app):
"""
ホール移動度データに対する非線形最小二乗フィッティングを実行する。
概要:
指定されたデータファイルからホール移動度データを読み込み、モデルに基づいて非線形最小二乗法でフィッティングを行う。
フィッティングの初期化、最適化、結果の表示、パラメータとフィッティングデータの保存、およびグラフ表示を行う。
詳細説明:
この関数は以下の主要なステップを実行します:
1. フィッティング設定のログ出力。
2. `tkFit`オブジェクトの初期化とデータファイルの読み込み。
3. フィッティング対象となるパラメータ（varname, unit, pk, optidなど）の設定。
4. パラメータファイルからの初期パラメータの読み込み。
5. フィッティング関数として`cal_mu_from_params`を設定。
6. 初期フィッティング結果の計算と表示。
7. 初期状態のグラフ描画。
8. `scipy.optimize.minimize`を使用した非線形最適化の実行。
9. 最適化された最終パラメータとフィッティングスコアの表示。
10. 各散乱機構からの移動度成分と、それらの重み（寄与度）の計算と表示。
11. パラメータファイル、フィッティングデータ、成分データをExcelファイルに保存。
12. 最終的なフィッティング結果と分解された移動度成分、重みをグラフとして表示。
引数:
:param app: tkApplicationのインスタンス。アプリケーション全体の制御とパラメータを保持。
戻り値:
None
"""
cparams = app.cparams
mu_params = app.mu_params
print("")
print("Fitting to mu-T Hall data by nonlinear least-squares method")
print("mode : ", cparams.mode)
print("infile : ", cparams.infile)
print(f"T range for fitting : {cparams.Tmin} - {cparams.Tmax} K")
print(f"T range for simulation: {cparams.Tcalmin} - {cparams.Tcalmax} K")
print("parameter file : ", cparams.parameterfile)
print("xlabel : ", cparams.xlabel)
print("ylabel : ", cparams.ylabel)
print("method : ", cparams.method)
print("tol : ", cparams.tol)
print("nmaxiter: ", cparams.nmaxiter)
print("h_diff : ", cparams.h_diff)
if '***' in cparams.method:
app.terminate("\nError: Choose method", pause = True)
if '***' in cparams.xlabel:
app.terminate("\nError: Choose xlabel", pause = True)
if '***' in cparams.ylabel:
app.terminate("\nError: Choose ylabel", pause = True)
fit = tkFit(tol = cparams.tol, nmaxiter = cparams.nmaxiter,
print_interval = cparams.print_interval, plot_interval = cparams.plot_interval)
print("")
print(f"Read data from [{cparams.infile}]")
fit.read_data(cparams.infile, xlabel = cparams.xlabel, ylabel = cparams.ylabel,
xmin = cparams.Tmin, xmax = cparams.Tmax, usage = usage)
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 = [ "Eop", "Aop", "AT0", "AT32", "VB", "s_phi"]
fit.unit = [ "", "", "", "", "eV", "eV"]
fit.pk = [mu_params.Eop, mu_params.Aop, mu_params.AT0, mu_params.AT32, mu_params.VB, mu_params.s_phi]
fit.optid = [ 0, 1, 1, 1, 1, 0]
fit.optid_llsq = [ 0, 1, 1, 1, 0, 0]
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)
method = cparams.method
optid_original = fit.optid.copy()
if cparams.mode == 'init':
method = 'cg'
for i in range(nvars):
if fit.optid_llsq[i] == 0:
fit.optid[i] = 0
fit.read_parameters(cparams.parameterfile, section = "mu",
keys = fit.varname, values = fit.pk, optid = fit.optid)
fit.func = lambda T, *pk: cal_mu_from_params(T, *pk)
# print("")
# fit.print_attributes()
# exit()
yini = fit.cal_ylist(fit.pk)
optpk = fit.extract_parameters()
fini = fit.minimize_func(optpk, x = fit.x, y = fit.y)
print("")
print(f"Initial function: fmin={fini:10.4g}")
print(f"{'i':4} {'T':10} {'mu':12} {'mu(cal)':12}")
for i in range(fit.ndata):
print(f"{i:4} {fit.x[i]:10.4g} {fit.y[i]:12.4g} {yini[i]:12.4g}")
#========================================
# plot
#========================================
fig, axes = plt.subplots(1, 3, figsize = cparams.figsize)
# axes = [axes]
# axes = axes.flatten()
title = app.replace_path(cparams.infile, template = "{filename}")
plt.title(title, fontsize = cparams.fontsize)
fit.initial_plot(data_axis = axes[0], error_axis = axes[1], yini = yini, label_ini = 'initial',
fmin = fini, fig = fig,
fontsize = cparams.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("")
print(f"Converged at iteration: {res.nit}")
else:
print(f"Function did not converge")
mu_params.Eop, mu_params.Aop, mu_params.AT0, mu_params.AT32, mu_params.VB, mu_params.s_phi = [*pfin]
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}")
# Recover optid to save to parameter file
if cparams.mode == 'init':
fit.optid = optid_original
#========================================
# Final result
#========================================
yfin = fit.cal_ylist(pfin)
print("")
print("Final result")
# 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} {yini[i]:12.4g} {yfin[i]:12.4g}")
fit.print_scores(heading = "\nScores between y(input) and y(fit)", y1 = fit.y, y2 = yfin)
#=============================
# 貢献度の計算
#=============================
yKgb = []
ymu_ingrain = []
ymuop = []
ymuT0 = []
ymuT32 = []
print("")
print("Components of mobility")
print(f"{'T(K)':6} {'mu,obs(cm2/Vs)':10} {'mu,tot':10} {'Kgb':10} {'mu,in-grain':10} {'mu,op':10} "
+ f"{'mu,T0':10} {'mu,T3/2':10}")
for iT in range(len(fit.x)):
T = fit.x[iT]
mu_tot, mu_KGB, mu_ingrain, mu_op, mu_T0, mu_T32 = cal_mu_components(T,
mu_params.Eop, mu_params.Aop, mu_params.AT0, mu_params.AT32,
mu_params.VB, mu_params.s_phi)
yKgb.append(mu_KGB)
ymu_ingrain.append(mu_ingrain)
ymuop.append(mu_op)
ymuT0.append(mu_T0)
ymuT32.append(mu_T32)
print(f"{T:6.3g} {fit.y[iT]:10.4g} {mu_tot:10.4g} {mu_KGB:10.4g} {mu_ingrain:10.4g} "
+ f"{mu_op:10.4g} {mu_T0:10.4g} {mu_T32:10.4g}")
ywmugb = []
ywmuop = []
ywmuT0 = []
ywmuT32 = []
print("")
print(f"{'T(K)':6} {'w,gb':10} {'w,op':10} {'w,T0':10} {'w,T32':10}")
for iT in range(len(fit.x)):
T = fit.x[iT]
totinv = 1.0 / abs(ymuop[iT]) + 1.0 / abs(ymuT0[iT]) + 1.0 / abs(ymuT32[iT])
ywmugb.append(1.0 - yKgb[iT])
ywmuop.append(yKgb[iT] / ymuop[iT] / totinv)
ywmuT0.append(yKgb[iT] / ymuT0[iT] / totinv)
ywmuT32.append(yKgb[iT] / ymuT32[iT] / totinv)
print(f"{T:6.3g} {ywmugb[iT]:10.4g} {ywmuop[iT]:10.4g} {ywmuT0[iT]:10.4g} {ywmuT32[iT]:10.4g}")
print("")
print(f"Save parameters to {cparams.parameterfile}")
cparams.save_parameters (cparams.parameterfile, section = 'Preferences', sort_by_keys = False, update_commandline = False, IsPrint = False)
# mu_params.save_parameters(cparams.parameterfile, section = 'mu', sort_by_keys = False, update_commandline = False, IsPrint = False)
fit.save_parameters (cparams.parameterfile, section = 'mu', heading = f"Save parameters to {cparams.parameterfile}")
# , keys = fit.varname, values = fit.pk, optid = fit.optid)
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}-muT-components.xlsx"])
# outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-muT-fit.xlsx"])
# print("")
# print(f"Save results to [{outfile}]")
# fit.to_excel(outfile, ['T(K)', 'mu', 'mu(ini)', 'mu(fin)'], [fit.x, fit.y, yini, yfin])
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)"],
[fit.x, fit.y, yini, yfin])
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", "mu,T0", "mu,T32", "w,gb", "w,op", "w,T0", "w,T32"],
[fit.x, fit.y, yfin, yKgb, ymu_ingrain, ymuop, ymuT0, ymuT32, ywmugb, ywmuop, ywmuT0, ywmuT32])
#=============================
# グラフの表示
#=============================
print("")
print("Plot optimized")
ymin = min(fit.y)
yinmax = max(fit.y)
ymax = yinmax
for m in [ymu_ingrain, ymuop, ymuT0, ymuT32]:
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
axes[0].tick_params(labelsize = cparams.fontsize)
axes[1].tick_params(labelsize = cparams.fontsize)
axes[2].tick_params(labelsize = cparams.fontsize)
fit.finalize_plot(yfin, iter = res.nit, fmin = ffin)
plot_muT_decomposed(app, cparams, fit, axes[1], fit.x, fit.y, yfin, ymu_ingrain, ymuop, ymuT0, ymuT32,
colors = cparams.colors, markersize = 3.0, plot_event = fit.plot_event)
axes[1].set_ylim([ymin * 0.8, ymax * 1.1])
plot_muT_weight(app, cparams, fit, axes[2], fit.x, ywmugb, ywmuop, ywmuT0, ywmuT32,
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 sim(app):
"""
ホール移動度の温度依存性をシミュレーションし、観測データと比較する。
概要:
パラメータファイルから読み込んだモデルパラメータを使用して、ホール移動度の温度依存性をシミュレーションする。
観測データとの比較を行い、結果をファイルに保存し、温度および1000/Tに対してプロットする。
詳細説明:
この関数は以下の主要なステップを実行します:
1. シミュレーション設定のログ出力。
2. `tkFit`オブジェクトの初期化とデータファイルの読み込み。
3. シミュレーションに使用するパラメータ（varname, unit, pk, optidなど）の設定。
4. パラメータファイルからのパラメータの読み込み。
5. シミュレーション関数として`cal_mu_from_params`を設定。
6. シミュレーション結果の計算と表示。
7. 結果データ（温度、1000/T、観測移動度、シミュレーション移動度）をExcelファイルに保存。
8. 観測データとシミュレーション結果を比較するグラフ（μ-Tとμ-1000/T）を表示。
引数:
:param app: tkApplicationのインスタンス。アプリケーション全体の制御とパラメータを保持。
戻り値:
None
"""
cparams = app.cparams
mu_params = app.mu_params
cparams.outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-simulate.xlsx"])
print("")
print("Simulate mu-T dependence and compare with Hall data")
print("mode : ", cparams.mode)
print("infile : ", cparams.infile)
print(f"T range for simulation: {cparams.Tcalmin} - {cparams.Tcalmax} K")
print("outfile: ", cparams.outfile)
print("xlabel : ", cparams.xlabel)
print("ylabel : ", cparams.ylabel)
fit = tkFit(tol = cparams.tol, nmaxiter = cparams.nmaxiter,
print_interval = cparams.print_interval, plot_interval = cparams.plot_interval)
print("")
print(f"Read data from [{cparams.infile}]")
fit.read_data(cparams.infile, xlabel = cparams.xlabel, ylabel = cparams.ylabel,
xmin = cparams.Tmin, xmax = cparams.Tmax, usage = usage)
fit.varname = [ "Eop", "Aop", "AT0", "AT32", "VB", "s_phi"]
fit.unit = [ "", "", "", "", "eV", "eV"]
fit.pk = [mu_params.Eop, mu_params.Aop, mu_params.AT0, mu_params.AT32, mu_params.VB, mu_params.s_phi]
fit.optid = [ 0, 1, 1, 1, 1, 1]
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.0e4, 1.0, 1.0e2, 1.0e4]
fit.read_parameters(cparams.parameterfile, section = "mu",
keys = fit.varname, values = fit.pk, optid = fit.optid)
fit.func = lambda T, *pk: cal_mu_from_params(T, *pk)
print("")
print("Paramters:")
print(" # of variables:", len(fit.pk))
fit.print_variables()
print("")
print("Input data:")
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}")
yini = fit.cal_ylist(fit.pk)
Tinv = [1000.0 / fit.x[i] for i in range(fit.ndata)]
optpk = fit.extract_parameters()
fini = fit.minimize_func(optpk)
print("")
print(f"Simulation result: fmin={fini:10.4g}")
print(f"{'i':4} {'T':10} {'1000/T':10} {'mu':12} {'mu(cal)':12}")
for i in range(fit.ndata):
print(f"{i:4} {fit.x[i]:10.4g} {Tinv[i]:10.4g} {fit.y[i]:12.4g} {yini[i]:12.4g}")
print("")
print(f"Save results to [{cparams.outfile}]")
fit.to_excel(cparams.outfile, ['T(K)', '1000/T(K-1)', 'mu', 'mu(sim)'], [fit.x, Tinv, fit.y, yini])
#========================================
# plot
#========================================
fig, axes = plt.subplots(1, 2, figsize = cparams.figsize)
# axes = [axes]
# axes = axes.flatten()
axis = axes[0]
ax2 = axes[1]
plot_event = tkPlotEvent(plt)
# fit.initial_plot(data_axis = axes[0], error_axis = axes[1], yini = yini, label_ini = 'initial',
# fmin = fini, fig = fig,
# fontsize = cparams.fontsize)
axis.tick_params(labelsize = cparams.fontsize)
ax2.tick_params(labelsize = cparams.fontsize)
input_data = axis.plot(fit.x, fit.y, label = "input", linestyle = '', marker = 'o', markersize = 5.0)
sim_data = axis.plot(fit.x, yini, label = "simulate", linestyle = '-', linewidth = cparams.linewidth, color = 'red')
axis.set_xlabel(fit.xlabel, fontsize = cparams.fontsize)
axis.set_ylabel(fit.ylabel, fontsize = cparams.fontsize)
axis.legend(fontsize = cparams.fontsize)
plot_event.add_data({"label": "input", "plot_type": "plot", "axis": axis, "data": input_data})
plot_event.add_data({"label": "simulate", "plot_type": "plot", "axis": axis, "data": sim_data})
input_Arr_data = ax2.plot(Tinv, fit.y, label = "input", linestyle = '', marker = 'o', markersize = 5.0)
sim_Arr_data = ax2.plot(Tinv, yini, label = "simulate", linestyle = '-', linewidth = cparams.linewidth, color = 'red')
ax2.set_xlabel("1000/T (K$^{-1}$)", fontsize = cparams.fontsize)
ax2.set_ylabel(fit.ylabel, fontsize = cparams.fontsize)
ax2.set_yscale('log')
ax2.legend(fontsize = cparams.fontsize)
plot_event.add_data({"label": "input", "plot_type": "plot", "axis": axis, "data": input_Arr_data})
plot_event.add_data({"label": "simulate", "plot_type": "plot", "axis": axis, "data": sim_Arr_data})
plot_event.register_event(fig)
plt.tight_layout()
plt.pause(0.001)
print("")
# print(f"Save parameters to {cparams.parameterfile}")
cparams.save_parameters (cparams.parameterfile, section = 'Preferences', sort_by_keys = False, update_commandline = False, IsPrint = False)
# mu_params.save_parameters(cparams.parameterfile, section = 'mu', sort_by_keys = False, update_commandline = False, IsPrint = False)
fit.save_parameters (cparams.parameterfile, section = 'mu', heading = f"Save parameters to {cparams.parameterfile}")
# , keys = fit.varname, values = fit.pk, optid = fit.optid)
app.terminate("", usage = usage, pause = True)
[ドキュメント]
def main():
"""
スクリプトのメインエントリポイント。
概要:
アプリケーションとパラメータの初期化、コマンドライン引数によるパラメータの更新、
ログファイル設定、パラメータファイル読み込み、および指定されたモードに応じた処理を実行する。
詳細説明:
この関数は、プログラムの開始点として以下の処理を行います。
1. `initialize`関数を呼び出して、`tkApplication`と`tkParams`オブジェクトを準備します。
2. コマンドライン引数を解析し、`cparams`の初期設定を上書きします。
3. ログファイルのパスを設定し、標準出力をログファイルにもリダイレクトします。
4. パラメータファイルのパスを決定し、既存のパラメータがあれば読み込みます。
5. 再度コマンドライン引数を解析し、読み込んだパラメータをさらに上書きします。
6. `cparams`の現在のパラメータ設定をコンソールに出力します。
7. `cparams.mode`の値に応じて、`fit_mu_T`関数（フィッティング）または`sim`関数（シミュレーション）を呼び出します。
8. サポートされていないモードが指定された場合はエラーメッセージを表示します。
9. 処理の終了時にアプリケーションを終了し、一時停止します。
戻り値:
None
"""
app, cparams = initialize()
update_vars(app, cparams)
cparams.logfile = app.replace_path(cparams.infile)
cparams.outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-out.xlsx"])
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"])
# app.mu_params.read_parameters(cparams.parameterfile, section = "mu",
# ignore_keys = ["logfile", "outfile", "parameterfile", "parameterbkfile"])
# check args again to use given parameters
update_vars(app, cparams)
cparams.print_parameters(heading = "\ncparams:")
if cparams.mode == 'init':
fit_mu_T(app)
elif cparams.mode == 'fit':
fit_mu_T(app)
elif cparams.mode == 'sim':
sim(app)
else:
app.terminate("Error in main: Invalide mode [{}]".format(cparams.mode), usage = usage, pause = True)
app.terminate("", usage = usage, pause = True)
if __name__ == "__main__":
main()