光電子分光スペクトル解析ドキュメント

プログラムの動作

ef_fit.py は、光電子分光 (PES) スペクトルデータを解析するためのPythonプログラムです。このプログラムは、フェルミ・ディラック (Fermi-Dirac) 関数と線形状態密度 (DOS) を組み合わせた物理モデルを、ガウス型の装置関数で畳み込み、測定されたスペクトルにフィッティングすることで、モデルのパラメータを決定します。

主な機能と解決する課題は以下の通りです。

PESスペクトル解析: 測定されたPESスペクトルに対し、物理モデル（フェルミ関数、DOS、装置関数）に基づいたフィッティングを行います。
パラメータ決定: フィッティングにより、フェルミ準位 (EF)、バックグラウンド (BG)、装置関数幅 (w)、信号強度 (I0)、DOSの線形傾き (DOSa) といった重要な物理パラメータを最適化して決定します。
柔軟なモード:
- fit モード: 測定データにモデルをフィッティングし、最適なパラメータを探索します。最適化の進行状況はリアルタイムでプロットされ、終了後には結果がExcelファイルに保存されます。
- sim モード: 指定されたパラメータでモデルスペクトルを計算し、測定データと比較して表示します。これにより、パラメータ変更がスペクトルに与える影響を視覚的に確認できます。
- init モード: ユーザーが指定したフィット範囲に基づいて、EF、BG、I0、w の初期値を自動で推定し、結果をパラメータCSVファイルに保存します。これにより、フィッティングの開始点を見つける手間が省けます。
パラメータ管理: パラメータはCSVファイル ([入力ファイル名]_parameters.csv) で管理され、プログラム実行時に読み込まれます。コマンドライン引数で上書きすることも可能です。これにより、実験条件に応じたパラメータの再利用や微調整が容易になります。
多様な入力ファイル形式: .csv, .xlsx, .txt (特定のHAXPESデータ形式またはシンプルなXY形式) のスペクトルデータを読み込むことができます。
リアルタイムプロット: フィッティングプロセス中のモデルスペクトルの変化をリアルタイムで視覚化し、最適化の収束状況を把握できます。

このプログラムは、PESスペクトルの定量的な解析を支援し、材料の電子状態に関する物理的洞察を得ることを目的としています。

原理

ef_fit.py は、光電子分光スペクトルが以下の物理モデルによって記述されるという前提に基づいています。

\[I_{model}(E) = I_0 \cdot (D(E, DOS_a) \cdot f(E, EF, T)) * G(E, k_G) + BG\]

ここで、各項は以下の通りです。

フェルミ・ディラック関数 ($f(E, EF, T)$): 電子が特定のエネルギー準位 $E$ を占有する確率を表します。 $$f(E, EF, T) = \frac{1}{\exp\left(\frac{E - EF}{k_B T}\right) + 1}$$
- $E$: エネルギー
- $EF$: フェルミ準位（電子占有確率が0.5となるエネルギー）
- $T$: 絶対温度
- $k_B$: ボルツマン定数 ($8.617333 \times 10^{-5}$ eV/K)
状態密度 (DOS, $D(E, DOS_a)$): このプログラムでは、単純な線形関数としてモデル化されています。 $$D(E, DOS_a) = 1.0 - DOS_a \cdot E$$
- $DOS_a$: 状態密度の線形傾きパラメータ。
装置関数 ($G(E, k_G)$): 実験装置のエネルギー分解能に起因するスペクトルの広がりを表します。ガウス関数でモデル化されています。 $$G(dE, k_G) = \exp(-k_G \cdot dE^2)$$
- $dE$: エネルギー差
- $k_G$: ガウス関数の幅を決定するパラメータ。
- $k_G$ と半値全幅 (FWHM) $w$ の関係は以下の通りです。 $$k_G = \frac{\ln(2)}{(0.5w)^2}$$ または、逆変換として $$w = 2.0 \sqrt{\frac{\ln(2)}{k_G}}$$ プログラム内部では、w (FWHM) をユーザーが指定し、それを kG に変換して計算に用います。畳み込み (*) は、numpy.convolve を用いて効率的に実装されています。
$I_0$: ピーク信号の強度をスケーリングする係数。
$BG$: バックグラウンド強度。

フィッティングアルゴリズム: fit モードでは、scipy.optimize.minimize 関数を用いて、測定データ $I$ とモデルスペクトル $I_{model}$ の間の残差二乗和 $S$ を最小化します。

\[S = \sum_{i} (I_i - I_{model,i})^2\]

最適化は、Nelder-Mead法 (デフォルト) を含む複数のアルゴリズムで実行可能です。

パラメータの変換と制約:

