crystal_XRD.py 技術ドキュメント
プログラムの動作
crystal_XRD.py は、指定された結晶構造（格子定数と原子サイト情報）とX線波長に基づき、理論的なX線回折パターンを計算するPythonスクリプトです。このプログラムは、実空間および逆空間の格子情報（格子ベクトル、計量テンソル、単位胞体積）を計算し、指定された最大回折角 (\(2\theta_{max}\)) までのすべての可能な\(hkl\)指数に対する回折角 (\(2\theta\)) と格子面間隔 (\(d\)) を算出します。最終的に、計算された回折ピークは\(2\theta\)値の昇順にソートされ、標準出力に表示されます。
主な機能:
結晶の格子定数から実空間の格子ベクトル、計量テンソル、単位胞体積を計算します。
実空間の格子情報から逆空間の格子定数、格子ベクトル、計量テンソル、単位胞体積を計算します。
指定されたX線波長と\(2\theta_{max}\)に基づき、回折可能な\(hkl\)反射を網羅的に探索します。
各\(hkl\)反射について、逆格子空間での距離\(G\)を計算し、それを用いて格子面間隔\(d\)を決定します。
格子面間隔\(d\)とX線波長からブラッグの法則を用いて回折角\(2\theta\)を計算します。
計算された\(2\theta\)と\(d\)、対応する\(hkl\)指数を\(2\theta\)の昇順にソートして出力します。
解決する課題: このプログラムは、既知の結晶構造からX線回折実験で観測されるべき回折ピークの位置を予測する際に役立ちます。これにより、実験データの解析や、特定の結晶相の同定のための理論的基盤を提供します。
原理
crystal_XRD.py は、結晶学の基本原理とブラッグの法則に基づいて回折パターンを計算します。
1. 格子定数と格子ベクトル
結晶は6つの格子定数 \((a, b, c, \alpha, \beta, \gamma)\) で定義されます。これらの定数から、実空間の基本格子ベクトル \(\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3\) が導出されます。
プログラムでは tkcrystalbase モジュールの cal_lattice_vectors 関数を使用してこれらを計算します。例えば、直交座標系で表現されるこれらのベクトルは以下のようになります。
ここで \(V_{unit}\) は単位胞の体積です。
2. 計量テンソルと単位胞体積
実空間の計量テンソル \(g_{ij}\) は、格子ベクトル間の内積によって定義されます。
プログラムでは tkcrystalbase モジュールの cal_metrics 関数で計算されます。
単位胞の体積 \(V\) は、格子ベクトルのスカラー三重積から計算されます。
プログラムでは tkcrystalbase モジュールの cal_volume 関数で計算されます。
3. 逆格子ベクトルと逆格子計量テンソル
回折現象は逆格子空間で記述されます。実空間の格子ベクトルから逆格子ベクトル \(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3\) が導出されます。これらのベクトルは、実格子ベクトルに対して以下の関係を満たします（ここでは係数 \(2\pi\) を含まない定義を採用し、後の \(G=1/d\) と整合させます）。
プログラムでは tkcrystalbase モジュールの cal_reciprocal_lattice_vectors 関数で計算されます。
逆格子計量テンソル \(g^*_{ij}\) は、逆格子ベクトル間の内積として定義されます。
これも tkcrystalbase モジュールの cal_metrics 関数（逆格子パラメータ適用時）で計算されます。
4. 逆格子空間における距離と格子面間隔
任意のリフレクション \(hkl\) は、逆格子ベクトル \(\mathbf{G}_{hkl}\) で表されます。
このベクトルの大きさ \(G = |\mathbf{G}_{hkl}|\) は、対応する格子面間隔 \(d\) の逆数に等しくなります。
プログラムでは tkcrystalbase モジュールの distance 関数がこの \(G\) の値を計算します。この \(G\) は、原点 \((0,0,0)\) から逆格子点 \((h,k,l)\) までの距離であり、逆格子計量テンソル \(g^*_{ij}\) を用いて以下のように計算されます。
ここで \(h_i\) は \((h,k,l)\) です。
5. ブラッグの法則と回折角
X線回折の原理はブラッグの法則で記述されます。
ここで \(d\) は格子面間隔、\(n\) は回折次数（通常 \(n=1\) と仮定）、\(\lambda\) はX線波長、\(\theta\) はブラッグ角です。
プログラムでは \(d = 1/G\) を用いて \(\sin\theta = \frac{\lambda}{2d}\) を計算し、それから \(2\theta = 2 \arcsin\left(\frac{\lambda}{2d}\right)\) を導出します。
計算された \(2\theta\) が指定された Q2max を超える場合はスキップされます。
Gmin は、逆格子空間の原点に近い（つまり \(d\) が非常に大きい）計算上の「回折」を除外するために使用されます。
必要な非標準ライブラリとインストール方法
本プログラムの実行には、以下の非標準ライブラリが必要です。
numpy: 数値計算を効率的に行うためのライブラリ。matplotlib: グラフ描画のためのライブラリ。mpl_toolkits.mplot3dはmatplotlibの一部として3Dプロット機能を提供しますが、本プログラムでは直接使用されていません。tkcrystalbase: 結晶学計算に関するユーティリティ関数を提供するカスタムモジュール。このモジュールは別途提供されるか、プログラムと同じディレクトリ、またはPythonの検索パスに配置されている必要があります。
numpy と matplotlib は pip コマンドでインストールできます。
pip install numpy matplotlib
必要な入力ファイル
本プログラムは、コマンドラインからのファイル入力を直接は受け付けません。必要な結晶構造のデータ（格子定数、原子サイト情報、X線波長、\(2\theta_{max}\)などの計算パラメータ）は、プログラムのソースコード内に直接記述されています。
tkcrystalbase.py: 結晶計算のユーティリティ関数が含まれているPythonモジュールです。このファイルがcrystal_XRD.pyと同じディレクトリに存在するか、Pythonのモジュール検索パス上に存在する必要があります。
ユーザーは、lattice_parameters 変数や wl (X線波長) などの値をソースコード内で編集することで、異なる結晶や条件での計算を行うことができます。
生成される出力ファイル
crystal_XRD.py は、いかなるファイルも生成しません。すべての計算結果は標準出力 (stdout) にテキスト形式で表示されます。出力には、結晶の格子情報、逆格子情報、そして計算された各\(hkl\)反射に対する\(2\theta\)角、格子面間隔\(d\)、および対応する\(hkl\)指数が含まれます。
コマンドラインでの使用例 (Usage)
crystal_XRD.py は、引数をとらずに直接Pythonインタープリタで実行されます。
python crystal_XRD.py
コマンドラインでの具体的な使用例
以下は、プログラムのデフォルト設定（NaCl型の立方晶構造、格子定数 \(a=5.62\) Å、X線波長 \(\lambda=1.5405\) Å、最大 \(2\theta=150\) 度）での実行例とその結果の説明です。
実行コマンド:
python crystal_XRD.py
実行結果（抜粋）:
Lattice parameters: [5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
Lattice vectors:
ax: ( 5.62, 0, 0) A
ay: ( 0, 5.62, 0) A
az: ( 0, 0, 5.62) A
Metric tensor:
gij: ( 31.58, 0, 0) A
( 0, 31.58, 0) A
( 0, 0, 31.58) A
Volume: 177.3 A^3
Unit cell volume: 177.3 A^3
Reciprocal lattice parameters: [0.17793594306049822, 0.17793594306049822, 0.17793594306049822, 90.0, 90.0, 90.0]
Reciprocal lattice vectors:
Rax: ( 0.1779, 0, 0) A^-1
Ray: ( 0, 0.1779, 0) A^-1
Raz: ( 0, 0, 0.1779) A^-1
Reciprocal lattice metric tensor:
Rgij: ( 0.0317, 0, 0) A^-1
( 0, 0.0317, 0) A^-1
( 0, 0, 0.0317) A^-1
Reciprocal unit cell volume: 0.00564 A^-3
hkl range: 3 3 3
Diffraction angle, d, h, k, l:
2Q= 23.7081 d= 3.766863 ( 1 1 1)
2Q= 27.4208 d= 3.243500 ( 2 0 0)
2Q= 39.0664 d= 2.298288 ( 2 2 0)
2Q= 46.2201 d= 1.961026 ( 3 1 1)
2Q= 49.6105 d= 1.836102 ( 2 2 2)
2Q= 57.2647 d= 1.609501 ( 4 0 0)
2Q= 65.7337 d= 1.419163 ( 3 3 1)
2Q= 69.3093 d= 1.354026 ( 4 2 0)
2Q= 76.9404 d= 1.238478 ( 4 2 2)
2Q= 83.3364 d= 1.157835 ( 5 1 1)
2Q= 83.3364 d= 1.157835 ( 3 3 3)
2Q= 90.1783 d= 1.082750 ( 4 4 0)
2Q= 95.6989 d= 1.026731 ( 5 3 1)
2Q= 99.7126 d= 0.988220 ( 6 0 0)
2Q= 99.7126 d= 0.988220 ( 4 4 2)
2Q= 104.9922 d= 0.941912 ( 6 2 0)
2Q= 112.5925 d= 0.880482 ( 5 3 3)
2Q= 112.5925 d= 0.880482 ( 6 2 2)
2Q= 119.8631 d= 0.829033 ( 4 4 4)
2Q= 126.3986 d= 0.787680 ( 5 5 1)
2Q= 126.3986 d= 0.787680 ( 7 1 1)
2Q= 131.6258 d= 0.757820 ( 6 4 0)
2Q= 135.4057 d= 0.737119 ( 7 3 1)
2Q= 138.8340 d= 0.719602 ( 6 4 2)
2Q= 141.9701 d= 0.704987 ( 5 5 3)
2Q= 141.9701 d= 0.704987 ( 7 3 3)
説明:
出力の冒頭では、設定された格子定数から計算された実空間の格子ベクトル、計量テンソル、および単位胞体積が表示されます。続いて、これらの実空間情報から導出された逆格子パラメータ、逆格子ベクトル、逆格子計量テンソル、および逆格子単位胞体積が表示されます。
hkl range は、計算対象となる\(h, k, l\)指数の最大絶対値を示します。
その後に続く Diffraction angle, d, h, k, l: のセクションでは、ブラッグの法則と結晶構造に基づいて計算されたすべての回折ピークが、\(2\theta\)の昇順にリストされます。各行は、回折角 \(2\theta\)（度）、格子面間隔 \(d\)（Å）、および対応するミラー指数 \((h, k, l)\) を示します。この例では、立方晶の主要な回折ピークがリストされており、同じ\(2\theta\)を持つ複数の\(hkl\)反射（例: (333) と (511)）も適切にソートされて表示されています。
プログラムは、Q2max (150度) で指定された最大角を超えるピークや、ブラッグの法則で物理的に不可能なピーク (\(\sin\theta \ge 1.0\)) は出力しません。