lsq_latt2.py 技術ドキュメント

プログラムの動作

lsq_latt2.py は、X線回折データなどの回折データから結晶の格子定数を最小二乗法を用いて計算するためのPythonスクリプトです。このプログラムは、様々な結晶格子系に対応しており、観測データの適切な前処理、重み付け、そして最小二乗フィッティングを通じて、最終的な逆格子定数および実格子定数を導出します。また、計算の精度を向上させるために、外れ値の自動除去機能も組み込まれています。

主な機能は以下の通りです。

多種類の格子系への対応: 三斜晶、単斜晶、直方晶、正方晶、六方晶、立方晶、菱面体晶など、主要な結晶格子系における格子定数計算が可能です。
柔軟なデータ処理: 入力データに応じて、回折位置の変換（例: $2\theta$ から $ \sin^2\theta $ へ）や、観測データへの重み付け方法（例: 一定重み、測定誤差に基づく重み）を選択できます。
補正項の適用: 測定における系統誤差（例: 試料位置ずれ、偏光因子）を考慮するための複数の補正項を最小二乗モデルに組み込むことができます。
重み付き最小二乗フィッティング: 観測データに信頼度に応じた重みを適用し、格子パラメータの最適な値を決定します。
外れ値の自動除去: 統計的な基準に基づいて、測定誤差が大きいと考えられる観測データを自動的に識別し、計算から除外することで、フィットのロバスト性を高めます。
格子定数の導出と誤差伝播: フィットされたパラメータから、逆格子および実格子の長さ、角度、体積、そしてそれぞれの標準偏差を計算します。標準偏差の計算には誤差伝播の法則が用いられます。
詳細なレポートとプロットの生成: 計算結果は整形されたテキストレポートファイルとして出力され、各ジョブの回折データ $Q$ (散乱ベクトル) 対 $1/d$ を示すプロットが画像ファイルとして生成されます。
複数ジョブの一括処理: 一つの入力ファイル内に複数の計算ジョブを記述することで、連続して処理を行うことができます。

このスクリプトは、結晶構造解析における格子定数精密化の強力なツールとして機能します。

原理

lsq_latt2.py は、X線回折などの回折実験から得られるデータを用いて、結晶の格子定数を最小二乗法により精密化します。その中心的な原理は以下の通りです。

最小二乗法による格子パラメータの決定

回折データは、ブラッグの法則 $2d \sin\theta = n\lambda$ に従い、格子面間隔 $d$ と回折角 $ \theta $ 、波長 $ \lambda $ の関係を示します。ここでは、シンプル化のため次数 $n=1$ とします。結晶格子における格子面間隔 $d$ は、ミラー指数 $ (hkl) $ と逆格子定数 $ (a^*, b^*, c^*, \alpha^*, \beta^*, \gamma^*) $ の関係によって表現されます。一般に、三斜晶系の逆格子定数とミラー指数 $ (hkl) $ の関係は以下の式で与えられます。

\[ \frac{1}{d^2} = h^2 a^{*2} + k^2 b^{*2} + l^2 c^{*2} + 2kl b^*c^* \cos\alpha^* + 2lh c^*a^* \cos\beta^* + 2hk a^*b^* \cos\gamma^* \]

この式を、未知のパラメータ $ \beta $ （逆格子パラメータの二乗や関連項）と既知の観測値 $y$ （通常は $1/d^2$ )、デザイン行列 $X$ （ミラー指数 $h, k, l$ に基づく項）を使った線形モデルとして表現します。

\[ y = X \beta + \epsilon \]

ここで、$ \epsilon $ は誤差項です。プログラムでは、観測された回折位置 $raw\_position$ から、選択された $position\_mode$ に応じて目的変数 $y$ が計算されます。例えば、position_mode=1（sin(theta)）が選択された場合、$y$ は $ (4 / \lambda^2) \times \sin^2\theta $ として計算され、これは $1/d^2$ に比例します。

重み付き最小二乗法: 各観測データの信頼度を考慮するため、重み付き最小二乗法が適用されます。各観測 $i$ に対して重み $w_i$ が与えられ、重み行列 $W = \text{diag}(w_1, w_2, \dots, w_n)$ が構築されます。最適なパラメータ $ \hat{\beta} $ は、以下の正規方程式を解くことで得られます。