w (FWHM) と I0 (信号強度) は常に正の値を取るべき物理量であるため、最適化時にはこれらのパラメータの対数 (log(w), log(I0)) を用いて探索が行われます。これにより、パラメータが負になることを防ぎ、広い範囲での探索を可能にします。最適化後、結果は指数関数で物理量に戻されます。
各パラメータには、物理的に妥当な範囲 (min, max) が設定されており、最適化中にこの範囲を逸脱した場合には、目的関数に大きなペナルティ項が加算されます。これにより、探索が物理的に意味のある領域に留まるように誘導されます。

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

ef_fit.py の実行には、以下のPython非標準ライブラリが必要です。

chardet: ファイルのエンコーディングを検出するために使用されます。
numpy: 数値計算、特に配列操作や数学関数に利用されます。
pandas: データ構造とデータ解析ツールを提供し、CSV/Excelファイルの読み書きに使用されます。
matplotlib: データのプロットと可視化に使用されます。
scipy: 科学技術計算ライブラリで、特にここではscipy.optimize.minimizeによる最適化処理に利用されます。

これらのライブラリは、Pythonのパッケージマネージャ pip を使ってインストールできます。コマンドプロンプトやターミナルで以下のコマンドを実行してください。

pip install chardet numpy pandas matplotlib scipy

必要な入力ファイル

ef_fit.py は、以下の種類の入力ファイルに対応しています。

スペクトルデータファイル: 測定された光電子分光スペクトルデータを含むファイル。
- ファイル形式: .csv, .xlsx, .txt
- データ構造:
  - .csv または .xlsx ファイルの場合: 1列目にエネルギー (eV)、2列目に測定強度が含まれていることを想定しています。ヘッダー行は自動的にスキップされます。
  - .txt ファイルの場合: プログラムは以下の3種類の形式を自動で判別しようとします。
    1. Spring-8 HAXPES 形式: ファイルの1行目に [Info] タグがあり、Dimension 1 size=、Dimension 2 size= などのメタデータ、そして [Data 1] タグの後にデータが続く形式。1列目がエネルギー、それ以降が測定強度（複数の強度データがある場合は平均化されます）。
    2. シンプルな2列XY形式 (Nデータあり): ファイルの2行目にデータ点数を示す整数があり、その後にエネルギーと強度がスペース区切りで並ぶ形式。
    3. シンプルな2列XY形式 (Nデータなし): エネルギーと強度がスペース区切りで並ぶ形式。
パラメータCSVファイル (オプション): フィッティングパラメータの初期値や制約設定を記述したCSVファイル。
- ファイル名: [スペクトルデータファイル名ベース]_parameters.csv
  - 例: spectrum.txt の場合、spectrum_parameters.csv
- データ構造: Parameter, optid, convert, eps, min, max, kpenalty, Value の列を持つDataFrame形式。
  - Parameter: パラメータ名 (例: EF, BG, w, I0, DOSa)
  - optid: 最適化対象とするか (1: 対象, 0: 固定)
  - convert: 最適化時の変換方法 (例: log は対数変換、空欄はそのまま)
  - eps: 収束判定に使用する許容誤差
  - min, max: パラメータの物理的下限/上限
  - kpenalty: 範囲外ペナルティの係数
  - Value: パラメータの現在の値（初期値、または前回のフィッティング結果）
- 動作: このファイルが存在しない場合、プログラムは内部のデフォルト値でパラメータDataFrameを初期化します。存在する場合はその内容を読み込み、コマンドライン引数で指定された値があれば、それが最優先で適用されます。

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

ef_fit.py の実行により、以下のファイルと出力が生成されます。

パラメータCSVファイル:
- ファイル名: [入力ファイル名ベース]_parameters.csv
  - 例: spectrum.txt を入力した場合、spectrum_parameters.csv が生成または更新されます。
- 内容: fit モードまたは init モードでプログラムを実行すると、最適化された（または推定された）パラメータの Value がこのファイルに保存されます。これにより、次回の実行時にこれらの値を初期値として利用できます。
フィッティング結果Excelファイル (fitモード):
- ファイル名: [入力ファイル名ベース]_fit.xlsx
  - 例: spectrum.txt を入力した場合、spectrum_fit.xlsx が生成されます。
- 内容:
  - Energy (eV): スペクトルデータのエネルギースケール。
  - Measured: 測定されたスペクトルの強度。
  - Initial: フィッティング開始前の初期モデルスペクトルの強度。
  - Final: フィッティング後の最終モデルスペクトルの強度。
