crystal_XRD.py 技術ドキュメント

プログラムの動作

本プログラム crystal_XRD.py は、与えられた結晶の格子定数とX線波長に基づき、X線回折(XRD)のBragg角(2θ)と格子面間隔(\(d\))を計算し、対応するミラー指数 \((hkl)\) と共に一覧表示することを目的としています。

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

  • 結晶の格子定数から、実空間および逆空間の格子ベクトル、メトリックテンソル、単位胞体積を計算し、詳細な情報を標準出力に表示します。

  • 指定されたX線波長と最大2θ値の範囲内で、可能な全てのミラー指数 \((hkl)\) に対するBragg角と格子面間隔を計算します。

  • 計算された回折ピークは、2θの昇順にソートされて標準出力に一覧表示されます。

このプログラムは、特定の結晶構造を持つ物質について、その理論的なX線回折パターンにおけるピークの位置を予測するのに役立ちます。これにより、実験的に得られたXRDデータとの比較や、未知試料の構造同定の基礎情報として利用できます。

原理

プログラムは主に以下の物理式と幾何学的計算に基づいています。

  1. 格子定数と格子ベクトル: 実空間の格子定数 \((a, b, c, \alpha, \beta, \gamma)\) から、単位胞の基本格子ベクトル \(\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3\) を計算します。一般的な斜交晶の場合には、以下の変換が用いられます(tkcrystalbase モジュール内の cal_lattice_vectors 関数で実装されていると推測されます)。 $\( \begin{aligned} \mathbf{a}_1 &= (a, 0, 0) \\ \mathbf{a}_2 &= (b \cos\gamma, b \sin\gamma, 0) \\ \mathbf{a}_3 &= (c \cos\beta, c (\cos\alpha - \cos\beta \cos\gamma) / \sin\gamma, c \sqrt{1 - \cos^2\alpha - \cos^2\beta - \cos^2\gamma + 2 \cos\alpha \cos\beta \cos\gamma} / \sin\gamma) \end{aligned} \)$

  2. メトリックテンソル: 実空間のメトリックテンソル \(g_{ij}\) は、格子ベクトルを用いて以下のように定義されます(cal_metrics 関数)。 $\( g_{ij} = \mathbf{a}_i \cdot \mathbf{a}_j \)$

  3. 単位胞の体積: 単位胞の体積 \(V\) は、格子ベクトルのスカラー三重積で計算されます(cal_volume 関数)。 $\( V = |\mathbf{a}_1 \cdot (\mathbf{a}_2 \times \mathbf{a}_3)| \)$

  4. 逆格子ベクトルと逆格子定数: 実空間の格子ベクトル \(\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3\) から、逆格子ベクトル \(\mathbf{a}_1^*, \mathbf{a}_2^*, \mathbf{a}_3^*\) を計算します(cal_reciprocal_lattice_vectors 関数)。 $\( \begin{aligned} \mathbf{a}_1^* &= \frac{\mathbf{a}_2 \times \mathbf{a}_3}{V} \\ \mathbf{a}_2^* &= \frac{\mathbf{a}_3 \times \mathbf{a}_1}{V} \\ \mathbf{a}_3^* &= \frac{\mathbf{a}_1 \times \mathbf{a}_2}{V} \end{aligned} \)\( 逆格子定数 \)(a^, b^, c^, \alpha^, \beta^, \gamma^)$ は、これらの逆格子ベクトルの長さと角度から導出されます(cal_reciprocal_lattice_parameters 関数)。

  5. 逆格子メトリックテンソルと逆格子空間の距離: 逆格子メトリックテンソル \(g^{ij}\) は、逆格子ベクトルを用いて以下のように定義されます(cal_metrics 関数)。 $\( g^{ij} = \mathbf{a}_i^* \cdot \mathbf{a}_j^* \)\( ミラー指数 \)(h, k, l)\( に対応する逆格子ベクトル \)\mathbf{G}{hkl}\( は、以下で与えられます。 \)\( \mathbf{G}_{hkl} = h \mathbf{a}_1^* + k \mathbf{a}_2^* + l \mathbf{a}_3^* \)\( 原点からの逆格子ベクトル \)\mathbf{G}{hkl}\( の長さ \)G_{hkl}\( は、逆格子メトリックテンソル \)g^{ij}\( を用いて計算されます(`distance` 関数で実装)。 \)\( G_{hkl} = ||\mathbf{G}_{hkl}|| = \sqrt{\sum_{i,j=1}^3 (h_i)(h_j) g^{ij}} \)\( ここで \)(h_1, h_2, h_3) = (h, k, l)$ です。

  6. 格子面間隔: ミラー指数 \((hkl)\) に対応する格子面間隔 \(d_{hkl}\) は、逆格子ベクトルの長さの逆数として定義されます。 $\( d_{hkl} = \frac{1}{G_{hkl}} \)$

  7. Braggの法則: X線回折のピークはBraggの法則に従って生じます。 $\( n \lambda = 2 d_{hkl} \sin\theta \)\( ここで、\)n\( は回折次数(通常は1と見なされます)、\)\lambda\( はX線波長、\)d_{hkl}\( は格子面間隔、\)\theta\( はBragg角(入射角)です。 プログラムでは \)\sin\theta = \lambda / (2 d_{hkl})\( を計算し、そこから \)2\theta\( を導出しています。 \)\( 2\theta = 2 \arcsin\left(\frac{\lambda}{2 d_{hkl}}\right) \)\( \)\theta$ はラジアンで計算され、最終的に todeg 定数(tkcrystalbase モジュールからインポート)を用いて度数に変換されます。同様に、torad 定数も度からラジアンへの変換に使用されます。

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

