test_lsq_latt2.py

プログラムの動作

test_lsq_latt2.py は、結晶学における格子定数決定のための最小二乗フィットプログラム lsq_latt2.py の機能を検証するためのテストスクリプトです。このプログラムは、以下の手順で lsq_latt2.py の動作を総合的に評価します。

理論的な回折ピークの生成: 三斜晶系、単斜晶系、斜方晶系、正方晶系、立方晶系、三方晶系、六方晶系といった、様々な結晶系の理想的な格子定数と指定されたX線波長を用いて、理論的な回折ピーク（Miller指数 \((hkl)\) と \(2\theta\) 値）を計算します。
ノイズの付加: 計算された \(2\theta\) 値にランダムなガウスノイズを付加し、実際の実験データに近い状況をシミュレートします。
入力ファイルの生成: 生成されたノイズ付き回折ピークデータから、lsq_latt2.py が処理できる形式の入力ファイルを動的に作成します。この入力ファイルは一時ディレクトリ内に保存されます。
lsq_latt2.py の実行: 作成された入力ファイルを引数として、指定された lsq_latt2.py スクリプトを subprocess として実行します。
出力の解析: lsq_latt2.py が生成した出力ファイルから、フィットによって得られた格子定数を解析・抽出します。
結果の評価: 抽出された格子定数を、元の理想的な格子定数と比較し、その誤差が許容範囲内であるかを判定します。
レポートの生成: 各テストケースの結果（成功/失敗、誤差の詳細、推測される問題点）を標準出力に表示し、テスト全体のサマリーを出力します。

このプログラムは、lsq_latt2.py が様々な結晶系、ノイズレベル、回折ピーク数に対してロバストかつ正確に動作することを確認することを目的としています。

原理

test_lsq_latt2.py は、X線回折における結晶学の基本原理と線形代数の概念を用いて、回折ピークデータを生成します。

計量テンソル (Metric Tensor): 結晶の格子定数 \(a, b, c, \alpha, \beta, \gamma\) から、実格子 (direct lattice) の計量テンソル \(G\) が計算されます。

\[\begin{split} G = \begin{pmatrix} a^2 & ab\cos\gamma & ac\cos\beta \\ ab\cos\gamma & b^2 & bc\cos\alpha \\ ac\cos\beta & bc\cos\alpha & c^2 \end{pmatrix} \end{split}\]

この \(G\) は、格子ベクトル間の内積を表現します。
逆格子計量テンソル (Reciprocal Metric Tensor): 回折現象を記述するために、逆格子 (reciprocal lattice) の概念が用いられます。逆格子計量テンソル \(G^*\) は実格子計量テンソル \(G\) の逆行列として定義されます。

\[ G^* = G^{-1} \]

プログラムでは inv3 関数により3x3行列の逆行列を計算しています。逆行列の計算にはサラスの公式が用いられます。
格子面間隔 (d-spacing): 特定のMiller指数 \((hkl)\) で指定される格子面間隔 \(d_{hkl}\) は、逆格子ベクトル \(\mathbf{s}_{hkl} = h\mathbf{a}^* + k\mathbf{b}^* + l\mathbf{c}^*\) の長さの逆数として得られます。ここで、\(\mathbf{a}^*, \mathbf{b}^*, \mathbf{c}^*\) は逆格子基本ベクトルです。 \(d_{hkl}^{-2}\) は、Miller指数ベクトル \(\mathbf{v} = (h, k, l)^T\) と逆格子計量テンソル \(G^*\) を用いて以下のように計算されます。

\[ d_{hkl}^{-2} = \mathbf{v}^T G^* \mathbf{v} = h^2a^{*2} + k^2b^{*2} + l^2c^{*2} + 2hka^*b^*\cos\gamma^* + 2hla^*c^*\cos\beta^* + 2klb^*c^*\cos\alpha^* \]

