lsq_latt2.py Technical Documentation
プログラムの動作
lsq_latt2.py は、X線回折や電子回折などの回折データから、最小二乗法を用いて結晶の格子定数を高精度に決定するためのPythonプログラムです。主に、以下の目的と機能を提供し、結晶学における構造解析の基礎となる正確な格子定数決定という課題を解決します。
目的と主な機能
格子定数計算: 回折ピークの位置（回折角など）とミラー指数（hkl）から、結晶の逆格子および直接格子定数（\(a, b, c, \alpha, \beta, \gamma\)）とその標準偏差を算出します。
結晶系への対応: 三斜晶系、単斜晶系、斜方晶系、正方晶系、六方晶系、菱面体晶系（Hexagonal setting）、菱面体晶系（Rhombohedral setting）など、複数の結晶系に対応したモデルで計算を行います。
回折位置の変換と重み付け: 入力された生の位置データ（例えば2\(\theta\)）を最小二乗法に適した形式に変換し、測定値の信頼性に応じた重みを適用します。
系統誤差の補正: いくつかの一般的な系統誤差（例: 吸収、アライメント、プロファイル形状に関連する誤差）を補正する機能を提供します。
外れ値の検出と除去: 統計的に有意な外れ値を自動的に検出し、計算から除外して結果のロバスト性を向上させます。
詳細なレポート出力: 算出された格子定数、その標準偏差、個々の観測値に対する残差などを含む詳細なレポートを生成します。
解決する課題
このプログラムは、実験的に得られた回折データに内在する測定誤差や系統誤差を考慮し、統計的に最も確からしい格子定数を導き出すことで、結晶構造解析における精密な格子パラメータの決定を可能にします。特に、複数の回折ピークデータから一貫性のある格子定数を決定する際に有効です。
原理
lsq_latt2.py は、回折データに基づいて格子定数を決定するために、重み付き最小二乗法と誤差伝播の原理を適用します。
1. 最小二乗法の基本モデル
回折データにおける格子定数決定の基本的な関係は、ブラッグの法則と結晶学的な式に基づいています。逆格子空間におけるベクトル \(d^*_{hkl}\) の長さは、実空間における面間隔 \(d_{hkl}\) の逆数に比例します。 $\(d^*_{hkl} = \frac{1}{d_{hkl}} = \frac{2 \sin\theta}{\lambda}\)\( ここで、\)\lambda\( はX線（または電子線）の波長、\)\theta\( はブラッグ角です。 したがって、\)1/d^2_{hkl}\( は \)\sin^2\theta\( に比例します。この関係を利用して、一般に最小二乗法のモデルは次のように構築されます。 \)\(Y = \frac{4}{\lambda^2} \sin^2\theta = \frac{1}{d^2_{hkl}}\)\( この \)Y\( を従属変数とし、ミラー指数 \)(h, k, l)\( の関数として表現される格子定数パラメータに対する線形モデルを立てます。 一般に、三斜晶系の逆格子ベクトル \)d^_{hkl}\( の二乗は次のように表されます。 \)\( \frac{1}{d^2_{hkl}} = a^{*2}h^2 + b^{*2}k^2 + c^{*2}l^2 + 2b^*c^*kl\cos\alpha^* + 2c^*a^*lh\cos\beta^* + 2a^*b^*hk\cos\gamma^* \)\( ここで、\)a^, b^, c^\( は逆格子定数、\)a^{*2}, b^{*2}, c^{2}, 2b^c^\cos\alpha^, \dots\( は最小二乗法で決定されるパラメータとなります。 このプログラムでは、\)Y\( の値を `transformed_position` 関数で計算し、右辺の式を `lattice_design_row` 関数で生成される設計行列 \)X\( の行として表現します。 つまり、 \)\( Y_i = \mathbf{X}_i \boldsymbol{\beta} + \epsilon_i \)\( ここで、\)Y_i\( は \)i\( 番目の観測に対する変換された位置、 \) \mathbf{X}_i \( は \)i\( 番目の観測に対する設計行列の行、\)\boldsymbol{\beta}\( は格子定数に関連するパラメータのベクトル、\)\epsilon_i$ は残差です。
2. 重み付き最小二乗法
測定の不確かさ（誤差）は観測値によって異なるため、重み \(w_i\) を導入して重み付き最小二乗法が適用されます。重みは通常、測定値の分散の逆数として与えられます。
目的関数は以下のようになります。
$\( S(\boldsymbol{\beta}) = \sum_{i=1}^{N} w_i (Y_i - \mathbf{X}_i \boldsymbol{\beta})^2 = (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{W} (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta}) \)\(
ここで、\)\mathbf{Y}\( は従属変数ベクトル、\)\mathbf{X}\( は設計行列、\)\mathbf{W}\( は対角に重み \)w_i\( を持つ重み行列です。
この目的関数を最小化するパラメータ \)\boldsymbol{\beta}\( の推定値 \)\hat{\boldsymbol{\beta}}\( は、正規方程式を解くことで得られます。
\)\( \hat{\boldsymbol{\beta}} = (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W} \mathbf{Y} \)$
この計算は weighted_least_squares 関数で行われます。
3. 系統誤差の補正
correction_terms 関数では、様々な実験的な系統誤差を考慮するための補正項が設計行列に追加されます。これらの補正項は、例えば角度依存の吸収補正や試料の位置ずれによる誤差をモデル化します。
4. 外れ値の除去
計算された残差が標準偏差の OUTLIER_SIGMA（デフォルトは3.0）倍を超える場合、その観測値は外れ値として識別され、最小二乗計算から除外されます。その後、残りのデータで計算が繰り返され、結果のロバスト性が向上します。
5. 誤差伝播
最小二乗法によって決定されたパラメータ \(\hat{\boldsymbol{\beta}}\) とその共分散行列から、最終的な直接格子定数と逆格子定数、およびそれらの標準偏差が導出されます。このプログラムでは、SigmaPropagator クラスを用いて誤差伝播を計算します。
これは、関数 \(f(x_1, x_2, \dots, x_n)\) の標準偏差 \(\sigma_f\) を、各入力変数 \(x_i\) の標準偏差 \(\sigma_{x_i}\) から近似的に求める手法です。具体的には、各入力変数 \(x_i\) を微小量 \(\epsilon_i = \text{PERTURB_SCALE} \cdot \sigma_{x_i}\) だけずらしたときの \(f\) の変化量を用いて数値微分を近似し、以下の式に基づいて分散を計算します。
$\( \sigma_f^2 \approx \sum_{i=1}^n \left( \frac{f(x_1, \dots, x_i+\epsilon_i, \dots, x_n) - f(x_1, \dots, x_n)}{\epsilon_i} \right)^2 \sigma_{x_i}^2 \)$
この手法は、パラメータ間の相関を直接考慮しない簡略化された誤差伝播です。
必要な非標準ライブラリとインストール方法
lsq_latt2.py は、以下の非標準ライブラリを利用しています。
numpy: 数値計算、特に行列演算とベクトル演算に広く使用されます。最小二乗法の正規方程式の解法や共分散行列の計算に不可欠です。
インストール方法
numpy は pip を使って簡単にインストールできます。コマンドプロンプトやターミナルで以下のコマンドを実行してください。
pip install numpy
必要な入力ファイル
入力ファイルはプレーンテキスト形式で、一つ以上のジョブデータを含めることができます。各ジョブは以下の構造で記述されます。
<ジョブタイトル>
<結晶系フラグ> <補正フラグ1> <補正フラグ2> <補正フラグ3> <補正フラグ4> <補正フラグ5> <重みモード> <位置モード> <プロファイルモード>
<初期波長> <初期半径>
<h1 k1 l1 生の回折位置1>
[0点シフト1 (位置モード3の場合)]
[シグマまたは重み1 (重みモード0, 1, 3, 4の場合)]
<h2 k2 l2 生の回折位置2>
[0点シフト2 (位置モード3の場合)]
[シグマまたは重み2 (重みモード0, 1, 3, 4の場合)]
...
<hN kN lN 生の回折位置N>
[0点シフトN (位置モード3の場合)]
[シグマまたは重みN (重みモード0, 1, 3, 4の場合)]
<反射データ終了コード (1000)>
<次のジョブタイトル> (または最終ジョブ終了コード 0)
<次のジョブのパラメータ...>
...
<最終ジョブ終了コード (0)>
ジョブタイトル: そのジョブを一意に識別する文字列。
結晶系フラグ (lattice_system): 結晶系を指定する整数 (1-8)。
1: 三斜晶系 (Triclinic)
2: 単斜晶系 (Monoclinic, \(\alpha^* = \gamma^* = 90^\circ\))
3: 単斜晶系 (Monoclinic, \(\beta^* = \gamma^* = 90^\circ\))
4: 斜方晶系 (Orthorhombic)
5: 正方晶系 (Tetragonal)
6: 立方晶系 (Cubic)
7: 菱面体晶系 (Trigonal/Rhombohedral, Hexagonal setting)
8: 六方晶系 (Hexagonal)
補正フラグ (correction_flags): 5つの整数 (
0または1)。特定の系統誤差補正を有効にするかどうか。重みモード (weight_mode): 重み計算方法を指定する整数 (0-4)。
0: \(\sigma(\theta)\) を使用
1: \(\sigma(2\theta)\) を使用
2: 一定の重み (1.0)
3: 入力された重みを使用
4: 入力された \(\sigma\) を使用し、重みは \(1/\sigma^2\)
位置モード (position_mode): 生の回折位置データを変換する方法を指定する整数 (1-6)。
1: \(\sin(\theta)\)
2: 生の回折位置をそのまま使用
3: 0点シフトを考慮した \(\sin(\theta_{true})\)
4: \(\sin(\theta/2)\)
5: \(1/d\)
6: \(x/\sqrt{1+x^2}\)
プロファイルモード (profile_mode): 現在のコードでは使用されていません。
初期波長 (wavelength): 回折測定に使用された波長。
初期半径 (radius): 回折計の半径など、幾何学的な定数。
h k l 生の回折位置: ミラー指数と生の回折位置（通常は2\(\theta\)値など）。
hが2000の場合、その行は波長と半径の更新を指定します。hが1000の場合、そのジョブの反射データ入力の終了を示します。
0点シフト (zero_shift):
position_modeが3の場合のみ、h k l raw_position行の直後に来る単一の数値。シグマまたは重み (sigma_or_weight):
weight_modeが0, 1, 3, 4の場合のみ、h k l raw_position行の直後に来る単一の数値。
入力ファイル例 (sample_input.txt):
Example Job - Triclinic System
1 0 0 0 0 0 2 1 1
1.540598 100.0
1 0 0 30.0
0 1 0 35.0
0 0 1 40.0
1 1 0 45.0
1 0 1 50.0
0 1 1 55.0
1 1 1 60.0
1000
0
この例では、三斜晶系 (1) で補正なし (0)、重みは一定 (2)、位置は \(\sin(\theta)\) 変換 (1) を使用しています。波長は1.540598、半径は100.0で、7つの反射データが与えられています。
生成される出力ファイル
出力ファイルは、コマンドラインで指定されたパスにテキスト形式で生成されます。各ジョブの計算結果が詳細なレポート形式で記述されます。
ファイル名
コマンドライン引数 output_file で指定されます（例: output.txt）。
内容
出力ファイルには、以下の情報が含まれます。
プログラムヘッダ: プログラム名とバージョン情報。
ジョブタイトル: 入力ファイルから読み取られたジョブのタイトル。
設定パラメータ: 結晶系、補正フラグ、重みモード、位置モード、プロファイルモードの数値。
使用波長: 計算に使用されたX線波長。
最小二乗パラメータ: 最小二乗法によって決定されたモデルパラメータ \(P(I)\) と、それらの標準偏差 \(SIGP(I)\)。
観測値の詳細:
観測番号 (No.)
ミラー指数 (H K L)
観測された面間隔 \(D_{obs}\) (Dobs.)
計算された面間隔 \(D_{cal}\) (Dcal.)
生の回折位置 (BO)
計算された2\(\theta\)値 (BC)
残差 (Resid.)
生の回折位置と計算された2\(\theta\)値の差 (\(d(2Q)\))
外れ値の報告: 外れ値として計算から除外された観測値があれば、その詳細が示されます。
逆格子定数 (Reciprocal cell constant):
逆格子軸の長さ (\(a^*, b^*, c^*\)) とその標準偏差。
逆格子角 (\(\alpha^*, \beta^*, \gamma^*\)) とその標準偏差。
逆格子角のコサイン値 (\(\cos\alpha^*, \cos\beta^*, \cos\gamma^*\)) とその標準偏差。
逆格子体積 (\(V^*\)) とその標準偏差。
直接格子定数 (Direct cell constant):
直接格子軸の長さ (\(a, b, c\)) とその標準偏差。
直接格子角 (\(\alpha, \beta, \gamma\)) とその標準偏差。
直接格子角のコサイン値 (\(\cos\alpha, \cos\beta, \cos\gamma\)) とその標準偏差。
直接格子体積 (\(V\)) とその標準偏差。
フッター: 計算終了を示すメッセージ。
コマンドラインでの使用例 (Usage)
lsq_latt2.py は、コマンドライン引数を使用して実行します。
基本的な実行コマンド
python lsq_latt2.py <input_file> <output_file> [options]
<input_file>: 必要な入力ファイル のセクションで説明されている形式の入力データファイルへのパスを指定します。<output_file>: 生成される出力ファイル のセクションで説明されている計算結果レポートを保存するファイルへのパスを指定します。
オプション
--silent: プログラムの起動時および終了時に標準出力に表示されるバナーメッセージを抑制します。
コマンドラインでの具体的な使用例
前述の「必要な入力ファイル」セクションで示した sample_input.txt を使用して、プログラムを実行する例を示します。
sample_input.txt の内容
Example Job - Triclinic System
1 0 0 0 0 0 2 1 1
1.540598 100.0
1 0 0 30.0
0 1 0 35.0
0 0 1 40.0
1 1 0 45.0
1 0 1 50.0
0 1 1 55.0
1 1 1 60.0
1000
0
実行コマンド
sample_input.txt を入力として、report.txt という名前の出力ファイルを生成します。
python lsq_latt2.py sample_input.txt report.txt
実行時の標準出力
--silent オプションを付けない場合、以下のようなメッセージが標準出力に表示されます。
Least Square Calculation for Lattice Constant
Refactored Python version
Input : sample_input.txt
Output: report.txt
Lattice Constant Calculation was finished!
report.txt の実行結果（抜粋）
生成される report.txt ファイルの内容は、以下のような形式になります。具体的な数値は、入力データと計算精度によって異なります。
Least Square Calculation for Lattice Constant
Refactored Python version
Example Job - Triclinic System
HLS LE(1) LE(2) LE(3) LE(4) LE(5) IW IO IP
1 0 0 0 0 0 2 1 1
Wave Length =1.540598
P(I) SIGP(I)
0.0270501 0.0000030
0.0356598 0.0000037
0.0463970 0.0000042
0.0000000 0.0000000
0.0000000 0.0000000
0.0000000 0.0000000
H K L Dobs. Dcal. BO BC Resid. d(2Q)
( 1) 1 0 0 2.5694 2.5695 30.000 30.001 -0.00000 -0.00140
( 2) 0 1 0 2.2963 2.2963 35.000 35.000 0.00000 0.00018
( 3) 0 0 1 2.0620 2.0619 40.000 40.001 0.00000 -0.00062
( 4) 1 1 0 1.8797 1.8798 45.000 44.999 0.00000 0.00064
( 5) 1 0 1 1.7317 1.7317 50.000 50.000 0.00000 -0.00012
( 6) 0 1 1 1.6111 1.6111 55.000 55.000 0.00000 -0.00004
( 7) 1 1 1 1.5126 1.5126 60.000 60.000 0.00000 0.00005
Reciprocal cell constant
a* b* c*
0.1644685( 0.0000091) 0.1888383( 0.0000098) 0.2154009( 0.0000098)
alpha* beta* gamma*
90.000000( 0.000000) 90.000000( 0.000000) 90.000000( 0.000000)
cos a* cos b* cos c*
0.000000( 0.000000) 0.000000( 0.000000) 0.000000( 0.000000)
V*=0.0067069116(0.0000000000)
Direct cell constant
a b c
6.08053( 0.00022) 5.29562( 0.00019) 4.64239( 0.00016)
alpha beta gamma
90.0000( 0.0000) 90.0000( 0.0000) 90.0000( 0.0000)
cos a cos b cos c
0.0000( 0.0000) 0.0000( 0.0000) 0.0000( 0.0000)
V=163.51347895(0.00000000)
Lattice Constant Calculation was finished!
出力結果の説明:
このレポートには、各回の観測データに対する計算結果と残差、そして最終的に決定された逆格子と直接格子の定数が、それぞれの標準偏差とともに記載されます。特に、Resid. 列は最小二乗法モデルからのずれを示し、d(2Q) 列は観測された生の位置と計算された2\(\theta\)値の差を示します。格子定数は、P(I) で示される最小二乗パラメータから derive_cell_constants 関数によって導出され、その不確かさも誤差伝播によって計算されています。
この例では、入力がシンプルであったため、残差は非常に小さく、単斜晶系を仮定しているためか、角度は全て90度として計算されています。