\[ (X^T W X) \hat{\beta} = X^T W y \]

ここで $X^T$ はデザイン行列 $X$ の転置行列です。

デザイン行列の構築: 格子系に応じてデザイン行列 $X$ の構造が異なります。例えば、立方晶系では $1/d^2 \propto (h^2+k^2+l^2)a^{*2}$ となり、パラメータは $a^{*2}$ のみです。各格子系には固有の結晶学的制約があり、これらがデザイン行列の構築に反映されます (lattice_design_row 関数)。また、測定誤差を補償するために、追加の補正項がデザイン行列にカラムとして追加されることがあります (correction_terms 関数)。

外れ値の自動除去

最小二乗フィッティングのロバスト性を高めるため、統計的に外れ値と判断される観測データは自動的に除去されます。フィットが実行された後、各観測データに対する残差 $ \text{residuals}_i = y_i - \hat{y}_i $ が計算されます。また、観測データの重み $w_i$ に基づいて、各観測の標準偏差 $ \sigma_{\text{obs},i} $ が推定されます。もし $ |\text{residuals}_i| > \text{OUTLIER\_SIGMA} \times \sigma_{\text{obs},i} $ （デフォルトで OUTLIER_SIGMA = 3.0）を満たす観測データがあれば、それは外れ値として除去され、残りのデータセットで最小二乗計算が繰り返されます。このプロセスは外れ値がなくなるか、残りの観測データ数がパラメータ数以下になるまで繰り返されます。

格子定数の導出と誤差伝播

最小二乗法で得られたパラメータ $ \hat{\beta} $ は、直接的に逆格子定数（例: $a^{*2}$ や $2b^*c^*\cos\alpha^*$ ）に関連します。これらのパラメータから、逆格子の長さ $a^*, b^*, c^*$ 、角度 $ \alpha^*, \beta^*, \gamma^* $ 、および体積 $V^*$ が計算されます。さらに、これらの逆格子定数から、実格子の長さ $a, b, c$ 、角度 $ \alpha, \beta, \gamma $ 、および体積 $V$ が、結晶学的な関係式を用いて導出されます。例えば、実格子の体積 $V$ は逆格子の体積 $V^*$ の逆数です。

誤差伝播: 導出された格子定数それぞれの標準偏差は、SigmaPropagator クラスを用いた誤差伝播の法則により計算されます。これは数値微分（摂動法）に基づいており、ある関数 $f(x_1, x_2, \dots, x_n)$ の出力の標準偏差 $ \sigma_f $ は、入力変数 $x_i$ の標準偏差 $ \sigma_{x_i} $ と関数 $f$ の各変数に対する偏微分係数 $ \frac{\partial f}{\partial x_i} $ を用いて、次のように近似されます。

\[ \sigma_f \approx \sqrt{\sum_{i=1}^n \left( \frac{\partial f}{\partial x_i} \sigma_{x_i} \right)^2} \]

数値微分は、微小な摂動 $ \delta x_i $ を用いて $ \frac{\partial f}{\partial x_i} \approx \frac{f(x_1, \dots, x_i + \delta x_i, \dots, x_n) - f(x_1, \dots, x_i, \dots, x_n)}{\delta x_i} $ と計算されます。これにより、フィットパラメータの不確かさが最終的な格子定数にどのように伝播するかを評価します。

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

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

numpy: 高度な数値計算、特に配列操作、線形代数演算、統計計算に使用されます。
matplotlib: 計算結果を視覚化するためのプロット生成に使用されます。

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

pip install numpy matplotlib

必要な入力ファイル

入力ファイルはテキスト形式で、一つまたは複数の格子定数計算ジョブを含めることができます。各ジョブは以下の構造に従います。

