アレニウスプロットによる活性化エネルギー解析ツール

プログラムの動作

arrhenius_plot.py は、アレニウスプロットを用いて、温度依存性のある物理量（物性値）から活性化エネルギー ($E_a$) や前指数因子 ($P_0$) などのパラメータを解析するためのPythonスクリプトです。

主な機能:

データ読み込み: Excel (.xlsx) 形式のファイルから、温度および物性値のデータを読み込みます。
データ変換: 入力された温度データと物性値データを、アレニウスプロットに適した形式（例: $1000/T$ や $log_{10}(P)$）に変換します。変換の種類はコマンドライン引数で選択可能です。
フィッティング: 変換されたデータに対し、最小二乗法を用いて多項式フィッティング（1次、2次、3次、4次）を行います。デフォルトは1次の「simple Arrhenius」モデルです。
パラメータ算出: フィッティング結果から、$P_0$ (前指数因子) および $E_a$ (活性化エネルギー) を計算します。2次フィッティングの場合（「percolation」モデル）、活性化エネルギーのばらつき ($\sigma_\phi$) も算出します。
結果出力: 読み込んだ生データ、変換データ、フィッティング結果、計算されたプロット用データをExcelファイルに出力します。
グラフ表示: 以下の4種類のグラフを生成し、インタラクティブなウィンドウに表示します。
1. 元のX-Yプロット
2. 変換後のT-Pプロットとフィッティング曲線
3. アレニウスプロット ($1000/T$ vs $log_{10}(P)$) とフィッティング曲線
4. 温度に対する活性化エネルギー $E_a$ のプロット
ログ出力: プログラムの実行中に標準出力される情報が、指定されたログファイルに記録されます。

解決する課題:

温度によって挙動が変化する物理現象において、その背後にある活性化プロセスを定量的に評価し、活性化エネルギーなどの重要な物理量を簡便に導出することを目的としています。特に、データの前処理（変換）、フィッティング、結果の解析、および視覚化といった一連の作業を自動化し、研究者や技術者のデータ解析を支援します。

原理

このプログラムは、アレニウスの式とその派生モデルに基づいています。

1. アレニウスの式

物理量 $P$ の温度依存性がアレニウスの式に従う場合、以下の形で表されます。

\[P = P_0 \exp\left(-\frac{E_a}{k_B T}\right)\]

ここで、$P_0$ は前指数因子、$E_a$ は活性化エネルギー、$k_B$ はボルツマン定数、$T$ は絶対温度（ケルビン）です。

この式の両辺を常用対数に変換すると、以下の直線関係が得られます。

\[\log_{10}(P) = \log_{10}(P_0) - \frac{E_a}{k_B T \ln 10}\]

プログラムでは、横軸を $x = 1000/T$ (単位: K$^{-1}$)、縦軸を $y = \log_{10}(P)$ とし、これに対して最小二乗法で多項式フィッティングを行います。

2. フィッティングモデル

プログラムは以下のモデルをサポートします。

simple Arrhenius (1次) $y = c_1 x + c_0$ ここで、$x = 1000/T$, $y = \log_{10}(P)$ です。係数 $c_0$ と $c_1$ から $P_0$ と $E_a$ を算出します。
- 前指数因子 $P_0$: $$P_0 = 10^{c_0}$$
- 活性化エネルギー $E_a$ (単位: eV): $$E_a = -c_1 \cdot \frac{k_B}{e} \cdot \ln 10 \cdot 1000$$ ここで、$k_B$ はボルツマン定数 (J/K)、$e$ は電気素量 (C) です。$k_B/e$ により単位が J から eV に変換されます。
percolation (2次) $y = c_2 x^2 + c_1 x + c_0$ このモデルは、活性化エネルギーにばらつきがある場合などに用いられることがあります。この場合も上記と同様に $P_0$ と $E_a$ が算出されます。さらに、活性化エネルギーのばらつき ($\sigma_\phi$) も計算されます。
- 活性化エネルギーのばらつき $\sigma_\phi$ (単位: eV): $$\sigma_\phi = \sqrt{c_2 \cdot 2 \cdot \ln 10 \cdot \left(\frac{1000 \cdot k_B}{e}\right)^2}$$ 負の値になる場合は警告が表示されます。
3rd order (3次) $y = c_3 x^3 + c_2 x^2 + c_1 x + c_0$
4th order (4次) $y = c_4 x^4 + c_3 x^3 + c_2 x^2 + c_1 x + c_0$

高次フィッティングの場合、$P_0$ と $E_a$ は一次項と定数項から上記と同様に計算されます。

必要な非標準ライブラリとインストール方法

このプログラムを実行するには、以下の非標準ライブラリが必要です。

numpy
pandas
scipy
matplotlib

これらのライブラリは、Pythonのパッケージマネージャー pip を使ってインストールできます。

pip install numpy pandas scipy matplotlib

