xrd_fit.py 技術ドキュメント
プログラムの動作
xrd_fit.py は、X線回折 (XRD) シミュレーション、実験データのフィッティング、およびフリンジ解析による膜厚自動推定を行うためのPythonスクリプトです。このプログラムは、薄膜試料のXRDパターンから、材料の組成比 (\(x\))、緩和度 (relax)、膜厚 (thick) といった物理的構造パラメータを非破壊で推定することを目的としています。
主な機能は以下の通りです。
XRDパターンシミュレーション (
simモード): xrayutilitiesライブラリの動的理論に基づき、指定された構造パラメータ（組成、緩和度、膜厚）を用いてXRDパターンを計算し表示します。XRDデータ読み込み表示 (
readモード): テキストファイル形式の2θ-Intensityデータを読み込み、プロットして表示します。構造パラメータ最適化フィッティング (
fitモード): 実験データに対してシミュレーションパターンをフィッティングし、構造パラメータ（組成、緩和度、膜厚）を最適化します。複数の最適化アルゴリズム（ランダム探索、Nelder-Mead、BFGS、CG、粒子群最適化(PSO)）をサポートします。膜厚自動推定 (
guessモード): 実験データに現れるフリンジ（振動構造）を解析し、高速フーリエ変換(FFT)や線形回帰を用いて膜厚の候補値を自動で推定します。推定結果はパラメータファイルに保存され、フィッティングの初期値として利用できます。リアルタイムプロット: フィッティングや膜厚推定の進捗状況をMatplotlibのグラフでリアルタイムに表示します。
原理
XRDシミュレーションの原理
xrd_fit.py は、xrayutilities ライブラリを用いて、動的理論に基づいたXRDシミュレーションを実行します。
材料定義: 基板と膜それぞれの結晶構造情報 (CIFファイル) を読み込み、
xrayutilities.materials.Crystalオブジェクトを生成します。六方晶系の結晶を前提としており、弾性定数が未定義の場合はデフォルト値が適用されます。合金層: 膜は基板と膜材料の合金として扱われ、組成比 \(x\) が導入されます。
層スタック: 基板層と膜層を組み合わせて
xrayutilities.simpack.PseudomorphicStack001を作成します。膜は基板に対して擬似整合 (pseudomorphic) な成長をすると仮定し、緩和度 (relaxation) パラメータを考慮して格子定数を調整します。モデル構築: 上記のスタック情報、X線エネルギー（デフォルトはCu Kα1線）、回折ピークの広がり（
resolution_width）を設定し、xrayutilities.simpack.DynamicalModelを構築します。強度計算: 構築されたモデルに対し、指定された角度範囲と回折指数 (デフォルト: HKL = (0,0,2)) でシミュレーションを実行し、強度分布を計算します。
最適化フィッティングの原理
フィッティングは、シミュレーションで得られたXRDパターンと実験データとの間の残差（目的関数）を最小化することで行われます。
目的関数: 以下の残差 \(R\) を最小化します。 $\(R = \sum_{i} (Y_{exp, i} - Y_{sim, i})^2 + P\)\( ここで \)Y_{exp, i}\( は実験強度、\)Y_{sim, i}\( はシミュレーション強度です。 `--residual_scale` オプションで `log` を選択した場合、強度は対数変換されてから残差が計算されます。 \)\(R = \sum_{i} (\log(I_{exp, i} + \epsilon) - \log(I_{sim, i} + \epsilon))^2 + P\)\( \)\epsilon\( は数値計算上の安定性のための小さな値です。 \)P\( は物理的に妥当な範囲 (\)x \in [0,1]\(, \)relax \in [0,1]\(, \)thick > 0$) からパラメータが外れた場合に大きな値を返すペナルティ項です。
最適化アルゴリズム:
scipy.optimize.minimizeを利用した汎用的な最適化アルゴリズムと、粒子群最適化 (PSO) のカスタム実装、およびシンプルなランダム探索をサポートします。ランダム探索 (
random): パラメータをランダムに微小量変化させ、残差が改善した場合にそのパラメータを採択するシンプルな探索手法です。Nelder-Mead, BFGS, CG (
nelder-mead,bfgs,cg):scipy.optimize.minimizeが提供する局所最適化アルゴリズムです。初期値の与え方や探索効率に違いがあります。粒子群最適化 (PSO) (
pso): 複数の粒子（候補解）が解空間を探索し、自身の過去の最良解 (pbest) と群全体の最良解 (gbest) に基づいて移動するグローバル最適化手法です。 各粒子 \(p\) の位置 \(\mathbf{x}_p\) と速度 \(\mathbf{v}_p\) は、以下の式で更新されます。 $\(\mathbf{v}_p(t+1) = w \mathbf{v}_p(t) + c_1 r_1 (\mathbf{pbest}_p - \mathbf{x}_p(t)) + c_2 r_2 (\mathbf{gbest} - \mathbf{x}_p(t))\)\( \)\(\mathbf{x}_p(t+1) = \mathbf{x}_p(t) + \mathbf{v}_p(t+1)\)\( ここで、\)w\( は慣性重み、\)c_1\( は自己学習係数、\)c_2\( は社会学習係数、\)r_1, r_2\( は \)[0,1]$ の一様乱数です。 PSOは、一定の反復回数で改善が見られない場合（pso_stall_max）や、粒子の位置が十分に収束した場合（pso_spread_rtol）に早期停止します。
膜厚推定の原理 (guess モード)
膜厚推定モードでは、実験データに現れるフリンジ（LaueフリンジやPendellösungフリンジ）の周期性を解析して膜厚を求めます。
信号前処理:
指定された角度範囲 (
--guess_low,--guess_high) のデータを抽出。強度を対数変換（またはそのまま使用）。
Savitzky-Golayフィルター (
scipy.signal.savgol_filter) で滑らかなエンベロープを抽出し、元の信号からエンベロープを差し引くことでフリンジ成分（残差）を抽出します。
高速フーリエ変換 (FFT) による解析:
角度2θを波数 \(q_z\) に変換します。 $\(q_z = \frac{4\pi}{\lambda} \sin\theta\)\( ここで \)\lambda\( はX線波長、\)\theta$ は半角 (2θ の半分) です。
フリンジ成分に対してFFTを適用し、スペクトルのピークから主要なフリンジ周期 \(\Delta q_z\) を特定します。
膜厚 \(T\) は以下の関係から推定されます。 $\(T = \frac{2\pi}{\Delta q_z}\)$
フリンジピーク検出と線形回帰:
残差信号からフリンジピークを検出 (
scipy.signal.find_peaks)。ピーク間の距離を基準にフリンジが欠損している箇所でピーク群をクラスターに分割 (
cluster_gap_factor)。各クラスター内で、連続するフリンジピークに整数次数の
nを割り当てます。nと \(q_z\) の関係が線形であるという仮定に基づき、線形回帰 (numpy.polyfit) を行います。 $\(q_z = A \cdot n + B\)$この回帰の傾き \(A = \frac{\Delta q_z}{\Delta n}\) を用いて、膜厚 \(T\) を推定します。 $\(T = \frac{2\pi}{A}\)$
よりロバストな推定のため、外れ値を除外する2段階の線形回帰が適用されます。
必要な非標準ライブラリとインストール方法
xrd_fit.py の実行には、以下の非標準Pythonライブラリが必要です。
numpy: 数値計算matplotlib: グラフ描画xrayutilities: XRDシミュレーションscipy: 最適化、信号処理
これらのライブラリは、以下の pip コマンドでインストールできます。
pip install numpy matplotlib xrayutilities scipy
必要な入力ファイル
プログラムの動作モードによって、以下のファイルが必要となります。
XRD実験データファイル (
--infile)形式: 2列のテキストファイル (スペースまたはタブ区切り)。
内容:
1列目: 2θ角度 [deg]
2列目: X線回折強度 [arb. u.]
例 (
my_exp_data.txt):32.000 123.45 32.005 125.67 32.010 128.90 ...
--mode simで--infileを指定しない場合、プログラム内部でターゲットパラメータから模擬データが生成されます。
基板の結晶構造ファイル (
--substrate_file)形式: Crystallographic Information File (CIF)。
内容: 基板材料の結晶構造情報。
デフォルト:
GaN.cif
膜の結晶構造ファイル (
--film_file)形式: Crystallographic Information File (CIF)。
内容: 膜材料の結晶構造情報。
デフォルト:
AlN.cif
パラメータ設定ファイル (自動生成/読み込み)
形式: CSV (Comma Separated Values) ファイル。
内容: 最適化する構造パラメータ (
x,relax,thick) の値と、最適化対象とするかどうかのフラグ (optid)。ファイル名:
--infileが指定された場合は[infile_basename]-parameters.csv、そうでない場合はparameters.csv。このファイルは、プログラムが初回実行時にデフォルト値で自動生成するか、既存のファイルを読み込みます。手動で編集して初期値を設定することも可能です。
例 (
my_exp_data-parameters.csv):varname,value,optid best0:x,0.2000000000000000,1 best0:relax,0.0000000000000000,1 best0:thick,1500.0000000000000000,1 pso0:x,0.1987654321000000,1 pso0:relax,0.0012345678900000,1 pso0:thick,1498.7654321000000000,1
varnameは[セット名]:[パラメータ名]の形式で指定されます。デフォルトの最適化結果はbest0セットに保存されます。PSOモードではpso0,pso1などの名前で各粒子の初期位置が保存されることがあります。
生成される出力ファイル
プログラムは実行モードに応じて以下のファイルを生成します。
パラメータ設定ファイル (
[infile_basename]-parameters.csvまたはparameters.csv)内容: フィッティングの最適化結果（組成 \(x\)、緩和度
relax、膜厚thickの最終値）と、各パラメータが最適化対象であったかどうかのフラグが保存されます。PSOモードでは、複数の粒子（候補解）のパラメータもこのファイルに保存されます。形式: 「必要な入力ファイル」セクションで示されたCSV形式と同様です。
膜厚推定ピークデータファイル (
[infile_basename]-guess-peaks.csvまたはguess-peaks.csv)モード:
guessモード実行時に生成されます。内容: フリンジ解析で検出された各ピークのID、フリンジ次数 \(n\)、波数 \(q\)、2θ角度、および2段階線形回帰における使用状況と残差情報が保存されます。
形式: CSVファイル。
例:
cluster_id,n,q,2theta,used_stage1,used_stage2,resid_stage1,resid_stage2 0,0,2.23456789,33.12345,1,1,0.000123,-0.000012 0,1,2.24567890,33.23456,1,1,0.000056,0.000005 ...
膜厚推定クラスターデータファイル (
[infile_basename]-guess-clusters.csvまたはguess-clusters.csv)モード:
guessモード実行時に生成されます。内容: フリンジ解析で検出された各クラスター（フリンジ群）ごとの情報（ピーク数、FFTによる推定膜厚、線形回帰による推定膜厚、信頼度など）が保存されます。
形式: CSVファイル。
例:
cluster_id,npeaks,dq_est,thick_fft,thick_reg,confidence 0,10,0.01234567,510.123,510.567,10.0 1,5,0.01234567,510.123,511.123,5.0
リアルタイムプロットウィンドウ:
fitおよびguessモードでは、Matplotlibによるグラフウィンドウがリアルタイムで更新され、フィッティングの進捗やフリンジ解析の結果（信号、エンベロープ、残差、FFTスペクトル、n-qプロットなど）が表示されます。これらのウィンドウはファイルとして自動保存されませんが、ユーザーが手動で保存できます。
コマンドラインでの使用例 (Usage)
xrd_fit.py --help を実行すると、以下のヘルプメッセージが表示されます。
usage: xrd_fit.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)
コマンドラインでの具体的な使用例
以下の例では、my_exp_data.txt という実験データファイルが存在し、GaN.cif と AlN.cif がプログラムと同じディレクトリにあると仮定します。
シミュレーションモードで模擬データを生成し、プロットする
--infileを指定しない場合、プログラム内のデフォルトパラメータ (TARGET_X,TARGET_RELAX,TARGET_THICK) を使用してシミュレーションが実行されます。
python xrd_fit.py --mode sim
実行結果: Matplotlibウィンドウにシミュレーション結果のXRDパターンが表示されます。
読み込みモードで実験データをプロットする
python xrd_fit.py --mode read --infile my_exp_data.txt --yscale log
実行結果: Matplotlibウィンドウに
my_exp_data.txtのXRDパターンがY軸対数スケールで表示されます。
膜厚推定モードでフリンジ解析を実行する
python xrd_fit.py --mode guess --infile my_exp_data.txt --guess_low 34.0 --guess_high 36.0 --nguess_keep 3
実行結果:
ターミナルに推定された膜厚候補がリスト形式で表示されます。
Matplotlibウィンドウにフリンジ解析の各ステップ（信号、エンベロープ、残差、FFTスペクトル、n-qプロット）が表示されます。
my_exp_data-guess-peaks.csvとmy_exp_data-guess-clusters.csvが出力されます。my_exp_data-parameters.csvが生成または更新され、最も信頼性の高い膜厚候補がbest0:thickとして保存されます。
フィッティングモードでPSO最適化を実行する
--infileで指定されたデータに対して、パラメータx,relax,thickを最適化します。
python xrd_fit.py --mode fit --method pso --infile my_exp_data.txt --nmaxiter 500 --pso_nparticles 20
実行結果:
ターミナルに反復回数、現在の残差、各パラメータ値が表示されます。
Matplotlibウィンドウに実験データと現在のシミュレーションパターンのフィッティング状況がリアルタイムで表示されます。
フィッティング完了後、
my_exp_data-parameters.csvに最適化されたパラメータの最終値が保存されます。
フィッティングモードでNelder-Mead最適化を実行し、特定のパラメータを固定する
この例では、組成
xを初期値に固定し、relaxとthickのみを最適化します。
python xrd_fit.py --mode fit --method nelder-mead --infile my_exp_data.txt --fix x --nmaxiter 1000
実行結果:
fix xオプションにより、my_exp_data-parameters.csv内のxのoptidが0に設定され、最適化中に \(x\) の値は変化しません。それ以外の動作は上記のフィッティング例と同様です。