- このファイルは、測定データ、初期モデル、最終モデルを比較するために利用できます。
リアルタイムプロットウィンドウ:
- 内容:
  - fit モード: 測定データに加え、フィッティングの進行に伴い更新されるモデルスペクトルがリアルタイムで表示されます。これにより、最適化の収束状況を視覚的に確認できます。フィッティング終了後、最終モデルスペクトルが表示されます。
  - sim モード: 測定データ、指定されたパラメータで計算されたモデルスペクトル、さらにフェルミ関数、DOS関数、それらの積、装置関数が個別にプロットされます。
  - init モード: 測定データと、自動推定された初期パラメータに基づくモデルスペクトルが表示されます。推定されたフェルミ準位、バックグラウンド、最大強度が垂直線・水平線で示されます。
- プロットウィンドウはプログラム終了まで開いたままで、終了後には自動で閉じられます。

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

ef_fit.py は、以下の形式でコマンドラインから実行されます。

python ef_fit.py <filename> [--mode <mode>] [--method <method>] [--T <temp>] [--EF <ef_val>] [--BG <bg_val>] [--w <w_val>] [--I0 <i0_val>] [--DOSa <dosa_val>] [--Efit_min <min_val>] [--Efit_max <max_val>] [--Escale <scale_val>] [--nmaxiter <iterations>] [--print_interval <interval>]

引数の説明:

<filename> (必須): 入力スペクトルデータファイルのパス。(csv, xlsx, txtに対応)
--mode {fit, sim, init} (デフォルト: fit): 実行モードを指定します。
- fit: スペクトルデータをフィッティングします。
- sim: 指定されたパラメータでモデルスペクトルをシミュレーションします。
- init: フィッティングパラメータの初期値を推定します。
--method <method> (デフォルト: Nelder-Mead): 最適化アルゴリズムを指定します。scipy.optimize.minimize でサポートされるアルゴリズムが利用可能です（例: BFGS, L-BFGS-B, Powell など）。simplex は nelder-mead に内部で変換されます。
--T <temp> (型: float, デフォルト: 300.0): 温度 [K]。フェルミ・ディラック関数に影響します。
--EF <ef_val> (型: float, デフォルト: 0.0): フェルミ準位 [eV]。
--BG <bg_val> (型: float, デフォルト: 0.0): バックグラウンド強度。
--w <w_val> (型: float, デフォルト: 0.20): 装置関数の半値全幅 (FWHM) [eV]。
--I0 <i0_val> (型: float, デフォルト: 200000.0): 信号強度。
--DOSa <dosa_val> (型: float, デフォルト: 0.0): DOSの線形傾き。
- --EF, --BG, --w, --I0, --DOSa の各引数は、_parameters.csv ファイルに記載された初期値を上書きします。
--Efit_min <min_val> (型: str または float, デフォルト: None): フィッティングまたはシミュレーション範囲の下限エネルギー [eV]。
--Efit_max <max_val> (型: str または float, デフォルト: None): フィッティングまたはシミュレーション範囲の上限エネルギー [eV]。
--Escale <scale_val> (型: float, デフォルト: -1.0): 入力エネルギーに対するスケーリング係数。例えば、負の値を指定することでエネルギースケールを反転させることができます。
--nmaxiter <iterations> (型: int, デフォルト: 1000): 最適化の最大反復回数。
--print_interval <interval> (型: int, デフォルト: 5): 最適化ステップの途中経過をコンソールに表示する間隔。

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

ここでは、架空のスペクトルデータファイル sample_spectrum.txt が存在すると仮定して、各モードの使用例を示します。

1. 初期値推定モード (mode=init)

このモードでは、フィッティングの開始点となるパラメータを自動で推定します。これにより、手動で初期値を設定する手間を省けます。

python ef_fit.py sample_spectrum.txt --mode init --Efit_min -0.5 --Efit_max 0.5 --T 300

説明:

sample_spectrum.txt を入力ファイルとして使用します。
--mode init で初期値推定モードを指定します。
--Efit_min -0.5 --Efit_max 0.5 で、初期値推定を行うエネルギー範囲を -0.5 eV から 0.5 eV に設定します。
--T 300 で温度を 300 K に設定します。

実行結果: コンソールには、推定された EF, BG, w, I0, DOSa の値が表示されます。また、sample_spectrum_parameters.csv というファイルが生成（または更新）され、推定された初期値が保存されます。リアルタイムプロットウィンドウが開き、測定データと、推定された初期値に基づくモデルスペクトル、および推定されたフェルミ準位やバックグラウンド位置が視覚的に表示されます。

Status: Initialized default parameters.
# ... (ファイルの存在状況などによってメッセージは変わる可能性があります)

--- Guess initial parameters ---
Efit_min=-0.5, Efit_max=0.5, T=300K