ジョブタイトル行
- 任意の文字列を記述します。これはレポートのタイトルとして使用されます。
ジョブ設定行
- スペース区切りの 9つの整数 を記述します。
  1. lattice_system (格子系コード):
    - 1: Triclinic (三斜晶)
    - 2: Monoclinic (unique b axis) (単斜晶, b軸特異)
    - 3: Monoclinic (unique c axis) (単斜晶, c軸特異)
    - 4: Orthorhombic (斜方晶)
    - 5: Tetragonal (正方晶)
    - 6: Cubic (立方晶)
    - 7: Rhombohedral (Trigonal) (菱面体晶, 三方晶)
    - 8: Hexagonal (六方晶)
  2. correction_flag_1
  3. correction_flag_2
  4. correction_flag_3
  5. correction_flag_4
  6. correction_flag_5
    - 各補正フラグは 0 (OFF) または 1 (ON) を指定します。
  7. weight_mode (重みモード):
    - 0: $ \sigma(\theta) $ ベースの重み
    - 1: $ \sigma(2\theta) $ ベースの重み
    - 2: 一定重み (1.0)
    - 3: ユーザー指定重み
    - 4: $1/\sigma^2$
  8. position_mode (位置変換モード):
    - 1: $ \sin(\theta) $ (角度は度数)
    - 2: 生の位置データそのまま
    - 3: ゼロシフト補正付き $ \sin(\theta) $
    - 4: $ \sin(\theta/2) $
    - 5: $ \lambda/(2 \times d) $ に類似の変換
    - その他: 幾何条件依存の変換
  9. profile_mode (プロファイルモード): このスクリプトでは使用されませんが、入力形式の一部です。
波長と装置半径行
- スペース区切りの 2つの浮動小数点数 を記述します。
  1. wavelength (X線の波長)
  2. radius (回折装置の半径、または計算に使用されるスケールファクター)
反射データブロック
- 各反射は1行で記述されますが、特定の weight_mode や position_mode の場合、追加のデータ行が続きます。
- 基本行: スペース区切りの 4つの値
  1. h (ミラー指数 h)
  2. k (ミラー指数 k)
  3. l (ミラー指数 l)
  4. raw_position (生の回折位置データ、単位は position_mode による)
- position_mode=3 の場合: 基本行の直後に zero_shift 値（浮動小数点数）を記述する行が続きます。
- weight_mode が 0, 1, 3, 4 の場合: 基本行（または zero_shift 行）の直後に sigma_or_weight 値（浮動小数点数）を記述する行が続きます。これはモードに応じて標準偏差 $ \sigma $ または直接指定する重み $w$ となります。
- 特殊な h 値:
  - h = 2000: この行は現在のジョブの波長と半径を更新する特殊な指示です。この行の l と raw_position の位置に新しい wavelength と radius を記述し、その後の反射データに適用されます。例: 2000 0 1.5418 100.0 は波長を1.5418、半径を100.0に更新します。（コード上は h k l raw_position の形式で読み込まれるため、k と l は無視され、l と raw_position の値がそれぞれ wavelength と radius と解釈されます。）
  - h = 1000 (または h >= 1000 かつ h != 2000): この値を持つ行は、現在のジョブの反射データブロックの終了を示します。この行のデータ自体は解析されず、次の行からジョブ終了コード行が始まります。
ジョブ終了コード行
- 単一の整数を記述します。
  - 0: 全てのジョブが終了したことを示します。プログラムが終了します。
  - 1: 次のジョブが続くことを示します。プログラムはファイル内の次のジョブタイトル行を探します。

入力ファイルの例

# This is a comment, lines starting with # are ignored
Example Job 1 - Monoclinic with zero shift and sigma(2theta) weights
0 0 0 0 0 1 3 0  # Lattice: Monoclinic(b), Corr:None, Weight:sigma(2theta), Pos:Zero-shift sin(theta)
5406 100.0       # Initial wavelength and radius
0 0 20.0         # h k l raw_position (e.g., 2theta in mm)
25               # zero_shift for position_mode=3
05               # sigma_2theta for weight_mode=1
1 0 30.0
30
08
0 1.5418 100.0 # Update wavelength to 1.5418, radius to 100.0
0 1 40.0
20
06
1 0 35.0
28
07
0 0 0         # End of current job
                # Continue to next job

Example Job 2 - Orthorhombic with constant weights
0 0 0 0 0 2 1 0  # Lattice: Orthorhombic, Corr:None, Weight:Constant, Pos:sin(theta)
7093 50.0        # Initial wavelength and radius
0 0 30.5         # h k l raw_position (e.g., 2theta in degrees)
1 0 45.2
0 1 60.1
1 0 50.3
0 0 0         # End of current job
                # All jobs finished

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

lsq_latt2.py は、計算が完了すると以下の種類のファイルを出力します。

レポートファイル (.out)
プロットファイル (.png)

1. レポートファイル