本プログラムは以下の非標準ライブラリに依存します。

  • numpy: 数値計算のための基盤ライブラリ。

  • matplotlib: データの可視化ライブラリ。プログラムでは mpl_toolkits.mplot3d がインポートされていますが、実際にプロットは行われていません。

  • tkcrystalbase: プログラムと同じディレクトリに配置されるカスタムライブラリ。結晶学的な計算(格子ベクトル、メトリックテンソル、体積、距離計算など)を提供します。

numpy および matplotlibpip コマンドでインストールできます。

pip install numpy matplotlib

tkcrystalbase.py は、本プログラム crystal_XRD.py の実行パス内に存在する必要があります。このファイルは別途用意される必要があります。

必要な入力ファイル

本プログラムは、外部からの入力ファイルには依存しません。 必要なデータは全て、プログラムのソースコード内に直接記述された変数として定義されています。

  • lattice_parameters: 格子定数 [a, b, c, alpha, beta, gamma] を定義します(\(a, b, c\) はÅ単位、\(\alpha, \beta, \gamma\) は度単位)。 例: lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]

  • sites: 単位胞内の原子サイト情報がリストとして定義されていますが、現在のプログラムバージョンではこの変数は計算に直接使用されていません。将来的な構造因子計算などの拡張のために用意されている可能性があります。

  • wl: X線波長(Å単位)を定義します。 例: wl = 1.5405

  • Gmin: 逆格子空間における原点からの最小距離を定義します。これは、原点に近い(物理的に意味のない)ミラー指数を除外するために使用されます。

  • Q2max: 計算対象とする2θの最大値(度単位)を定義します。

これらのパラメータは、crystal_XRD.py ファイルの冒頭部分を編集することで変更できます。

生成される出力ファイル

本プログラムは、いかなるファイルも生成したり保存したりしません。 すべての計算結果は、標準出力(コンソール)にテキスト形式で表示されます。

出力される内容は以下の通りです。

  1. 初期設定として指定された格子定数の情報。

  2. 計算された実空間の格子ベクトル、メトリックテンソル、単位胞体積。

  3. 計算された逆格子定数、逆格子ベクトル、逆格子メトリックテンソル、逆格子単位胞体積。

  4. 計算に用いられたミラー指数 \((hkl)\) の探索範囲。

  5. 主要な出力: Bragg角(2θ)、格子面間隔(\(d\))、および対応するミラー指数 \((hkl)\) のリスト。このリストは、2θの昇順にソートされて表示されます。

出力の例(一部):