Fitting:
Escale  : -1.0
Efit_min: -0.5
Efit_max: 0.5

Guessed parameters:
  EF=0.005312
  BG=150.234567
  w=0.038000
  I0=180000.000000
  DOSa=0.000000
  w =0.038000 (kG=1211.751680)

Guessed parameters saved to [sample_spectrum_parameters.csv]

（プロットウィンドウに推定結果が表示されます）

2. フィッティングモード (mode=fit)

このモードでは、スペクトルデータに対してモデルをフィッティングし、最適なパラメータを探索します。

python ef_fit.py sample_spectrum.txt --mode fit --Efit_min -0.5 --Efit_max 0.5 --T 300 --nmaxiter 500 --print_interval 10

説明:

sample_spectrum.txt を入力ファイルとして使用します。
--mode fit でフィッティングモードを指定します。
--Efit_min -0.5 --Efit_max 0.5 で、フィッティングを行うエネルギー範囲を -0.5 eV から 0.5 eV に設定します。
--T 300 で温度を 300 K に設定します。
--nmaxiter 500 で最適化の最大反復回数を 500 に設定します。
--print_interval 10 で、10ステップごとに途中経過をコンソールに出力します。

実行結果: コンソールには、初期パラメータと、最適化の各ステップでのパラメータ値および誤差が表示されます。リアルタイムプロットウィンドウが開き、測定データと、最適化の進行に伴って更新されるモデルスペクトルがアニメーションで表示されます。フィッティング終了後、最終的な最適化パラメータがコンソールに出力され、sample_spectrum_parameters.csv が更新されます。また、sample_spectrum_fit.xlsx というExcelファイルが生成され、測定データ、初期モデル、最終モデルのデータが保存されます。

Fitting:
Escale  : -1.0
Efit_min: -0.5
Efit_max: 0.5
method  : Nelder-Mead
nmaxiter: 500

Initial parameters:
  EF=0.005312
  BG=150.234567
  w=0.038000
  I0=180000.000000
  DOSa=0.000000
  w =0.038000 (kG=1211.751680)

[Step    1] EF=0.00531, BG=150.23457, w=0.03800, I0=180000.00000, DOSa=0.00000, Error=2.500000e+09
[Step   10] EF=0.00498, BG=145.12345, w=0.03750, I0=181000.00000, DOSa=0.00000, Error=2.100000e+09
# ... (最適化の途中経過がprint_intervalごとに表示されます)
[Step  150] EF=0.00123, BG=120.98765, w=0.03012, I0=195000.00000, DOSa=-0.01000, Error=5.000000e+08

Optimization terminated successfully: All fitted parameters changed by less than their respective EPS value (Step 152).
Optimized parameters:
  EF=0.00123
  BG=120.98765
  w=0.03012
  I0=195000.00000
  DOSa=-0.01000
  w =0.03012 (kG=1910.876543)

Status: Parameters saved to [sample_spectrum_parameters.csv].
Fitting results saved to [sample_spectrum_fit.xlsx]

（プロットウィンドウに最終的なフィッティング結果が表示されます）

3. シミュレーションモード (mode=sim)

このモードでは、指定されたパラメータでモデルスペクトルを計算し、測定データと比較して表示します。フィッティングは行われません。

python ef_fit.py sample_spectrum.txt --mode sim --EF 0.001 --BG 120 --w 0.03 --I0 195000 --DOSa -0.01 --T 300 --Efit_min -0.5 --Efit_max 0.5

説明:

sample_spectrum.txt を入力ファイルとして使用します。
--mode sim でシミュレーションモードを指定します。
--EF, --BG, --w, --I0, --DOSa で、シミュレーションに使用するパラメータの値を直接指定します。
--T 300 で温度を 300 K に設定します。
--Efit_min -0.5 --Efit_max 0.5 で、シミュレーション範囲を設定します。

実行結果: コンソールには、シミュレーションに使用されたパラメータが表示されます。リアルタイムプロットウィンドウが開き、測定データ、指定されたパラメータに基づくモデルスペクトル、そして各物理モデル要素（DOS、フェルミ関数、畳み込み前の物理スペクトル、装置関数）が個別にプロットされます。これにより、各パラメータがスペクトル形状にどのように寄与しているかを視覚的に理解できます。

Fitting:
Escale  : -1.0
Efit_min: -0.5
Efit_max: 0.5

Simulation parameters:
使用パラメータ
  EF=0.001000
  BG=120.000000
  w=0.030000
  I0=195000.000000
  DOSa=-0.010000
  w =0.030000 (kG=1925.322896)

（プロットウィンドウにシミュレーション結果が表示されます）