デフォルトでは入力ファイル名に .out 拡張子が付加されたファイル（例: input.txt → input.out）が生成されます。このファイルには、各計算ジョブの詳細な結果が人間が読みやすい形式で記述されます。

レポートの内容は以下のセクションで構成されます。

ジョブ概要 (Job summary / ジョブ概要):
- ジョブのタイトル、格子系、参照波長、使用されたモード（重み、位置変換、プロファイル）、入力された反射数、採用された反射数、除外された反射数、自由度、重み付き残差二乗和 (Weighted SS)、y空間RMS残差 (RMS residual in y) など、計算の全体的な要約。
実格子定数 (Direct lattice constants / 実格子定数):
- 実格子の長さ $a, b, c$ 、角度 $ \alpha, \beta, \gamma $ (度数表記)、コサイン値 $ \cos\alpha, \cos\beta, \cos\gamma $ 、および体積 $V$ の計算値とそれぞれの標準偏差が記載されます。
逆格子定数 (Reciprocal lattice constants / 逆格子定数):
- 逆格子の長さ $a^*, b^*, c^*$ 、角度 $ \alpha^*, \beta^*, \gamma^* $ (度数表記)、コサイン値 $ \cos\alpha^*, \cos\beta^*, \cos\gamma^* $ 、および体積 $V^*$ の計算値とそれぞれの標準偏差が記載されます。
格子フィットパラメータ (Lattice fit parameters / 格子フィットパラメータ):
- 最小二乗法で直接フィットされたパラメータ（例: $a^{*2}$ , $b^{*2}$ , $c^{*2}$ 、または関連するクロス項）の値と、それぞれの標準偏差がリストされます。
補正モデル (Correction model / 補正モデル):
- 各補正項 (C1〜C5) について、その式、英語と日本語による推定される物理的意味が説明されます。有効化された補正項については、フィットされた係数とその標準偏差も記載されます。
反射一覧 (Reflection table / 反射一覧):
- 最小二乗法で採用された各反射（hkl）について、通し番号、観測された $d$ 値 (d_obs)、計算された $d$ 値 (d_calc)、生の回折位置 (raw_pos)、計算された $2\theta$ (2theta_calc)、y 空間での残差 (residual_y)、 $2\theta$ 空間でのずれ (delta_2theta)、およびその反射に適用された重み (weight) が表形式でリストされます。
除外された反射 (Rejected reflections / 除外された反射):
- 外れ値として最小二乗計算から除外された反射のリストです。外れ値判定基準（デフォルトでは $ |\text{residual}| > 3.0 \times \sigma_{\text{obs}} $ ）が示され、除外された各反射の通し番号、hkl、生の回折位置、および除外理由が記載されます。

2. プロットファイル

各ジョブに対して、Q-1/dプロットが1枚のPNG画像ファイルとして生成されます。ファイル名は、出力ファイル名をベースに、ジョブ番号とサニタイズされたジョブタイトルが追加されます（例: input_job01_sanitizedjobtitle.png）。

プロットの内容は以下の通りです。

散乱ベクトル $Q$ vs 逆格子距離 $1/d$ :
- 横軸に散乱ベクトル $Q = 2\pi/d$ 、縦軸に逆格子距離 $1/d$ をとります。
- 理論的な $1/d = Q / 2\pi$ の関係を示す線。
- 最小二乗法で計算された $1/d$ 値 ( $1/d_\text{calc}$ ) が誤差棒付きでプロットされます。
- 観測された $1/d$ 値 ( $1/d_\text{obs}$ ) がプロットされます。
- 右側のY軸には $ \sigma(1/d) $ ( $1/d$ の標準偏差) がプロットされます。
- 各プロットにはジョブタイトルが付与され、グリッドと凡例が表示されます。

これらの出力ファイルを通じて、計算の妥当性と結果を詳細に確認できます。

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

lsq_latt2.py は、以下の基本的なコマンドライン構文で実行されます。

python lsq_latt2.py [-h] [--silent] [--no-show] input_file [output_file]