ここで、\(a^*, b^*, c^*, \alpha^*, \beta^*, \gamma^*\) は逆格子定数です。プログラムでは d_spacing 関数がこの計算を行います。
ブラッグの法則 (Bragg's Law) と \(2\theta\): 格子面間隔 \(d\) とX線波長 \(\lambda\) から、回折角 \(2\theta\) はブラッグの法則 \(2d\sin\theta = \lambda\) に従って計算されます。

\[ \sin\theta = \frac{\lambda}{2d} \]

したがって、

\[ 2\theta = 2 \arcsin\left(\frac{\lambda}{2d}\right) \]

ここで計算される \(\theta\) はラジアン単位ですが、出力される \(2\theta\) は度単位に変換されます。プログラムでは two_theta_from_d 関数がこの計算を行います。

これらの原理に基づいて、特定の結晶系に対する理論的な回折ピーク位置が生成され、その後、ノイズを加えることで実際の実験データに近いシナリオがシミュレートされます。lsq_latt2.py は、このノイズ付き \(2\theta\) データを用いて、元の格子定数を最小二乗法で再構築しようと試みます。

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

test_lsq_latt2.py スクリプト自体は、Pythonの標準ライブラリのみを使用しており、追加の非標準ライブラリを必要としません。

ただし、このテストスクリプトが対象とする lsq_latt2.py は、その内部実装に応じて別途非標準ライブラリを必要とする可能性があります。lsq_latt2.py のドキュメントまたはソースコードを参照して、必要なライブラリを確認してください。

必要な入力ファイル

test_lsq_latt2.py が直接必要とする入力ファイルは、コマンドライン引数として与えられるテスト対象 lsq_latt2.py スクリプトのパスです。

test_lsq_latt2.py は、lsq_latt2.py の実行時に必要な入力ファイルを動的に生成します。これらのファイルは一時ディレクトリ内に作成され、通常はテスト終了後に削除されます（--keep-workdir オプションが指定されていない場合）。

生成される入力ファイルの形式は以下の通りです。

[任意のタイトル行]
[LS] 0 0 0 0 0 2 4 [IP]
[Wavelength (Å)] 0.0
[h1] [k1] [l1] [2theta1_noisy (deg)]
[h2] [k2] [l2] [2theta2_noisy (deg)]
...
1000 0 0 0
0

任意のタイトル行: テストケースの名前など。
LS: 結晶系を指定する整数コード (1: 三斜晶系, 2: 単斜晶系b軸, 3: 単斜晶系c軸, 4: 斜方晶系, 5: 正方晶系, 6: 立方晶系, 7: 三方晶系, 8: 六方晶系)。
IP: 処理モード（通常は 1）。
Wavelength (Å): X線波長をÅ単位で指定。
h k l: Miller指数（整数）。
2theta_noisy (deg): ノイズが加わった \(2\theta\) 値を度単位で指定（小数点以下6桁まで）。
1000 0 0 0: ピークデータリストの終了マーカー。
0: 入力終了マーカー。

このファイルは、各テストケースの一時ディレクトリ内に input.BZ.K の名前で生成されます。

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

test_lsq_latt2.py は、テスト対象の lsq_latt2.py スクリプトに、結果を出力するためのファイルパスを渡します。この出力ファイルは一時ディレクトリ内に作成され、通常は output.txt という名前になります。

test_lsq_latt2.py は、この output.txt ファイルの内容を解析し、フィットによって得られた格子定数を抽出します。解析対象となるのは、lsq_latt2.py の出力に含まれる以下のブロックです。

新フォーマットの場合:

[Direct lattice constants / 実格子定数]
 a             = [value] ± [sigma]
 b             = [value] ± [sigma]
 c             = [value] ± [sigma]
 alpha         = [value] ± [sigma]
 beta          = = [value] ± [sigma]
 gamma         = [value] ± [sigma]
 ...

旧フォーマットの場合:

Direct cell constant
 (some lines)
 a          b          c
 [value(sigma)] [value(sigma)] [value(sigma)]
 (some lines)
 alpha      beta       gamma
 [value(sigma)] [value(sigma)] [value(sigma)]
 ...

ここで [value] はフィットされた格子定数、[sigma] はその標準偏差です。test_lsq_latt2.py はこれらの値（主に [value] の部分）を抽出し、元の理論値と比較してテストの合否を判定します。

また、test_lsq_latt2.py はテストの実行中に、標準出力に各テストケースの実行結果、成功/失敗のステータス、詳細な誤差情報、およびテストのサマリーを出力します。

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

test_lsq_latt2.py は、以下の書式でコマンドラインから実行できます。

python test_lsq_latt2.py [lsq_script] [options]

lsq_script: （オプション）テスト対象の lsq_latt2.py スクリプトへのパスを指定します。省略された場合、デフォルトで ./lsq_latt2.py が使用されます。

オプション:

--wavelength <WAVELENGTH>: X線波長をオングストローム (Å) 単位で指定します。
- デフォルト値: 1.5406 (Cu Kα1)
--n-lines <N_LINES>: 各テストケースで生成する回折ピークの数を指定します。
- デフォルト値: 15
--noise-sigma <NOISE_SIGMA>: \(2\theta\) 値に加えるガウスノイズの標準偏差を度 (degree) 単位で指定します。
- デフォルト値: 0.02
--seed <SEED>: 乱数生成のシード値を指定します。これにより、ノイズの再現性が保証されます。
- デフォルト値: 1234
--keep-workdir: このフラグが指定されると、テスト中に生成された一時作業ディレクトリとその内容がテスト終了後も削除されずに保持されます。デバッグや詳細なログ確認に有用です。

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

1. 基本的な実行例 (デフォルトの `lsq_latt2.py` をテスト)

カレントディレクトリに lsq_latt2.py が存在する場合、最も簡単な実行方法です。

python test_lsq_latt2.py

実行結果の説明: このコマンドは、test_lsq_latt2.py が見つけることができる lsq_latt2.py を対象に、8つの異なる結晶系（三斜晶系、単斜晶系など）に対してテストを実行します。各テストケースでは、X線波長 1.5406 Å、15個の回折ピーク、および \(2\theta\) 値に 0.02 度のガウスノイズを適用します。テスト結果は各ケースごとに [PASS] または [FAIL] と共に、元の格子定数、フィットされた格子定数、および各パラメータの絶対誤差と相対誤差が報告されます。最後に、全体の成功/失敗数がまとめられて表示されます。

2. 特定の `lsq_latt2.py` パス、異なる波長、多くのピーク、高めのノイズ、および作業ディレクトリ保持オプション付きの実行例

/opt/my_scripts/lsq_latt2.py というパスに lsq_latt2.py が存在し、これをテストしたい場合、またテストパラメータを調整したい場合。

python test_lsq_latt2.py /opt/my_scripts/lsq_latt2.py --wavelength 0.7107 --n-lines 20 --noise-sigma 0.05 --keep-workdir

実行結果の説明: このコマンドは、/opt/my_scripts/lsq_latt2.py をテスト対象とします。X線波長はモリブデンKα1に相当する 0.7107 Å に設定され、各テストケースで生成される回折ピーク数は 20 に増えます。\(2\theta\) 値に加わるノイズの標準偏差は 0.05 度と、デフォルトよりも大きくなります。--keep-workdir オプションにより、テスト中に生成された一時ディレクトリ（input.BZ.K や output.txt などが含まれます）はテスト終了後も削除されずに残り、デバッグや結果の確認に利用できます。各ケースの結果の表示形式は基本的な実行例と同様ですが、指定されたパラメータが反映された結果となります。最終的に、テストのサマリーが表示され、一時ディレクトリがどこに保存されたかの情報も出力されます。