xrd_fit_xrayutilities.py 技術ドキュメント
プログラムの動作
xrd_fit_xrayutilities.py は、X線回折 (XRD) シミュレーション、フィッティング、および膜厚推定のためのPythonスクリプトです。xrayutilities ライブラリの動的理論に基づき、薄膜サンプルのXRDパターンを解析します。
主な機能は以下の通りです。
XRDシミュレーション (
--mode sim): 指定された構造パラメータ（組成、緩和度、膜厚）に基づいて、理論的なXRDパターンを計算しプロットします。入力ファイルが指定された場合、実験データとの比較も可能です。データ読み込み (
--mode read): テキストファイルから2theta-intensityデータを読み込み、そのパターンをプロットします。構造パラメータの最適化 (
--mode fit): 実験的なXRDパターンにシミュレーションパターンをフィッティングさせ、組成、緩和度、膜厚などの構造パラメータを最適化します。以下の最適化手法をサポートしています。ランダム探索
SciPyの最適化アルゴリズム (Nelder-Mead, BFGS, CG)
粒子群最適化 (PSO) 最適化の進行状況はリアルタイムでプロットに表示されます。
膜厚の自動推定 (
--mode guess): 実験データからフリンジパターンを解析し、高速フーリエ変換 (FFT) と線形回帰を用いて膜厚を自動的に推定します。推定された膜厚候補はCSVファイルに出力され、フィッティングの初期値として利用できます。
このプログラムは、薄膜の非破壊的な構造解析において、実験データと理論モデルを比較し、物理的なパラメータを決定するという課題を解決します。
原理
XRDシミュレーションの原理
本プログラムのXRDシミュレーションは、xrayutilities ライブラリの動的理論に基づいて行われます。動的理論は、X線が結晶中で多重散乱される効果を考慮に入れるため、特に薄膜やエピタキシャル層のような高分解能なXRDパターンを正確に記述できます。
材料モデル:
xu.materials.Crystal.fromCIFを用いて、基板と膜の結晶構造情報 (CIFファイル) からCrystalオブジェクトを構築します。膜は通常、基板と異なる組成を持つ合金としてxu.materials.material.Alloyでモデル化されます。層構造と歪み:
xu.simpack.Layerオブジェクトで各層（基板、膜）の厚さや緩和度を定義します。膜層の緩和度 (relax) は、膜が基板との格子不整合をどの程度緩和しているかを示します (0: 完全擬似格子、1: 完全緩和)。これらの層はxu.simpack.PseudomorphicStack001を使用して積層構造を形成します。このスタックは、膜と基板間の格子整合や歪みを考慮したモデルです。動的理論モデル: 最終的に
xu.simpack.DynamicalModelオブジェクトが、上記で定義されたスタック、X線エネルギー、回折指数 (HKL, 例: (0,0,2))、および分解能幅を用いて動的理論シミュレーションを実行するためのモデルを構築します。X線エネルギーはCu Kα1波長から計算されます（Energy [eV] = 12390.0 / Wavelength [Å]）。強度計算:
dyn_model.simulate()メソッドにより、指定されたオメガ角 (omega) 範囲での回折強度分布が計算されます。
最適化（フィッティング）の原理
フィッティングは、シミュレーションで得られたXRDパターンと実験データとの間の差（残差）を最小化する目的関数に基づいています。
目的関数: シミュレーション強度 \(I_{\text{sim}}\) と実験強度 \(I_{\text{target}}\) の残差は、L2ノルム（二乗和）として計算されます。
--residual_scaleオプションに応じて、線形スケールまたは対数スケールで計算されます。対数スケールは、強度の低いフリンジ部分のフィッティング精度を向上させるのに役立ちます。 $\( \text{Residual} = \sum_{i} \left( \log\left(\max(I_{\text{target},i}, \epsilon)\right) - \log\left(\max(I_{\text{sim},i}, \epsilon)\right) \right)^2 \quad \text{(if residual\_scale = "log")} \)\( \)\( \text{Residual} = \sum_{i} \left( I_{\text{target},i} - I_{\text{sim},i} \right)^2 \quad \text{(if residual\_scale = "linear")} \)\( ここで、\)\epsilon\( は数値安定性のための微小な正の値（\)10^{-30}$）。ペナルティ関数: フィッティングパラメータ（組成 \(x\)、緩和度 \(relax\)、膜厚 \(thick\)）が物理的に妥当な範囲（\(0 \le x \le 1\)、\(0 \le relax \le 1\)、\(thick > 0\)）を逸脱した場合、目的関数に大きなペナルティが加算され、最適化がこれらの制約内で解を見つけるように誘導されます。
最適化アルゴリズム:
ランダム探索: パラメータをランダムに微小変化させ、目的関数が改善すればその変化を採用するという、単純な局所探索手法です。
SciPyの最適化アルゴリズム:
scipy.optimize.minimizeを利用し、Nelder-Mead (勾配フリー)、BFGS (準ニュートン法)、CG (共役勾配法) などの標準的な手法を使用します。粒子群最適化 (PSO): 複数の「粒子」が解空間を探索する群知能アルゴリズムです。各粒子は自身の最良解 (
pbest) と群全体の最良解 (gbest) に基づいて速度と位置を更新します。これにより、局所最適解に陥りにくいグローバルな探索が可能です。
膜厚推定（フリンジ解析）の原理
--mode guess では、XRDフリンジパターンから膜厚を自動的に推定します。
信号抽出: まず、指定された角度範囲のデータが切り出されます。強度は
--yscaleに基づき対数または線形スケールに変換されます。次に、scipy.signal.savgol_filter(Savitzky-Golayフィルター) を用いて、なだらかなエンベロープ（バックグラウンドやブラッグピークの形状）が抽出され、元の信号からエンベロープを差し引くことでフリンジ成分のみの残差信号が分離されます。FFTによる初期推定: 抽出された残差信号はQ空間に変換され、
numpy.fft.rfft(高速フーリエ変換) にかけられます。フリンジはQ空間で周期的な変動として現れるため、FFTスペクトル中の主要な周波数成分はフリンジのQ空間における周期 \(\Delta q\) に対応します。膜厚 \(T\) は以下の関係から推定されます。 $\( T = \frac{2\pi}{\Delta q} \)\( FFTから得られた最も支配的な \)\Delta q$ が、フリンジ次数割り当ての基準として使用されます。ピーク検出とクラスター化:
scipy.signal.find_peaksを使用して、残差信号からフリンジピークが検出されます。検出されたピークは、隣接するピーク間のQ値の差が、FFTで推定された \(\Delta q\) にcluster_gap_factorを掛けた値よりも大きい場合に、個別の「クラスター」として分割されます。これにより、フリンジが途切れている場合でも個別に解析できます。線形回帰による精密膜厚推定: 各クラスター内で、検出されたフリンジピークのQ値 \(Q_i\) に、FFTで推定された \(\Delta q\) を基準として整数フリンジ次数 \(n_i\) が割り当てられます。フリンジのQ値と次数の間には近似的に線形関係があります。 $\( Q_n \approx Q_0 + n \cdot \Delta q \)\( この関係に対し2段階の線形回帰が適用されます。まず全てのピークで回帰し、その後、残差が標準偏差の2倍を超える外れ値を除外して再回帰を行います。最終的な回帰直線の傾き \)\frac{\Delta Q}{\Delta n}\( から、より精密な膜厚 \)T = \frac{2\pi}{|\text{傾き}|}$ が算出されます。
必要な非標準ライブラリとインストール方法
このプログラムの実行には、以下の非標準Pythonライブラリが必要です。
numpy: 数値計算matplotlib: グラフ描画xrayutilities: X線回折シミュレーション (動的理論)scipy: 最適化アルゴリズム、信号処理 (Savitzky-Golayフィルター、ピーク検出)
これらのライブラリは、pip コマンドを使用してインストールできます。
pip install numpy matplotlib xrayutilities scipy
必要な入力ファイル
プログラムの実行モードによって、以下の入力ファイルが必要となります。
XRDデータファイル (
--infile)形式: テキストファイル
内容: 2列の数値データ。1列目は2theta角度 [deg]、2列目はX線強度 [arb.u.] です。ヘッダー行は不要です。
例 (
my_xrd_data.txt):32.000 120.5 32.005 130.1 32.010 145.8 ...
備考:
--mode simで--infileが指定されない場合、内部的に模擬データが生成されます。
基板のCIFファイル (
--substrate_file)形式: CIF (Crystallographic Information File)
内容: 基板材料の結晶構造情報（格子定数、原子位置、空間群など）。
xrayutilitiesで読み込み可能な形式である必要があります。デフォルト:
GaN.cif
膜のCIFファイル (
--film_file)形式: CIF (Crystallographic Information File)
内容: 膜材料の結晶構造情報。
デフォルト:
AlN.cif
パラメータファイル (自動読み込み/生成)
形式: CSVファイル
内容: 最適化される構造パラメータ（x, relax, thick）の現在の値と、最適化対象であるかどうかのフラグ (
optid) を含みます。ファイル名:
--infileがpath/to/data.txtの場合、path/to/data-parameters.csvとなります。--infileが指定されない場合はparameters.csvです。備考: このファイルはプログラムが自動的に生成・更新します。ユーザーが手動で作成することも可能ですが、ヘッダー行 (
varname,value,optid) と各パラメータの形式（例:best0:x,0.05,1）に注意が必要です。
生成される出力ファイル
プログラムは、実行モードと引数に応じて以下のファイルを出力します。
パラメータファイル:
ファイル名:
--infileがpath/to/data.txtの場合、path/to/data-parameters.csvとなります。--infileが指定されない場合はparameters.csvです。内容: 最適化された、または推定された最新の構造パラメータ (
x,relax,thick) の値と、それらが最適化対象 (optid=1) か固定 (optid=0) かを示すフラグをCSV形式で保存します。PSOモードで実行された場合、複数のPSO粒子候補のパラメータも保存されます（例:pso0:x,0.048,1）。例:
varname,value,optid best0:x,0.0512345678901234,1 best0:relax,0.0123456789012345,1 best0:thick,1512.345678901234,1 pso0:x,0.0501,1 pso0:relax,0.0115,1 pso0:thick,1510.89,1 ...
膜厚推定モード (
--mode guess) 関連ファイル:フリンジピーク情報ファイル:
ファイル名:
--infileがpath/to/data.txtの場合、path/to/data-guess-peaks.csvとなります。--infileが指定されない場合はguess-peaks.csvです。内容: 膜厚推定モードで検出された個々のフリンジピークの詳細情報（所属クラスターID、フリンジ次数
n、Q値、2theta角、線形回帰で採用されたかどうかのフラグ、各回帰ステップでの残差など）をCSV形式で保存します。
フリンジクラスター情報ファイル:
ファイル名:
--infileがpath/to/data.txtの場合、path/to/data-guess-clusters.csvとなります。--infileが指定されない場合はguess-clusters.csvです。内容: 膜厚推定モードで解析された各フリンジクラスターのサマリー（クラスターID、ピーク数、FFTによる推定Δq、FFTと線形回帰による推定膜厚、信頼度スコアなど）をCSV形式で保存します。
プロット:
各モード（
sim,read,guess,fit）の実行中に、matplotlibを用いたグラフウィンドウがGUIで表示されます。これらのプロットは、デフォルトではファイルに保存されません。
コマンドラインでの使用例 (Usage)
プログラムの基本的な使い方と引数は、--help オプションで確認できます。
python xrd_fit_xrayutilities.py --help
出力例:
usage: xrd_fit_xrayutilities.py [-h] [--mode {sim,read,fit,guess}] [--method {random,nelder-mead,bfgs,cg,pso}] [--infile INFILE] [--substrate_file SUBSTRATE_FILE] [--film_file FILM_FILE] [--fix FIX] [--nmaxiter NMAXITER] [--tol TOL] [--yscale {linear,log}] [--residual_scale {linear,log}] [--guess_low GUESS_LOW] [--guess_high GUESS_HIGH] [--nsmooth_points NSMOOTH_POINTS] [--nsmooth_order NSMOOTH_ORDER] [--nguess_keep NGUESS_KEEP] [--cluster_gap_factor CLUSTER_GAP_FACTOR] [--pso_nparticles PSO_NPARTICLES] [--pso_w PSO_W] [--pso_c1 PSO_C1] [--pso_c2 PSO_C2] [--pso_stall_max PSO_STALL_MAX] [--pso_spread_rtol PSO_SPREAD_RTOL]
動的理論に基づくXRDシミュレーション、フィッティング、および膜厚推定
options:
-h, --help show this help message and exit
--mode {sim,read,fit,guess}
実行モード: sim (シミュレーション), read (読込), fit (最適化), guess (膜厚推定) (デフォルト: sim)
--method {random,nelder-mead,bfgs,cg,pso}
最適化手法。fitモード時に有効 (デフォルト: pso)
--infile INFILE 2theta-intensityデータを含むテキストファイルへのパス (デフォルト: )
--substrate_file SUBSTRATE_FILE
基板のCIFファイルパス (デフォルト: GaN.cif)
--film_file FILM_FILE
膜のCIFファイルパス (デフォルト: AlN.cif)
--fix FIX 固定するパラメータ（カンマ区切り、例: "x,relax"） (デフォルト: )
--nmaxiter NMAXITER 最大反復回数 (デフォルト: 1000)
--tol TOL 収束判定の許容誤差 (デフォルト: 1e-07)
--yscale {linear,log}
プロットのY軸スケール (デフォルト: "log")
--residual_scale {linear,log}
残差計算に使用するスケール (デフォルト: "log")
--guess_low GUESS_LOW
guessモードでの解析開始角度 [deg] (デフォルト: 16.0)
--guess_high GUESS_HIGH
guessモードでの解析終了角度 [deg] (デフォルト: 17.3)
--nsmooth_points NSMOOTH_POINTS
guessモードでの平滑化ウィンドウ点数 (デフォルト: 31)
--nsmooth_order NSMOOTH_ORDER
guessモードでの平滑化多項式次数 (デフォルト: 3)
--nguess_keep NGUESS_KEEP
guessモードで保持する膜厚候補の数 (デフォルト: 5)
--cluster_gap_factor CLUSTER_GAP_FACTOR
フリンジクラスター判定用のギャップ係数 (デフォルト: 1.6)
--pso_nparticles PSO_NPARTICLES
PSOの粒子数 (デフォルト: 12)
--pso_w PSO_W PSOの慣性重み (デフォルト: 0.72)
--pso_c1 PSO_C1 PSOの自己学習係数 (デフォルト: 1.49)
--pso_c2 PSO_C2 PSOの社会学習係数 (デフォルト: 1.49)
--pso_stall_max PSO_STALL_MAX
改善がない場合にPSOを停止する最大反復数 (デフォルト: 15)
--pso_spread_rtol PSO_SPREAD_RTOL
粒子の分散に基づく停止条件の相対許容誤差 (デフォルト: 0.1)
コマンドラインでの具体的な使用例
以下の例では、デフォルトのCIFファイル (GaN.cif, AlN.cif) を使用していると仮定します。
XRDシミュレーションの実行 指定されたパラメータ（デフォルト値を使用、または
parameters.csvから読み込み）でXRDパターンをシミュレートし、プロットします。python xrd_fit_xrayutilities.py --mode sim
実行結果: GaN基板上のAlN膜のXRDパターン（組成
x=0.0、緩和度relax=1.0、膜厚1500.0 Å）がシミュレートされ、ログスケールの強度でプロット表示されます。
既存のXRDデータファイルの読み込みと表示
my_xrd_data.txtというファイルからXRDデータを読み込み、プロットします。python xrd_fit_xrayutilities.py --mode read --infile my_xrd_data.txt --yscale linear
実行結果:
my_xrd_data.txtの2theta-intensityデータが読み込まれ、線形スケールの強度でプロット表示されます。
フリンジ解析による膜厚の自動推定
my_xrd_data.txtの指定された角度範囲のフリンジから膜厚を推定し、結果をCSVファイルに保存します。python xrd_fit_xrayutilities.py --mode guess --infile my_xrd_data.txt --guess_low 34.0 --guess_high 36.0 --nguess_keep 3
実行結果:
my_xrd_data.txtの34.0度から36.0度の範囲のデータがフリンジ解析されます。Savitzky-Golayフィルターによるエンベロープ除去、FFT解析、フリンジピーク検出、線形回帰が実行され、推定された膜厚候補がコンソールに表示されます。また、my_xrd_data-guess-peaks.csvとmy_xrd_data-guess-clusters.csvが生成され、最も信頼性の高い膜厚がmy_xrd_data-parameters.csvに保存されます。分析の各ステップ（信号、残差、FFTスペクトル、n-qプロット）を示すグラフが複数表示されます。
Nelder-Mead法を用いた組成固定でのフィッティング
my_xrd_data.txtに対して、組成 (x) を固定し、緩和度 (relax) と膜厚 (thick) をNelder-Mead法で最適化します。python xrd_fit_xrayutilities.py --mode fit --method nelder-mead --infile my_xrd_data.txt --fix "x" --nmaxiter 500
実行結果: フィッティングの進行状況がライブプロットでリアルタイムに表示され、
nmaxiterで指定された最大500回反復します。最適化後、最終的なパラメータがmy_xrd_data-parameters.csvに保存されます。コンソールには、各反復での残差とパラメータ値が出力されます。
粒子群最適化 (PSO) を用いたフィッティング
my_xrd_data.txtに対して、組成、緩和度、膜厚の全てをPSOで最適化します。粒子数を20、最大反復回数を200に設定します。python xrd_fit_xrayutilities.py --mode fit --method pso --infile my_xrd_data.txt --pso_nparticles 20 --nmaxiter 200
実行結果: PSOアルゴリズムにより、複数の粒子が解空間を探索しながら最適なパラメータを探索します。フィッティングの進行状況はライブプロットで表示され、コンソールには群全体の最良解 (
gbest) の残差とパラメータ値が反復ごとに表示されます。最終的なパラメータはmy_xrd_data-parameters.csvに保存され、PSO粒子群の最良なものもあわせて保存されます。