xrd_fit.py 技術ドキュメント
プログラムの動作
xrd_fit.py は、X線回折（XRD）実験データを解析し、複数の結晶相の寄与を評価するためのPythonプログラムです。主に、結晶構造情報ファイル（CIF）から理論的なX線回折パターンをシミュレーションし、これを実験データと比較・フィッティングすることで、多相混合物の定量分析や回折ピークの帰属を支援します。
主な機能:
データ読み込み: テキストファイル形式の実験XRDデータおよびCIF形式の結晶構造ファイルを読み込みます。
XRDパターンシミュレーション: CIFファイルから指定されたX線波長、2θ範囲、ピークの半値幅（FWHM）、ガウス-ローレンツ混合比に基づいて理論的なX線回折パターンを計算します。
データ正規化: 観測されたXRDデータとシミュレーションされたCIFデータを0から1の範囲に正規化します。
グラフ表示:
plot: 入力XRDデータのみをプロットします。plot_all: 入力XRDデータと、各CIFファイルから計算された回折パターンおよび回折線を個別のサブプロットに表示します。plot_one: 入力XRDデータと、各CIFファイルから計算された回折線を単一のグラフに重ねて表示します。plot_one2: 入力XRDデータと、各CIFファイルから計算された回折パターンを単一のグラフに重ねて表示します。
相関分析:
overwrap: 各CIF由来のXRDパターン間のオーバーラップ（相関）を視覚的にチェックします。CIFcorrelation: 各CIF由来のXRDパターン間の相関係数を計算し、表示します。correlation: 入力XRDデータと各CIF由来のXRDパターン間の相関係数を計算し、スマアリングのFWHMを変化させた場合の相関の変化をプロットします。
LASSO回帰フィッティング:
fitモードでは、入力XRDデータを、複数のCIF由来XRDパターンとLegendre多項式による背景成分の線形結合でフィッティングします。LASSO（Least Absolute Shrinkage and Selection Operator）回帰を用いることで、不必要な相の寄与をゼロに近づけるL1正則化が適用され、スパースな解（少数の主要な相に絞られた解）を得ることができます。フィッティング結果として、各相および背景成分の係数、MAE、MSE、RMSEなどの評価指標が出力されます。
解決する課題:
多相混合物のXRDパターン解析における、各相の相対寄与の推定。
実験パターンにおける未知のピークの帰属。
背景成分の分離と除去。
各結晶相の存在量の定量的評価。
原理
本プログラムは、X線回折の物理的原理と統計的な機械学習手法を組み合わせてXRDデータの解析を行います。
X線回折パターンの計算:
ブラッグの法則: 結晶格子からのX線回折は、ブラッグの法則 \(n\lambda = 2d \sin\theta\) に従います。ここで、\(n\) は回折次数、\(\lambda\) はX線波長、\(d\) は結晶面間隔、\(2\theta\) は回折角です。
ピーク形状: 回折ピークは、装置的要因や結晶子のサイズ・歪みなどにより広がります。本プログラムでは、
GaussLorentz関数を使用してピーク形状をモデル化します。これは、ガウス関数とローレンツ関数の線形結合または畳み込みとして表現されます。 $\( I(x) = C_G \exp\left(-\frac{(x - x_0)^2}{2\sigma^2}\right) + C_L \frac{1}{1 + \frac{(x - x_0)^2}{\gamma^2}} \)\( ここで、\)I(x)\( は強度、\)x\( は角度、\)x_0\( はピーク中心、\)\sigma\( はガウス成分の幅、\)\gamma\( はローレンツ成分の幅、\)C_G, C_L$ はそれぞれの係数です。Gfraction(ガウス比率) パラメータでガウス成分とローレンツ成分の混合比を調整します。構造因子: CIFファイルから読み取られた結晶構造情報（原子の種類と位置、空間群など）に基づき、各結晶面 (hkl) の構造因子 \(F_{hkl}\) を計算します。回折強度は \(|F_{hkl}|^2\) に比例します。
背景のモデル化:
XRDパターンには、試料ホルダーや非晶質成分、散乱などによる背景（バックグラウンド）が含まれます。これをモデル化するために、Legendre多項式が使用されます。Legendre多項式は、特定の区間で直交性を持つ多項式であり、背景の緩やかな変化を表現するのに適しています。 $\( P_n(x) = \frac{1}{2^n n!} \frac{d^n}{dx^n} (x^2 - 1)^n \)$ プログラムでは、指定された次数（
BG_order）までのLegendre多項式が背景の基底関数として用いられます。
データ正規化と畳み込み:
異なるデータセット間での比較を容易にするため、全てのXRDパターン（実験、CIFシミュレーション、背景成分）は0から1の範囲に正規化されます。
convolution関数により、指定されたFWHMとガウス比率を持つピーク形状（GaussLorentz）を用いて、XRDパターンにスマアリング処理を施すことができます。これにより、実験データにおけるピークのブロードニング効果をシミュレーションパターンに反映させることが可能です。
LASSO回帰フィッティング:
実験的に観測されたXRD強度 \(y_{obs}\) を、各結晶相の理論パターン \(y_{cif,i}\) と背景多項式 \(P_j\) の線形結合でモデル化します。 $\( y_{obs} = \sum_{i=1}^{N_{phases}} c_i y_{cif,i} + \sum_{j=0}^{BG_{order}} b_j P_j + \epsilon \)\( ここで、\)c_i\( は各結晶相の寄与係数、\)b_j\( は背景多項式の係数、\)\epsilon$ は誤差です。
本プログラムでは、LASSO (Least Absolute Shrinkage and Selection Operator) 回帰を使用します。これは、L1正則化項を含む線形回帰手法です。 $\( \min_{\mathbf{w}} \left( \frac{1}{2n_{samples}} ||\mathbf{X}\mathbf{w} - \mathbf{y}_{obs}||_2^2 + \alpha ||\mathbf{w}||_1 \right) \)\( ここで、\)\mathbf{X}\( は各CIFパターンと背景多項式を特徴量とする行列、\)\mathbf{y}_{obs}\( は観測強度ベクトル、\)\mathbf{w}\( は係数ベクトル（\)c_i\( と \)b_j\( を含む）、\)\alpha$ は正則化パラメータです。
L1正則化項 \(\alpha ||\mathbf{w}||_1\) は、係数の絶対値の合計を最小化しようとします。これにより、重要でない係数が厳密にゼロになる傾向があり、特徴選択（この場合、存在しない相を「選択しない」）の効果があります。パラメータ \(\alpha\) を調整することで、フィッティングの厳しさを制御できます。
相関係数:
2つのXRDパターン \(A\) と \(B\) の相関係数は、一般的に以下のように計算されます。 $\( \text{Corr}(A, B) = \frac{\sum (A_i - \bar{A})(B_i - \bar{B})}{\sqrt{\sum (A_i - \bar{A})^2} \sqrt{\sum (B_i - \bar{B})^2}} \)$ ただし、本プログラムでは正規化されたパターンに対して
np.dot(xlist[i], xlist[j]) / norm[i] / norm[j]の形式で計算されており、これは原点を通る線形回帰の相関係数またはコサイン類似度に近いです。パターンが非負で正規化されているため、この値が高いほど類似性が高いことを示します。
必要な非標準ライブラリとインストール方法
xrd_fit.py の実行には、以下のPython非標準ライブラリが必要です。
numpy: 数値計算
scipy: 科学技術計算（特殊関数
legendre、補間interp1dなど）matplotlib: グラフ描画
pandas: データ分析（直接は使用されていないがインポートされている）
scikit-learn (sklearn): 機械学習（LASSO回帰
Lasso、評価指標mean_absolute_error,mean_squared_error）tklib: カスタムライブラリ（
tkApplication,tkFile,tkUtils,tkVariousData,tksci.tksci(Gaussian, Lorentzian, GaussLorentz),tksci.tkmatrix,tkcrystal.tkxrd,tkgraphic.tkplotevent）
tklib 以外のライブラリは、pip コマンドでインストールできます。
pip install numpy scipy matplotlib pandas scikit-learn
tklib は、このプログラムの作者または配布元が提供するカスタムライブラリです。通常、pip で直接インストールすることはできません。以下のいずれかの方法で入手・設定する必要があります。
ソースコードを入手:
tklibライブラリのソースコードを別途入手し、xrd_fit.pyと同じディレクトリ、またはPythonのパスが通っているディレクトリに配置します。開発元から提供されるインストール手順に従う:
tklibの開発元が提供するインストール手順（例えば、GitHubからのクローンとセットアップスクリプトの実行など）に従ってください。
必要な入力ファイル
xrd_fit.py は、以下の2種類の入力ファイルを必要とします。
実験X線回折データファイル
ファイル名: コマンドライン引数
input_pathで指定します。ワイルドカード（例:data/*.txt）も使用可能です。形式: 通常は2列以上の数値データを含むテキストファイル（
.txtなど）。プログラム内部のプラグインメカニズムにより、様々なフォーマットに対応できるよう設計されていますが、基本的には1列目に回折角（2θ）、2列目以降に強度データが記述されている形式が期待されます。例:
2theta Intensity 10.000 150 10.020 155 ... 120.000 200
結晶構造情報ファイル (CIFファイル)
ファイル名: コマンドライン引数
cif_dirで指定します。ディレクトリパスやワイルドカード（例:cif_data/*.cif）も使用可能です。プログラムは指定されたパス内のすべてのCIFファイルを読み込みます。形式: Crystallographic Information File (CIF) 形式。これは結晶構造データを記述するための標準的なテキストベースのフォーマットです。
例: CIFファイルは構造情報、格子定数、原子位置、空間群などの詳細な情報を含みます。
生成される出力ファイル
xrd_fit.py は、以下の出力ファイルを生成します。
ログファイル:
ファイル名: 入力ファイル
input_pathのファイル名に基づいて自動生成されます。例えば、sample.txtが入力ファイルの場合、sample-out.txtという名前で同じディレクトリに保存されます。内容: プログラムの実行中の標準出力（コンソールに表示される内容）が全て記録されます。これには、パラメータ設定、ファイル読み込み状況、相関計算結果、LASSOフィッティングの過程と最終的な係数、エラーメッセージなどが含まれます。
グラフ表示:
各モード (
plot,plot_all,plot_one,plot_one2,overwrap,correlation,fit) の実行に応じて、matplotlibを用いたグラフが画面に表示されます。これらのグラフは、プログラムが終了するまで表示され、ユーザーの入力（"Press ENTER to terminate>>"）を待って閉じられます。
グラフ画像としてファイルに自動保存する機能は含まれていませんが、表示されたグラフウィンドウから手動で画像として保存することが可能です（matplotlibの標準機能）。
コマンドラインでの使用例 (Usage)
xrd_fit.py は、コマンドライン引数で動作モードとパラメータを指定して実行します。
python xrd_fit.py mode input_path cif_dir wavelength Q2min Q2max Q2step fwhm Gfraction fwhm_smear Gfraction_smear yscale BG_order LASSO_alpha
引数説明:
mode: 実行モードを指定します。以下のいずれか。plot: 入力XRDデータのみをプロット。overwrap: CIF XRDパターン間のオーバーラップをチェック。CIFcorrelation: CIF XRDパターン間の相関係数を計算。correlation: 入力XRDデータとCIF XRDパターン間の相関係数を計算。plot_all: 各CIFパターンと入力データを個別のグラフで表示。plot_one: 各CIFパターンと入力データを一つのグラフに回折線として重ねて表示。plot_one2: 各CIFパターンと入力データを一つのグラフにパターンとして重ねて表示。fit: LASSO回帰によるフィッティングを実行。
input_path: 実験XRDデータファイルのパス（ワイルドカード可）。cif_dir: CIFファイルが格納されているディレクトリまたはファイルパス（ワイルドカード可）。wavelength: X線の波長。CuKa,MoKaなどの一般的なX線源を指定。Q2min: 2θの最小値（度）。Q2max: 2θの最大値（度）。Q2step: 2θのステップ幅（度）。fwhm: CIFパターンシミュレーション時のピークの半値幅（FWHM、度）。Gfraction: CIFパターンシミュレーション時のガウス関数の割合 (0.0: Lorentzian, 1.0: Gaussian)。fwhm_smear: 相関解析やフィット前にデータを畳み込む際の追加のスマアリングFWHM（度）。Gfraction_smear: 追加スマアリング時のガウス関数の割合。yscale: グラフのY軸スケール (linearまたはlog)。BG_order: 背景フィッティングに使用するLegendre多項式の次数。LASSO_alpha: LASSO回帰の正則化パラメータ \(\alpha\)。
コマンドラインでの具体的な使用例
ここでは、sample.txt という実験XRDデータファイルと、cif_data/ ディレクトリ内のすべてのCIFファイル (*.cif) を使用する例を示します。
plot_oneモードで入力データとCIF由来の回折線を重ねて表示する例:python xrd_fit.py plot_one sample.txt cif_data/*.cif CuKa 20.0 80.0 0.02 0.2 0.5 0.0 0.0 linear 3 0.0
mode:plot_oneinput_path:sample.txtcif_dir:cif_data/*.cifwavelength:CuKaQ2min:20.0Q2max:80.0Q2step:0.02fwhm:0.2(CIFシミュレーションのFWHM)Gfraction:0.5(CIFシミュレーションのガウス比率)fwhm_smear:0.0(追加スマアリングなし)Gfraction_smear:0.0(追加スマアリングなし)yscale:linearBG_order:3LASSO_alpha:0.0
実行結果の説明:
sample.txtの観測XRDパターンがプロットされ、その上にcif_data/内の各CIFファイルから計算された理論的な回折ピーク位置と強度が回折線として重ねて表示されます。これにより、どの結晶相が観測ピークに寄与しているかを視覚的に確認できます。Y軸は線形スケールです。fitモードでLASSO回帰によるフィッティングを実行する例:python xrd_fit.py fit sample.txt cif_data/*.cif CuKa 20.0 80.0 0.02 0.2 0.5 0.1 0.5 linear 3 0.01
mode:fitinput_path:sample.txtcif_dir:cif_data/*.cifwavelength:CuKaQ2min:20.0Q2max:80.0Q2step:0.02fwhm:0.2(CIFシミュレーションのFWHM)Gfraction:0.5(CIFシミュレーションのガウス比率)fwhm_smear:0.1(追加スマアリングFWHM)Gfraction_smear:0.5(追加スマアリングのガウス比率)yscale:linearBG_order:3(3次Legendre多項式で背景をフィッティング)LASSO_alpha:0.01(LASSO回帰の正則化パラメータ)
実行結果の説明: プログラムは、まずLASSO回帰の正則化パラメータ \(\alpha\) を変化させた場合のRMSEと各係数の変化をプロットしたグラフ（
fig_lasso）を表示します。これにより、適切な \(\alpha\) の選択を支援します。 次に、指定された \(\alpha=0.01\) でLASSO回帰が実行されます。標準出力には、背景多項式の各次数および各CIFファイルに対応する相の寄与係数が表示されます。また、MAE、MSE、RMSEといったフィッティングの評価指標も出力されます。 最後に、観測XRDデータ、フィッティング結果、フィッティングされた背景、および各結晶相のフィッティング済みパターンを重ねて表示したグラフが表示されます。これにより、フィッティングの品質と各成分の寄与度を視覚的に評価できます。ログファイルsample-out.txtにもこれらの情報が詳細に記録されます。