input_file (必須): 入力データファイルへのパスを指定します。このファイルは、「必要な入力ファイル」セクションで説明されている形式で記述されている必要があります。
output_file (オプション): 計算結果を書き込む出力レポートファイルへのパスを指定します。この引数を省略した場合、プログラムは input_file の拡張子を .out に置き換えたファイル名（例: data.txt の場合は data.out）をデフォルトの出力ファイル名として使用します。
--silent (オプション): このフラグを指定すると、プログラムの開始時および終了時に表示されるバナーメッセージや、各ジョブのサマリーのコンソール出力が抑制されます。レポートファイルへの出力は通常通り行われます。
--no-show (オプション): このフラグを指定すると、各ジョブの結果として生成されるプロット画像が自動的に新しいウィンドウで表示されなくなります。画像ファイル自体は指定された出力ディレクトリ（またはデフォルトの場所）に保存されます。
-h, --help (オプション): ヘルプメッセージを表示し、プログラムを終了します。

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

ここでは、架空の入力ファイル monoclinic_sample.txt を用いた具体的な使用例と、その実行結果について説明します。

`monoclinic_sample.txt` の内容例

# 単斜晶 (Monoclinic, unique b) のテストデータ
Monoclinic Test Job 1
0 0 0 0 0 1 1 0  # Monoclinic(b), No Corr, Weight:sigma(2theta), Pos:sin(theta)
5406 100.0       # Wavelength and Radius
0 0 25.0         # h k l 2theta(deg)
1                # sigma_2theta
1 0 35.0
15
0 1 45.0
2
1 0 30.0
12
0 0 0         # End of current job
                # All jobs finished

この入力ファイルは、単斜晶系（b軸特異）の格子定数計算ジョブを1つ定義しています。補正項は使用せず、 $ \sigma(2\theta) $ に基づく重み付けを行い、 $2\theta$ から $ \sin(\theta) $ へ変換して処理します。

基本的な実行コマンドと結果

以下のコマンドでプログラムを実行します。

python lsq_latt2.py monoclinic_sample.txt

実行時のコンソール出力 (例)

Least Square Calculation for Lattice Constant
Refactored Python version
Input : monoclinic_sample.txt
Output: monoclinic_sample.out

[Job summary / ジョブ概要]
Lattice system / 格子系                : Monoclinic(unique b) / 単斜晶(b軸特異)
Reference wavelength / 参照波長        : 1.5406000
Weight mode / 重みモード              : 1 (sigma(2theta)-based weight / sigma(2theta) ベース重み)
Position mode / 位置変換モード        : 1 (sin(theta) from angle in degree / 角度(deg)から sin(theta) へ変換)
Profile mode / プロファイルモード      : 0
Input observations / 入力反射数        : 4
Used observations / 採用反射数         : 4
Rejected observations / 除外反射数     : 0
Degrees of freedom / 自由度            : 0
Weighted SS / 重み付き残差二乗和       : 0.0000000e+00
RMS residual in y / y空間RMS残差       : 0.0000000e+00
... (以下、詳細な計算結果が続く) ...
Lattice Constant Calculation was finished!
Saved plot files:
  monoclinic_sample_job01_Monoclinic_Test_Job_1.png

生成されるファイル

monoclinic_sample.out: このファイルは、上記「生成される出力ファイル」セクションで説明されている形式の詳細なレポートを含みます。計算された実格子・逆格子定数、それらの標準偏差、フィットされたパラメータ、各反射のフィット結果、および除外された外れ値（この例ではなし）などが網羅的に記述されます。
monoclinic_sample_job01_Monoclinic_Test_Job_1.png: 散乱ベクトル $Q$ 対 $1/d$ のプロット画像ファイルです。計算された $1/d$ (誤差棒付き) と観測された $1/d$ が比較され、理論曲線も描画されます。この画像はデフォルトで新しいウィンドウに表示されます。

その他の使用例

出力ファイルを指定し、プロットウィンドウを非表示にする:
```
python lsq_latt2.py monoclinic_sample.txt custom_report.txt --no-show
```
このコマンドは、custom_report.txt という名前のレポートファイルを生成し、プロット画像は custom_report_job01_Monoclinic_Test_Job_1.png として保存されますが、プロットウィンドウは表示されません。
バナーメッセージとプロットウィンドウの両方を抑制する:
```
python lsq_latt2.py monoclinic_sample.txt --silent --no-show
```
この場合、コンソールには Saved plot files: の情報のみが表示され、その他のメッセージやプロットウィンドウの表示は抑制されます。レポートファイルとプロットファイルは通常通り生成されます。