Lattice parameters: [5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
...
Diffraction angle, d, h, k, l:
  2Q=    15.7533  d=  5.620000  (  1   0   0)
  2Q=    15.7533  d=  5.620000  (  0   1   0)
  2Q=    15.7533  d=  5.620000  (  0   0   1)
  2Q=    22.3364  d=  3.974987  (  1   1   0)
  2Q=    22.3364  d=  3.974987  (  1   0   1)
  ...

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

本プログラム crystal_XRD.py はコマンドライン引数を一切取りません。実行するだけで、内部に設定されたパラメータに基づいて計算を行い、結果を標準出力に出力します。

基本的な実行コマンドは以下の通りです。

python crystal_XRD.py

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

crystal_XRD.py ファイルと、それに依存する tkcrystalbase.py ファイルが同じディレクトリにある状態で、コマンドラインまたはターミナルから以下のコマンドを実行します。

python crystal_XRD.py

実行結果の説明:

(デフォルト設定である lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0] (立方晶) と wl = 1.5405 (Cu Kα1線相当) の場合)

プログラムはまず、入力された格子情報と、それから導出される実空間および逆空間の格子定数、格子ベクトル、メトリックテンソル、体積といった詳細な情報を表示します。

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.3006 A^3

Unit cell volume:    177.3006 A^3
Reciprocal lattice parameters: [0.1779359430604982, 0.1779359430604982, 0.1779359430604982, 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.03166,        0,        0) A^-1
        (        0,  0.03166,        0) A^-1
        (        0,        0,  0.03166) A^-1
Reciprocal unit cell volume:  0.005640 A^-3

hkl range: 3 3 3

続いて、Braggの法則に基づいて計算されたX線回折ピークのリストが表示されます。このリストは2θの値が小さいものから順にソートされており、各ピークに対して以下の情報が含まれます。

  • 2Q: Bragg角 2θ (度)

  • d: 格子面間隔 \(d\) (Å)

  • h, k, l: 対応するミラー指数

デフォルト設定では、立方晶の格子定数が指定されているため、面心立方構造(NaCl型など)の回折ピークに似たパターンが出力されます(例:(100) や (110) など)。