注記: このプログラムは tklib という独自のライブラリに依存しています。tklib は一般的なPyPIパッケージではないため、上記の pip install コマンドではインストールできません。tklib は、このプログラムが配布される環境で別途提供・配置される必要があります。

必要な入力ファイル

プログラムは、温度と物性値のデータを含むファイルを必要とします。

ファイル形式: Excel (.xlsx) ファイルが推奨されます。tklib.tkVariousData クラスが読み込み可能なCSVなどのテキストファイルも利用できる可能性があります。
データ構造: 少なくとも2列のデータが含まれている必要があります。
- 1列目: 温度に関連するデータ (デフォルトの列ラベルは 'T(K)')
- 2列目: 物性値に関連するデータ (デフォルトの列ラベルは 'P')

これらの列ラベルは、コマンドライン引数で自由に指定できます。

例: data.xlsx

T(K)	P
300	1e-5
320	2e-5
340	5e-5
360	1e-4
380	2e-4
400	4e-4

生成される出力ファイル

プログラムは以下のファイルを生成します。

結果データファイル: {入力ファイル名本体}-out.xlsx
- ファイル名例: data-out.xlsx (入力ファイルが data.xlsx の場合)
- 内容:
  - 元の入力データ (XとY)
  - 変換後の温度 (T(K)) と物性値 (P)
  - アレニウスプロット用の変換データ ($1000/T$ と $log_{10}(P)$)
  - フィッティングモデルから計算された $log_{10}(P)$ の値 (log10(P)(cal))
  - プロット用に計算された滑らかな温度範囲 ($T$, $1000/T$) における物性値 ($P(cal)$) と $log_{10}(P)$ の値 (log10(P)(cal))。
  - 計算された温度範囲における活性化エネルギー ($E_a$(eV))
ログファイル: {入力ファイル名}
- ファイル名例: data.xlsx (入力ファイルが data.xlsx の場合)
- 内容: プログラムの実行中に標準出力されるメッセージ（パラメータ設定、データ読み込み状況、フィッティング結果、計算された$P_0, E_a, \sigma_\phi$ の値、スコアなど）が記録されます。
- 注記: プログラムの現在の実装では、ログファイル名が入力ファイル名と同一となるため、入力ファイルの内容が上書きされて、ログの内容に置き換わる可能性が非常に高いです。 重要な入力ファイルを使用する際は、事前にバックアップを取るか、ログファイルを別のファイル名 (-logfile オプションがあれば設定するなど) に指定するよう注意してください。このプログラムにはログファイル名を指定するオプションが提供されていません。

グラフ画像: プログラムは実行時にグラフウィンドウを表示しますが、自動的に画像ファイルとして保存する機能は含まれていません。表示されたグラフは、手動でスクリーンショットを撮るか、グラフウィンドウの機能を利用して保存する必要があります。

コマンドラインでの使用例 (Usage)

arrhenius_plot.py は、以下の基本的な書式でコマンドラインから実行できます。

python arrhenius_plot.py [infile] [options]

引数:

infile (位置引数、省略可能):
- 入力データファイル名。
- デフォルト値: Hall-T.xlsx
- 例: data.xlsx
-model (オプション引数、文字列):
- フィッティングに用いるモデルを指定します。
- 選択肢: 'simple Arrhenius' (デフォルト), 'percolation', '3rd order', '4th order'
- 例: -model="percolation"
-Tlabel (オプション引数、文字列):
- 入力ファイル内の温度関連データの列ラベルを指定します。
- デフォルト値: 'T(K)'
- 例: -Tlabel="Temperature"
-Plabel (オプション引数、文字列):
- 入力ファイル内の物性値関連データの列ラベルを指定します。
- デフォルト値: 'P'
- 例: -Plabel="Conductivity"
-Ttype (オプション引数、文字列):
- 入力ファイル内の温度データの種類を指定します。データは指定された形式で読み込まれ、その後ケルビン温度 (T(K)) に変換され、さらに $1000/T$ に変換されます。
- 選択肢: 'T(K)' (デフォルト), 'T(C)' (摂氏からケルビンに変換), '1/T', '1000/T'
- 例: -Ttype="T(C)"
-Ptype (オプション引数、文字列):
- 入力ファイル内の物性値データの種類を指定します。データは指定された形式で読み込まれ、その後元の物性値 (P) に変換され、さらに $log_{10}(P)$ に変換されます。
- 選択肢: 'P' (デフォルト), 'log10(P)' (10のP乗に変換), 'log_e(P)' (exp(P)に変換)
- 例: -Ptype="log10(P)"
-xmin, -xmax (オプション引数、浮動小数点数):
- フィッティングに用いるX軸データ ($1000/T$) の最小値と最大値を指定します。この範囲外のデータはフィッティングから除外されます。
- デフォルト値: * (無制限)
- 例: -xmin=2.5 -xmax=3.0
-Tmin, -Tmax (オプション引数、浮動小数点数):
- フィッティングに用いる元の温度データ (ケルビン) の最小値と最大値を指定します。この範囲外のデータはフィッティングから除外されます。
- デフォルト値: * (無制限)
- 例: -Tmin=300 -Tmax=450
-Tcalmin, -Tcalmax (オプション引数、文字列):
- フィッティング結果から物性値を計算しプロットする際の温度範囲を指定します。
- デフォルト値: * (入力データの最小/最大温度)
- 例: -Tcalmin=250 -Tcalmax=500

