"""
アレニウスプロット解析ツール
概要:
このスクリプトは、温度依存性データ（P vs T）をアレニウスプロットに変換し、活性化エネルギーなどの物理量を解析します。
指定されたモデル（例: シンプルアレニウス、パーコレーションなど）に基づいて多項式フィッティングを行い、
結果をプロットしてExcelファイルに保存します。
詳細説明:
入力ファイルから温度と特性値のデータを読み込み、ユーザーが指定したデータ形式（T(K), T(C), 1/T, 1000/T, P, log10(P), log_e(P)）に変換します。
フィッティング範囲を指定してデータを抽出し、1次以上の多項式フィッティング（通常は1000/Tに対するlog10(P)）を実行します。
フィッティング結果から、プレファクターP0、活性化エネルギーEa、およびモデルによってはその分布の標準偏差σ_phiなどを算出します。
計算されたデータとフィッティングカーブを複数のグラフに表示し、イベントハンドラでインタラクティブな操作を可能にします。
すべての結果はログファイルとExcelファイルに出力されます。
関連リンク:
:doc:`arrhenius_plot_usage`
"""
import sys
import os
import signal
import builtins
import numpy as np
from numpy import exp, log, sin, cos, tan, arcsin, arccos, arctan, sinh, cosh, tanh, sqrt, abs
import pandas as pd
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from tklib.tkutils import replace_path, getarg, getintarg, getfloatarg, pint, pfloat
from tklib.tksci.tksci import kB, e
from tklib.tkvariousdata import tkVariousData
from tklib.tkapplication import tkApplication
from tklib.tkparams import tkParams
from tklib.tksci.tksci import asin, acos, atan, degcos, degsin, degtan, degacos, degasin, degatan, arcsin, arccos, arctan
from tklib.tksci.tksci import factorial, log10, gamma, combination, eVTonm, nmToeV, Bn
from tklib.tksci.tksci import h, h_bar, hbar, e, kB, NA, c, pi, pi2, torad, todeg, basee
from tklib.tksci.tksci import me, mp, mn
from tklib.tksci.tksci import u0, e0, a0, R, F, g, G
from tklib.tksci.tksci import HartreeToeV, RyToeV, KToeV, eVToK, JToeV, eVToJ, Debye
from tklib.tksci.tkFit import tkFit
from tklib.tkgraphic.tkplotevent import tkPlotEvent
from tklib.tksci.tkFit import tkFit
#=======================================
# プログラム開始
#=======================================
[ドキュメント]
def initialize(app):
"""
アプリケーションの初期パラメータを設定します。
概要:
`tkParams` オブジェクトを初期化し、解析に必要な各種設定値のデフォルト値を設定します。
詳細説明:
入力ファイル名、解析モデル、データラベル、データタイプ、フィッティングおよび計算の範囲、
プロットの図サイズとフォントサイズなどのデフォルト値を設定し、アプリケーションの `cparams` 属性に格納します。
:param app: アプリケーションインスタンス。
:type app: tkApplication
:returns: なし
:rtype: None
"""
cparams = tkParams()
app.cparams = cparams
cparams.infile = ''
cparams.model = 'simple Arrhenius'
cparams.Tlabel = 'T(K)'
cparams.Plabel = 'P'
cparams.Ttype = 'T(K)'
cparams.Ptype = 'P'
cparams.xmin = -1.0e100
cparams.xmax = 1.0e100
cparams.Tmin = -1.0e100
cparams.Tmax = 1.0e100
cparams.Tcalmin = '*'
cparams.Tcalmax = '*'
cparams.figsize = [8, 8]
cparams.fontsize = 12
cparams.legend_fontsize = 10
[ドキュメント]
def update_vars(app):
"""
コマンドライン引数に基づいてアプリケーションのパラメータを更新します。
概要:
コマンドライン引数パーサーに引数を登録し、実際にコマンドラインから与えられた値でパラメータを更新します。
詳細説明:
入力ファイル名、解析モデル、T軸・P軸のデータラベル、T軸・P軸のデータタイプ、
フィッティングのX軸・T軸範囲、計算のT軸範囲、プロットのサイズなどの引数を `app.add_argument` を用いて定義します。
その後、`app.read_args` でこれらの引数を解析し、エラーが発生した場合はプログラムを終了します。
:param app: アプリケーションインスタンス。
:type app: tkApplication
:returns: なし
:rtype: None
"""
cparams = app.cparams
app.add_argument(opt = None, type = "str", var_name = 'infile', opt_str = "infile", desc = 'Input file',
defval = "Hall-T.xlsx", optional = True);
app.add_argument(opt = "-model", type = "str", var_name = 'model', opt_str = "-model='simple Arrhenius'",
desc = 'model: [simple Arrhenius|percolation]',
defval = "simple Arrhenius", optional = True);
app.add_argument(opt = "-Tlabel", type = "str", var_name = 'Tlabel', opt_str = "-Tlabel=T(K)",
desc = 'Label of T-related data in the input file',
defval = "T(K)", optional = True);
app.add_argument(opt = "-Plabel", type = "str", var_name = 'Plabel', opt_str = "-Plabel=P",
desc = 'Label of property-related data in the input file]',
defval = "P", optional = True);
app.add_argument(opt = "-Ttype", type = "str", var_name = 'Ttype', opt_str = "-Ttype=T(K)",
desc = 'Type of T-related data: [T(K)|T(C)|1/T|1000/T]',
defval = "T(K)", optional = True);
app.add_argument(opt = "-Ptype", type = "str", var_name = 'Ptype', opt_str = "-Ptype=P",
desc = 'Type of property-related data: [P|log10(P)|log_e(P)]',
defval = "P", optional = True);
app.add_argument(opt = "-xmin", type = "double", var_name = 'xmin', opt_str = "-xmin=-1.0e-100",
desc = 'Lower limit of T-related data (x) for fitting', defval = "-1.0e-100", optional = True);
app.add_argument(opt = "-xmax", type = "double", var_name = 'xmax', opt_str = "-xmax=1.0e-100",
desc = 'Upper limit of T-related data (x) for fitting', defval = "1.0e-100", optional = True);
app.add_argument(opt = "-Tmin", type = "double", var_name = 'Tmin', opt_str = "-Tmin=-1.0e-100",
desc = 'Lower limit of T for fitting', defval = "-1.0e-100", optional = True);
app.add_argument(opt = "-Tmax", type = "double", var_name = 'Tmax', opt_str = "-Tmax=1.0e-100",
desc = 'Upper limit of T for fitting', defval = "1.0e-100", optional = True);
app.add_argument(opt = "-Tcalmin", type = "str", var_name = 'Tcalmin', opt_str = "-Tcalmin='*'",
desc = 'Lower limit of T for calculation', defval = "*", optional = True);
app.add_argument(opt = "-Tcalmax", type = "str", var_name = 'Tcalmax', opt_str = "-Tcalmax='*'",
desc = 'Upper limit of T for calculation', defval = "*", optional = True);
args_opt, args_idx, args_vars = app.read_args(vars = cparams.__dict__, check_allowed_args = True)
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, pause = True)
# app.set_usage(usage_str)
[ドキュメント]
def execute(app):
"""
アレニウスプロット解析の主要な処理を実行します。
概要:
入力データの読み込み、データ変換、指定範囲でのフィッティング、結果計算、Excelファイルへの出力、およびグラフ描画を行います。
詳細説明:
1. ログファイル、出力ファイル、フィッティング結果ファイルのパスを設定し、標準出力をログファイルにリダイレクトします。
2. コマンドライン引数で受け取った数値パラメータ（xmin, xmax, Tmin, Tmax）を浮動小数点数に変換します。
3. 指定された入力ファイルからデータを読み込み、`Ttype` と `Ptype` の設定に従ってデータ（TとP）を変換します。
例えば、`Ttype='1000/T'` であれば1000/Tの形式に変換し、`Ptype='log10(P)'` であればlog10(P)の形式に変換します。
4. `xmin`, `xmax`, `Tmin`, `Tmax` の範囲内でフィッティング対象のデータ点を選別します。
5. `model` パラメータに応じて多項式の次数（1次から4次）を決定し、`numpy.polyfit` を用いて最小二乗フィッティングを実行します。
6. フィッティングの結果得られた多項式係数から、P0、活性化エネルギーEa、および2次以上のモデルの場合にはσ_phiを計算し出力します。
7. `tkFit` を使用して、元のデータとフィッティング結果の間のスコア（RMSDなど）を計算し表示します。
8. `Tcalmin` から `Tcalmax` までの範囲で、フィッティング結果に基づいてP(cal)、log10(P)(cal)、およびEa(T)を計算します。
9. 元のデータ、変換されたデータ、フィッティング結果、計算されたデータを含む結果をExcelファイルに出力します。
10. Matplotlibを用いて、以下の4つのプロットを含む図を作成し表示します。
- 元のX-Yプロット
- T-Pプロットとフィッティングカーブ
- アレニウスプロット (1000/T vs log10(P)) とフィッティングカーブ
- 1000/Tに対する活性化エネルギーEaのプロット
11. `tkPlotEvent` を使用して、プロットにマウスイベントハンドラを登録します。
12. プログラムを一時停止し、ユーザーの終了待ち状態にします。
:param app: アプリケーションインスタンス。
:type app: tkApplication
:returns: なし
:rtype: None
"""
cparams = app.cparams
cparams.logfile = app.replace_path(cparams.infile)
cparams.outfile = app.replace_path(cparams.infile, "{dirname}/{filebody}-out.xlsx")
cparams.outfitfile = app.replace_path(cparams.infile, "{dirname}/{filebody}-fit.xlsx")
print(f"Open logfile [{cparams.logfile}]")
app.redict(targets = ["stdout", cparams.logfile], mode = 'w')
cparams.xmin = pfloat(cparams.xmin)
cparams.xmax = pfloat(cparams.xmax)
cparams.Tmin = pfloat(cparams.Tmin)
cparams.Tmax = pfloat(cparams.Tmax)
print("")
print("#======================================================")
print("# Analyze activation energy etc by Arrhenius plot")
print("#======================================================")
print(f"infile : {cparams.infile}")
print(f"outfile: {cparams.outfile}")
print(f"logfile: {cparams.logfile}")
print(f"model : {cparams.model}")
print(f"Tlabel: {cparams.Tlabel}")
print(f"Plabel: {cparams.Plabel}")
print(f"Ttype : {cparams.Ttype}")
print(f"Ptype : {cparams.Ptype}")
print(f"Fitting x range: {cparams.xmin} - {cparams.xmax}")
print(f"Fitting T range: {cparams.Tmin} - {cparams.Tmax}")
if '***' in cparams.model:
app.terminate("Error: Choose model", pause = True)
print(f"")
print(f"Read [{cparams.infile}]")
datafile = tkVariousData(cparams.infile)
labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = app.usage)
label_x, xX = datafile.FindDataArray(cparams.Tlabel, flag = 'i')
label_y, yY = datafile.FindDataArray(cparams.Plabel, flag = 'i')
# print("X=", label_x, xX)
# print("Y=", label_y, yY)
if cparams.Ttype == 'T(K)':
T = xX
elif cparams.Ttype == 'T(C)':
T = [T + 273.15 for T in xX]
elif cparams.Ttype == '1/T':
T = [1.0 / T for T in xX]
elif cparams.Ttype == '1000/T':
T = [1000.0 / T for T in xX]
else:
app.terminate(f"\nInvalid Ttype [{cparams.Ttype}]", usage = app.usage, pause = True)
if cparams.Ptype == 'P':
P = yY
elif cparams.Ptype == 'log10(P)':
P = [np.power(10.0, P) for P in yY]
elif Ttype == 'log_e(P)':
P = [exp(P) for P in yY]
else:
app.terminate(f"\nInvalid Ptype [{cparams.Ptype}]", usage = app.usage, pause = True)
print("T(K)=", T)
print("P =", P)
T1000 = [1000.0 / v for v in T]
log10P = [log(v) / log(10.0) for v in P]
# print("1000/T (K^-1)=", T1000)
# print("log10(P) =", log10P)
x_fit = []
y_fit = []
for i in range(len(xX)):
if cparams.xmin <= xX[i] <= cparams.xmax and cparams.Tmin <= T[i] <= cparams.Tmax:
x_fit.append(T1000[i])
y_fit.append(log10P[i])
if cparams.model == 'percolation':
norder = 2
elif cparams.model == '3rd order':
norder = 3
elif cparams.model == '4th order':
norder = 4
else:
norder = 1
print("")
print(f"Least-squares fitting with {norder}-th order polynomial of 1000/T")
print("Data to be fitted:")
print(f" {'1000/T':12} {'log10(P)':12}")
for i in range(len(x_fit)):
print(f" {x_fit[i]:12.4g} {y_fit[i]:12.4g}")
ci = np.polyfit(x_fit, y_fit, norder)
print("Coefficients:")
for i in range(norder + 1):
idx = norder - i
print(f" c{i}={ci[idx]:12.4g}")
log10Pfit = np.poly1d(ci)(x_fit)
log10Pcal = np.poly1d(ci)(T1000)
# P(T) = P0 * exp[-e * Ea / kT + e^2 * s0^2 / 2 / kB^2 / T^2]
# log10P(T) = log10(P0) - log_10(e) * e / kB / 1000.0 * Ea * (1000 / T) + log_10(e) * e^2 / kB^2 / 1000^2 * s0^2 / 2 * (1000 / T)^2
# c1 = -log10(e) * e / kB / 1000 * Ea: Ea = -c1 * ln(10) * kB / e * 1000.0
# c2 = log10(e) * (e / kB / 1000)^2 /2 * s0^2: s0 = sqrt(c2 * ln(10) * (1000 kB/e)^2 * 2)
P0 = np.power(10.0, ci[norder])
Ea = -ci[norder-1] * kB / e * 1000.0 * log(10.0)
print("")
print(f" P0={P0:12.4g}")
print(f" Ea=phi_0={Ea:12.4g} eV")
if norder >= 2:
sigma_phi = ci[norder-2] * (1000.0 * kB / e)**2 * 2.0 * log(10.0)
if sigma_phi < 0.0:
print(f"\n ***Warning: sigma_phi^2 is negative: {sigma_phi:12.4g} eV^2")
else:
sigma_phi = sqrt(sigma_phi)
print(f" sigma_phi={sigma_phi:12.4g} eV")
fit = tkFit()
P_fit = [pow(10.0, log10Pfit[i]) for i in range(len(x_fit))]
fit.print_scores(heading = "\nScores between P(input) and P(fit)", y1 = P, y2 = P_fit)
fit.print_scores(heading = "\nScores between log10(P(input)) and log10(P(fit))", y1 = log10P, y2 = log10Pfit)
print("")
print( "Calculate data from the fitting result")
Tcalmin = pfloat(cparams.Tcalmin, defval = min(T))
Tcalmax = pfloat(cparams.Tcalmax, defval = max(T))
nT = 101
Tstep = (Tcalmax - Tcalmin) / (nT - 1)
print(f" T range: {Tcalmin:8.3f} - {Tcalmax:8.3f} K")
T_plot = [Tcalmin + i * Tstep for i in range(nT)]
x_plot = [1000.0 / T_plot[i] for i in range(nT)]
y_plot = np.poly1d(ci)(x_plot)
P_plot = [pow(10.0, y_plot[i]) for i in range(nT)]
diff_plot = [(y_plot[1] - y_plot[0]) / (x_plot[1] - x_plot[0])]
for i in range(1, nT-1):
diff_plot.append((y_plot[i+1] - y_plot[i-1]) / (x_plot[i+1] - x_plot[i-1]))
diff_plot.append((y_plot[nT-1] - y_plot[nT-2]) / (x_plot[nT-1] - x_plot[nT-1]))
Ea_plot = [-diff_plot[i] * kB / e * 1000.0 * log(10.0) for i in range(nT)]
"""
print("")
print("Data")
print(f"{'idx':5} {'X':12} {'Y':12} {'T':8} {'P':12} {'1000/X':12} {'log10(P)':12} {'log10(P)(cal)':12}")
for i in range(len(xX)):
print(f"{i:5} {xX[i]:12.4g} {yY[i]:12.4g} {T[i]:8.3g} {P[i]:12.4g} {T1000[i]:12.4g} {log10P[i]:12.4g} {log10Pcal[i]:12.4g}")
"""
print("")
print(f"Save to [{cparams.outfile}]")
fit.to_excel(cparams.outfile, [label_x, label_y, 'T(K)', 'P', '1000/T (K^-1)', 'log10(P)', 'log10(P)(cal)',
'', 'T (K)', '1000/T (K^-1)', 'P(cal)', 'log10(P)(cal)'],
[xX, yY, T, P, T1000, log10P, log10Pcal, [], T_plot, x_plot, P_plot, y_plot])
#=========================================
# plot
#=========================================
fig, axes = plt.subplots(2, 2, figsize = cparams.figsize)
plot_event = tkPlotEvent(plt)
axes = axes.flatten()
axes[0].tick_params(labelsize = cparams.fontsize)
axes[1].tick_params(labelsize = cparams.fontsize)
axes[2].tick_params(labelsize = cparams.fontsize)
axes[3].tick_params(labelsize = cparams.fontsize)
XY_data = axes[0].plot(xX, yY, linestyle = '', marker = 'o', markerfacecolor = 'black', markersize = 5.0)
axes[0].set_xlabel(label_x, fontsize = cparams.fontsize)
axes[0].set_ylabel(label_y, fontsize = cparams.fontsize)
TP_data = axes[1].plot(T, P, linestyle = '', marker = 'o', markerfacecolor = 'black', markersize = 5.0)
TP_fit_data = axes[1].plot(T_plot, P_plot, linestyle = '-', color = 'red', linewidth = 0.5)
axes[1].set_xlabel('$T$ (K)', fontsize = cparams.fontsize)
axes[1].set_ylabel('$P$', fontsize = cparams.fontsize)
arr_data = axes[2].plot(T1000, log10P, label = 'raw data', linestyle = '', marker = 'o', markerfacecolor = 'black', markersize = 5.0)
fit_data = axes[2].plot(x_fit, y_fit, label = 'raw data (fitted)', linestyle = '', marker = 'o', markerfacecolor = 'red', markersize = 8.0)
plot_data = axes[2].plot(x_plot, y_plot, label = 'fitted', color = 'red', linewidth = 0.5)
axes[2].set_xlabel('$1000/T$ (K$^{-1}$)', fontsize = cparams.fontsize)
axes[2].set_ylabel('$log_{10}(P)$', fontsize = cparams.fontsize)
axes[2].legend(fontsize = cparams.legend_fontsize)
Ea_data = axes[3].plot(x_plot, Ea_plot, label = 'Ea (eV)', linestyle = '-', linewidth = 0.5, color = 'black')
axes[3].set_xlabel('$1000/T$ (K$^{-1}$)', fontsize = cparams.fontsize)
axes[3].set_ylabel('$E_a$ (eV)', fontsize = cparams.fontsize)
axes[3].legend(fontsize = cparams.legend_fontsize)
all_data = datalist
plot_event.add_data({"label": "X-Y plot", "plot_type": "2D", "axis": axes[0], "data": XY_data,
"xlist": all_data, "xlabels": labels})
plot_event.add_data({"label": "T-P plot", "plot_type": "2D", "axis": axes[1], "data": TP_data,
"xlist": all_data, "xlabels": labels})
# plot_event.add_data({"label": "fitted", "plot_type": "2D", "axis": axes[1], "data": TP_fit_data,
# "xlist": all_data, "xlabels": labels})
plot_event.add_data({"label": "Arrhenius plot", "plot_type": "2D", "axis": axes[2], "data": arr_data,
"xlist": all_data, "xlabels": labels})
# plot_event.add_data({"label": "fitted", "plot_type": "2D", "axis": axes[2], "data": plot_data})
plot_event.add_data({"label": "Ea", "plot_type": "2D", "axis": axes[3], "data": Ea_data})
plot_event.register_event(fig, event = "button_press_event",
callback = lambda event: plot_event.onclick(event))
plt.tight_layout()
plt.pause(0.001)
# app.terminate("", usage = app.usage, pause = True)
app.terminate("", usage = None, pause = True)
#==========================================
# Main routine
#==========================================
[ドキュメント]
def main():
"""
アプリケーションのエントリポイントです。
概要:
`tkApplication` インスタンスを生成し、初期化、パラメータ更新、メイン処理の実行を行います。
詳細説明:
プログラムの実行フローを管理し、`initialize` 関数でデフォルトパラメータを設定し、
`update_vars` 関数でコマンドライン引数からパラメータを更新し、
`execute` 関数でアレニウスプロット解析の主要なタスクを実行します。
:returns: なし
:rtype: None
"""
app = tkApplication()
print(f"Initialize parameters")
initialize(app)
print(f"Update parameters by command-line arguments")
update_vars(app)
execute(app)
if __name__ == "__main__":
main()