Diffraction angle, d, h, k, l:
  2Q=    15.7533  d=  5.620000  (  1   0   0)
  2Q=    15.7533  d=  5.620000  (  0   1   0)
  2Q=    15.7533  d=  5.620000  (  0   0   1)
  2Q=    22.3364  d=  3.974987  (  1   1   0)
  2Q=    22.3364  d=  3.974987  (  1   0   1)
  2Q=    22.3364  d=  3.974987  (  0   1   1)
  2Q=    27.3995  d=  3.250766  (  1   1   1)
  2Q=    31.9686  d=  2.790757  (  2   0   0)
  2Q=    31.9686  d=  2.790757  (  0   2   0)
  2Q=    31.9686  d=  2.790757  (  0   0   2)
  2Q=    36.1951  d=  2.477742  (  2   1   0)
  2Q=    36.1951  d=  2.477742  (  1   2   0)
  2Q=    36.1951  d=  2.477742  (  2   0   1)
  2Q=    36.1951  d=  2.477742  (  0   2   1)
  2Q=    36.1951  d=  2.477742  (  1   0   2)
  2Q=    36.1951  d=  2.477742  (  0   1   2)
  2Q=    40.1650  d=  2.245999  (  2   1   1)
  2Q=    40.1650  d=  2.245999  (  1   2   1)
  2Q=    40.1650  d=  2.245999  (  1   1   2)
  2Q=    43.9930  d=  2.060136  (  2   2   0)
  2Q=    43.9930  d=  2.060136  (  2   0   2)
  2Q=    43.9930  d=  2.060136  (  0   2   2)
  2Q=    47.7029  d=  1.910321  (  3   0   0)
  2Q=    47.7029  d=  1.910321  (  0   3   0)
  2Q=    47.7029  d=  1.910321  (  0   0   3)
  2Q=    51.3140  d=  1.785897  (  2   2   1)
  2Q=    51.3140  d=  1.785897  (  2   1   2)
  2Q=    51.3140  d=  1.785897  (  1   2   2)
  2Q=    54.8398  d=  1.681144  (  3   1   0)
  2Q=    54.8398  d=  1.681144  (  1   3   0)
  2Q=    54.8398  d=  1.681144  (  3   0   1)
  2Q=    54.8398  d=  1.681144  (  0   3   1)
  2Q=    54.8398  d=  1.681144  (  1   0   3)
  2Q=    54.8398  d=  1.681144  (  0   1   3)
  2Q=    58.2903  d=  1.591557  (  3   1   1)
  2Q=    58.2903  d=  1.591557  (  1   3   1)
  2Q=    58.2903  d=  1.591557  (  1   1   3)
  2Q=    61.6749  d=  1.513567  (  2   2   2)
  2Q=    65.0028  d=  1.444652  (  3   2   0)
  2Q=    65.0028  d=  1.444652  (  2   3   0)
  2Q=    65.0028  d=  1.444652  (  3   0   2)
  2Q=    65.0028  d=  1.444652  (  0   3   2)
  2Q=    65.0028  d=  1.444652  (  2   0   3)
  2Q=    65.0028  d=  1.444652  (  0   2   3)
  2Q=    68.2844  d=  1.383186  (  3   2   1)
  2Q=    68.2844  d=  1.383186  (  2   3   1)
  2Q=    68.2844  d=  1.383186  (  3   1   2)
  2Q=    68.2844  d=  1.383186  (  1   3   2)
  2Q=    68.2844  d=  1.383186  (  2   1   3)
  2Q=    68.2844  d=  1.383186  (  1   2   3)
  2Q=    71.5284  d=  1.328328  (  3   2   2)
  2Q=    71.5284  d=  1.328328  (  2   3   2)
  2Q=    71.5284  d=  1.328328  (  2   2   3)
  2Q=    74.7431  d=  1.279261  (  3   3   0)
  2Q=    74.7431  d=  1.279261  (  3   0   3)
  2Q=    74.7431  d=  1.279261  (  0   3   3)
  2Q=    77.9362  d=  1.235123  (  3   3   1)
  2Q=    77.9362  d=  1.235123  (  3   1   3)
  2Q=    77.9362  d=  1.235123  (  1   3   3)
  2Q=    81.1143  d=  1.195048  (  3   3   2)
  2Q=    81.1143  d=  1.195048  (  3   2   3)
  2Q=    81.1143  d=  1.195048  (  2   3   3)
  2Q=    84.2825  d=  1.158308  (  4   0   0)
  2Q=    84.2825  d=  1.158308  (  0   4   0)
  2Q=    84.2825  d=  1.158308  (  0   0   4)
  2Q=    87.4452  d=  1.124233  (  4   1   0)
  2Q=    87.4452  d=  1.124233  (  1   4   0)
  2Q=    87.4452  d=  1.124233  (  4   0   1)
  2Q=    87.4452  d=  1.124233  (  0   4   1)
  2Q=    87.4452  d=  1.124233  (  1   0   4)
  2Q=    87.4452  d=  1.124233  (  0   1   4)
  2Q=    90.6063  d=  1.092289  (  4   1   1)
  2Q=    90.6063  d=  1.092289  (  1   4   1)
  2Q=    90.6063  d=  1.092289  (  1   1   4)
  2Q=    93.7691  d=  1.062143  (  4   2   0)
  2Q=    93.7691  d=  1.062143  (  2   4   0)
  2Q=    93.7691  d=  1.062143  (  4   0   2)
  2Q=    93.7691  d=  1.062143  (  0   4   2)
  2Q=    93.7691  d=  1.062143  (  2   0   4)
  2Q=    93.7691  d=  1.062143  (  0   2   4)
  2Q=    96.9372  d=  1.033534  (  3   3   3)
  2Q=    96.9372  d=  1.033534  (  4   2   1)
  2Q=    96.9372  d=  1.033534  (  2   4   1)
  2Q=    96.9372  d=  1.033534  (  4   1   2)
  2Q=    96.9372  d=  1.033534  (  1   4   2)
  2Q=    96.9372  d=  1.033534  (  2   1   4)
  2Q=    96.9372  d=  1.033534  (  1   2   4)
  2Q=   100.1139  d=  1.006249  (  4   2   2)
  2Q=   100.1139  d=  1.006249  (  2   4   2)
  2Q=   100.1139  d=  1.006249  (  2   2   4)
  2Q=   103.3023  d=  0.980145  (  4   3   0)
  2Q=   103.3023  d=  0.980145  (  3   4   0)
  2Q=   103.3023  d=  0.980145  (  4   0   3)
  2Q=   103.3023  d=  0.980145  (  0   4   3)
  2Q=   103.3023  d=  0.980145  (  3   0   4)
  2Q=   103.3023  d=  0.980145  (  0   3   4)
  2Q=   106.5057  d=  0.955097  (  4   3   1)
  2Q=   106.5057  d=  0.955097  (  3   4   1)
  2Q=   106.5057  d=  0.955097  (  4   1   3)
  2Q=   106.5057  d=  0.955097  (  1   4   3)
  2Q=   106.5057  d=  0.955097  (  3   1   4)
  2Q=   106.5057  d=  0.955097  (  1   3   4)
  2Q=   109.7289  d=  0.931005  (  4   3   2)
  2Q=   109.7289  d=  0.931005  (  3   4   2)
  2Q=   109.7289  d=  0.931005  (  4   2   3)
  2Q=   109.7289  d=  0.931005  (  2   4   3)
  2Q=   109.7289  d=  0.931005  (  3   2   4)
  2Q=   109.7289  d=  0.931005  (  2   3   4)
  2Q=   112.9754  d=  0.907797  (  4   4   0)
  2Q=   112.9754  d=  0.907797  (  4   0   4)
  2Q=   112.9754  d=  0.907797  (  0   4   4)
  2Q=   116.2483  d=  0.885408  (  4   4   1)
  2Q=   116.2483  d=  0.885408  (  4   1   4)
  2Q=   116.2483  d=  0.885408  (  1   4   4)
  2Q=   119.5508  d=  0.863778  (  4   4   2)
  2Q=   119.5508  d=  0.863778  (  4   2   4)
  2Q=   119.5508  d=  0.863778  (  2   4   4)
  2Q=   122.8858  d=  0.842858  (  3   3   4)
  2Q=   122.8858  d=  0.842858  (  3   4   3)
  2Q=   122.8858  d=  0.842858  (  4   3   3)
  2Q=   126.2562  d=  0.822596  (  4   4   3)
  2Q=   126.2562  d=  0.822596  (  4   3   4)
  2Q=   126.2562  d=  0.822596  (  3   4   4)
  2Q=   129.6644  d=  0.802952  (  5   0   0)
  2Q=   129.6644  d=  0.802952  (  0   5   0)
  2Q=   129.6644  d=  0.802952  (  0   0   5)
  2Q=   133.1132  d=  0.783888  (  5   1   0)
  2Q=   133.1132  d=  0.783888  (  1   5   0)
  2Q=   133.1132  d=  0.783888  (  5   0   1)
  2Q=   133.1132  d=  0.783888  (  0   5   1)
  2Q=   133.1132  d=  0.783888  (  1   0   5)
  2Q=   133.1132  d=  0.783888  (  0   1   5)
  2Q=   136.6052  d=  0.765369  (  5   1   1)
  2Q=   136.6052  d=  0.765369  (  1   5   1)
  2Q=   136.6052  d=  0.765369  (  1   1   5)
  2Q=   140.1430  d=  0.747363  (  5   2   0)
  2Q=   140.1430  d=  0.747363  (  2   5   0)
  2Q=   140.1430  d=  0.747363  (  5   0   2)
  2Q=   140.1430  d=  0.747363  (  0   5   2)
  2Q=   140.1430  d=  0.747363  (  2   0   5)
  2Q=   140.1430  d=  0.747363  (  0   2   5)
  2Q=   143.7299  d=  0.729841  (  5   2   1)
  2Q=   143.7299  d=  0.729841  (  2   5   1)
  2Q=   143.7299  d=  0.729841  (  5   1   2)
  2Q=   143.7299  d=  0.729841  (  1   5   2)
  2Q=   143.7299  d=  0.729841  (  2   1   5)
  2Q=   143.7299  d=  0.729841  (  1   2   5)
  2Q=   147.3644  d=  0.712773  (  5   2   2)
  2Q=   147.3644  d=  0.712773  (  2   5   2)
  2Q=   147.3644  d=  0.712773  (  2   2   5)

この出力により、指定された結晶構造に対する理論的なX線回折ピークの位置を簡単に確認できます。プログラムの動作や出力は、crystal_XRD.py ファイル内の lattice_parameters, wl, Q2max などの変数を編集することでカスタマイズ可能です。