コマンドラインでの具体的な使用例

ここでは、架空の入力ファイル sample_data.xlsx を使用した具体的な実行例を示します。

sample_data.xlsx の内容:

Temperature	Resistance
273.15	1.0E+06
283.15	8.0E+05
293.15	5.0E+05
303.15	3.0E+05
313.15	1.5E+05
323.15	8.0E+04
333.15	4.0E+04

上記のデータは温度がケルビン (Temperature)、物性値が抵抗 (Resistance) で与えられているとします。

実行コマンド:

python arrhenius_plot.py sample_data.xlsx -model="simple Arrhenius" -Tlabel="Temperature" -Plabel="Resistance" -Ttype="T(K)" -Ptype="P" -Tmin=280 -Tmax=320 -Tcalmin=270 -Tcalmax=340

実行結果の説明:

標準出力 (および sample_data.xlsx に上書きされるログ): プログラムはまず引数と設定値を表示し、入力ファイル sample_data.xlsx を読み込みます。データの変換 (T(K) から $1000/T$、Resistance から $log_{10}(Resistance)$) が行われ、指定された -Tmin=280 から -Tmax=320 の範囲のデータがフィッティングに使用されます。 $1000/T$ と $log_{10}(Resistance)$ に対して一次多項式フィッティングが行われ、係数 ($c_0, c_1$) が算出されます。その後、これらの係数から $P_0$ および $E_a$ の値が表示されます。フィッティングのスコア（元の物性値とフィッティング値、およびそれぞれの対数値）が報告されます。最後に、指定された計算範囲 (-Tcalmin=270, -Tcalmax=340) でプロット用のデータが生成されたことが表示されます。

出力例の一部:

#======================================================
# Analyze activation energy etc by Arrhenius plot
#======================================================
infile : sample_data.xlsx
outfile: sample_data-out.xlsx
logfile: sample_data.xlsx
model : simple Arrhenius
Tlabel: Temperature
Plabel: Resistance
Ttype : T(K)
Ptype : P
Fitting x range: -1e+100 - 1e+100
Fitting T range: 280.0 - 320.0

Read [sample_data.xlsx]
...
T(K)=[273.15, 283.15, 293.15, 303.15, 313.15, 323.15, 333.15]
P   =[1000000.0, 800000.0, 500000.0, 300000.0, 150000.0, 80000.0, 40000.0]

Least-squares fitting with 1-th order polynomial of 1000/T
Data to be fitted:
  1000/T       log10(P)
   3.531         5.903
   3.385         5.778
   3.298         5.602
   3.193         5.176
Coefficients:
  c0=  1.536e+01
  c1= -2.678e+00

  P0=  2.399e+15
  Ea=phi_0=     0.5368 eV

Scores between P(input) and P(fit)
  y1: P(input)
  y2: P(fit)
  R^2        = 0.99220967
  R          = 0.99609722
  RMSE       = 59286.09549320
  max_abs_err= 92839.06250000
  max_rel_err= 0.18567812

Scores between log10(P(input)) and log10(P(fit))
  y1: log10(P(input))
  y2: log10(P(fit))
  R^2        = 0.99849514
  R          = 0.99924729
  RMSE       = 0.04018788
  max_abs_err= 0.06385288
  max_rel_err= 0.01277058

Calculate data from the fitting result
  T range:  270.000 -  340.000 K

Save to [sample_data-out.xlsx]

生成されるExcelファイル (sample_data-out.xlsx): このファイルには、元のデータ、変換後のデータ、フィッティング結果、および計算されたプロット用データが複数の列にわたって保存されます。log10(P)(cal) の列にはフィッティング曲線から計算された値、P(cal) の列にはその逆変換値が格納され、グラフのプロットに使用されます。
表示されるグラフ: 4つのサブプロットを持つMatplotlibウィンドウが表示されます。
- 左上: 元の温度 (K) 対抵抗 (Ω) のデータポイント。
- 右上: 変換後のケルビン温度 (K) 対抵抗 (Ω) のデータポイントと、フィッティングから計算された滑らかな曲線。
- 左下: アレニウスプロット ($1000/T$ 対 $log_{10}(Resistance)$)。黒丸は生データ、赤丸はフィッティングに使用されたデータ、赤線はフィッティング曲線。
- 右下: $1000/T$ に対する活性化エネルギー $E_a$ (eV) のプロット。この例では一次フィッティングのため、$E_a$ は一定値として表示されます。
ユーザーはウィンドウを閉じるまで、これらのグラフを閲覧できます。