技術ドキュメント: crystal_distance.py

プログラムの動作

本プログラム crystal_distance.py は、結晶構造内の原子間距離を計算し、その結果を標準出力に表示するツールです。

目的: 結晶学において、原子間の相互作用を理解するためには正確な原子間距離の計算が不可欠です。本プログラムは、周期的な結晶構造における任意の2つの原子サイト間の距離を、指定された距離範囲内で網羅的に計算することを目的とします。

主な機能:

結晶情報の計算: 設定された格子パラメーター（格子定数と軸角）に基づいて、デカルト座標系における格子ベクトル、メトリックテンソル、単位胞体積、およびそれらの逆格子空間における対応する情報を計算します。
原子間距離の計算: プログラムにハードコードされた原子サイト情報と、定義された距離範囲（$r_{min}$ から $r_{max}$）を使用し、周期境界条件を考慮して原子間の距離を計算します。
結果の表示: 計算されたすべての原子間距離を昇順にソートし、どのサイトがどの並進ベクトルによって関連付けられているかとともに標準出力に出力します。

解決する課題: 結晶は周期的な構造を持つため、特定の原子からの距離を計算する際には、同一単位胞内の原子だけでなく、隣接する単位胞内の原子も考慮する必要があります。本プログラムは、格子ベクトルの整数倍による並進を考慮することで、この周期性を適切に処理し、正確な原子間距離を効率的に特定します。

原理

プログラムは、結晶学の基本的な概念と数式を用いて原子間距離を計算します。

格子ベクトル (Lattice Vectors) 結晶の単位胞は、格子定数 $a, b, c$ と軸角 $\alpha, \beta, \gamma$ によって特徴づけられます。これらのパラメーターから、デカルト座標系における格子ベクトル $\vec{a}, \vec{b}, \vec{c}$ が以下のように定義されます。

\[\vec{a} = (a, 0, 0)\]

\[\vec{b} = (b \cos\gamma, b \sin\gamma, 0)\]

\[\vec{c} = (c \cos\beta, c \frac{\cos\alpha - \cos\beta \cos\gamma}{\sin\gamma}, c \frac{\sqrt{1 - \cos^2\alpha - \cos^2\beta - \cos^2\gamma + 2\cos\alpha\cos\beta\cos\gamma}}{\sin\gamma})\]

これらの計算は、外部モジュール tkcrystalbase の cal_lattice_vectors 関数によって実行されます。
メトリックテンソル (Metric Tensor) メトリックテンソル $g_{ij}$ は、結晶学において分数座標系での距離計算をデカルト座標系に変換するために使用されます。これは、格子ベクトルの内積で構成される3x3の対称行列です。

\[\begin{split}g = \begin{pmatrix} \vec{a} \cdot \vec{a} & \vec{a} \cdot \vec{b} & \vec{a} \cdot \vec{c} \\ \vec{b} \cdot \vec{a} & \vec{b} \cdot \vec{b} & \vec{b} \cdot \vec{c} \\ \vec{c} \cdot \vec{a} & \vec{c} \cdot \vec{b} & \vec{c} \cdot \vec{c} \end{pmatrix} = \begin{pmatrix} a^2 & ab\cos\gamma & ac\cos\beta \\ ba\cos\gamma & b^2 & bc\cos\alpha \\ ca\cos\beta & cb\cos\alpha & c^2 \end{pmatrix}\end{split}\]

このテンソルは tkcrystalbase.cal_metrics 関数の一部として計算されます。
原子間距離 (Interatomic Distance) 結晶内の2つの原子サイトの分数座標をそれぞれ $\vec{p_0} = (x_0, y_0, z_0)$ と $\vec{p_1} = (x_1, y_1, z_1)$ とします。周期境界条件を考慮するため、サイト1は任意の整数ベクトル $\vec{n} = (ix, iy, iz)$ だけ並進させた位置にあると見なします。この場合、相対座標の差は $\Delta\vec{p} = \vec{p_1} + \vec{n} - \vec{p_0}$ となります。原子間距離 $d$ の二乗は、メトリックテンソル $g_{ij}$ を用いて次のように計算されます。

\[d^2 = \sum_{i=1}^{3}\sum_{j=1}^{3} \Delta p_i g_{ij} \Delta p_j\]

より具体的には、$\Delta p = (\Delta x, \Delta y, \Delta z)$ とすると、 $$d^2 = \Delta x^2 g_{11} + \Delta y^2 g_{22} + \Delta z^2 g_{33} + 2\Delta x \Delta y g_{12} + 2\Delta y \Delta z g_{23} + 2\Delta z \Delta x g_{31}$$ または、行列形式で表すと、 $$d^2 = (\vec{p_1} + \vec{n} - \vec{p_0})^T \cdot g \cdot (\vec{p_1} + \vec{n} - \vec{p_0})$$

この距離計算は tkcrystalbase.distance 関数によって実行されます。プログラムは、指定された距離範囲 $r_{min}$ から $r_{max}$ の間のすべての距離を収集し、ソートして出力します。周期境界条件を適用するために、並進ベクトル $(ix, iy, iz)$ は、格子定数と $r_{max}$ に基づいて決定された範囲 (nxmax, nymax, nzmax) でループされます。
逆格子情報 (Reciprocal Lattice Information) プログラムは、実空間の格子ベクトルから逆格子ベクトルを計算し、さらに逆格子パラメーター、逆格子メトリックテンソル、逆格子単位胞体積も計算します。これらはそれぞれ tkcrystalbase.cal_reciprocal_lattice_vectors, tkcrystalbase.cal_reciprocal_lattice_parameters, tkcrystalbase.cal_metrics, tkcrystalbase.cal_volume 関数によって計算され、標準出力に表示されます。

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

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

numpy: 高性能な数値計算（配列操作、線形代数など）を可能にするPythonライブラリです。
matplotlib: グラフ描画ライブラリです。特に mpl_toolkits.mplot3d は3Dプロット機能を提供しますが、このプログラムの現在のバージョンではプロットは実行されていません。
tkcrystalbase: 結晶学計算用のカスタムライブラリです。このファイル（tkcrystalbase.py）は、crystal_distance.py と同じディレクトリ、またはPythonのパスが通っている場所に存在する必要があります。

これらのライブラリは、pip コマンドを使用してインストールできます。

pip install numpy matplotlib

tkcrystalbase モジュールについては、通常、crystal_distance.py と同じディレクトリに tkcrystalbase.py ファイルとして配置する必要があります。

必要な入力ファイル

このプログラムは、外部の入力ファイルを必要としません。すべての設定（格子パラメーター、サイト情報、距離範囲）は、プログラムのソースコード内に直接ハードコードされています。

lattice_parameters:
- 形式: Pythonのリスト [a, b, c, alpha, beta, gamma]
- 内容: 単位胞の格子定数 $(a, b, c)$ はオングストローム (Å) 単位、軸角 $(\alpha, \beta, \gamma)$ は度数 (degree) 単位。
- 例: [5.62, 5.62, 5.62, 90.0, 90.0, 90.0] （塩化ナトリウムの立方晶を示す例）
sites:
- 形式: Pythonのリストのリスト。各内部リストは1つの原子サイト情報。
- 内容: ['原子名', 'サイトラベル', 原子番号, 原子量, 電荷, 半径, '色', numpy.array([x, y, z])]
  - x, y, z: 単位胞内の分数座標。
- 例: ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])]
rmin, rmax:
- 形式: 浮動小数点数
- 内容: 計算対象とする原子間距離の最小値 ($r_{min}$) および最大値 ($r_{max}$)。単位はオングストローム (Å)。rmin は同一サイトの判定にも使用されます。

これらの変数を変更するには、直接 crystal_distance.py のソースコードを編集する必要があります。

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

このプログラムは、外部ファイルを出力として生成しません。すべての計算結果は、標準出力 (stdout) に直接表示されます。

出力される情報には以下のものが含まれます。

格子パラメーター: 設定された格子定数と軸角。
格子ベクトル: デカルト座標系における単位胞の格子ベクトル ($\vec{a}, \vec{b}, \vec{c}$)。
メトリックテンソル: 結晶格子に対応するメトリックテンソル。
単位胞体積: 単位胞の体積。
逆格子パラメーター: 逆格子空間における格子定数と軸角。
逆格子ベクトル: 逆格子空間における格子ベクトル ($\vec{a}^*, \vec{b}^*, \vec{c}^*$)。
逆格子メトリックテンソル: 逆格子空間のメトリックテンソル。
逆格子単位胞体積: 逆格子単位胞の体積。
探索範囲: 周期境界条件で探索される単位胞の最大範囲 (nmax)。
原子間距離リスト: $r_{min}$ から $r_{max}$ の範囲で計算されたすべての原子間距離のリスト。各エントリは、距離、開始サイト、終了サイト、および相対する単位胞の並進ベクトル $(ix, iy, iz)$ を含みます。このリストは距離でソートされて表示されます。

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

このプログラムは、コマンドライン引数を受け付けません。実行はPythonインタープリタを使用して、スクリプトファイルを直接指定するだけです。

python crystal_distance.py

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

デフォルトのハードコードされた設定（立方晶の格子定数と塩化ナトリウム型構造のサイト）でプログラムを実行します。

プログラムの準備 crystal_distance.py と、このプログラムが依存する tkcrystalbase.py が同じディレクトリにあることを確認してください。（注: tkcrystalbase.py は提供されていないため、ここでは存在すると仮定します。）
実行コマンド
```
python crystal_distance.py
```

実行結果の説明 プログラムは以下のような情報を標準出力に出力します。まず、設定された格子パラメーター、それから計算された格子ベクトル、メトリックテンソル、単位胞体積が表示されます。次に、逆格子に関する同様の情報が表示されます。その後、周期境界条件で考慮される単位胞の最大探索範囲 (nmax) が示されます。これは、rmax の値と格子定数に基づいて決定されます。最後に、計算されたすべての原子間距離が昇順でリストアップされます。各行は「開始サイトのラベル (分数座標) - 終了サイトのラベル (分数座標) + (並進ベクトルix, iy, iz): 距離」の形式で表示されます。

例えば、以下のような出力が期待されます。

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

Unit cell volume:     177.6 A^3
Reciprocal lattice parameters: [1.1189, 1.1189, 1.1189, 90.0, 90.0, 90.0]
Reciprocal lattice vectors:
  Rax: (  1.1189,        0,        0) A^-1
  Ray: (       0,   1.1189,        0) A^-1
  Raz: (       0,        0,   1.1189) A^-1
Reciprocal lattice metric tensor:
  Rgij: (  1.2519,        0,        0) A^-1
        (       0,   1.2519,        0) A^-1
        (       0,        0,   1.2519) A^-1
Reciprocal unit cell volume:   0.00563 A^-3

nmax: 1 1 1

Interatomic distances:
  Na1 (   0.0000,    0.0000,    0.0000) - Cl1 (   0.5000,    0.0000,    0.0000) + ( 0,  0,  0): dis =     2.8100 A
  Na1 (   0.0000,    0.0000,    0.0000) - Cl3 (   0.0000,    0.0000,    0.5000) + ( 0,  0,  0): dis =     2.8100 A
  Na1 (   0.0000,    0.0000,    0.0000) - Cl4 (   0.0000,    0.5000,    0.0000) + ( 0,  0,  0): dis =     2.8100 A
  Cl1 (   0.5000,    0.0000,    0.0000) - Na1 (   0.0000,    0.0000,    0.0000) + ( 0,  0,  0): dis =     2.8100 A
  Cl3 (   0.0000,    0.0000,    0.5000) - Na1 (   0.0000,    0.0000,    0.0000) + ( 0,  0,  0): dis =     2.8100 A
  Cl4 (   0.0000,    0.5000,    0.0000) - Na1 (   0.0000,    0.0000,    0.0000) + ( 0,  0,  0): dis =     2.8100 A
  Na1 (   0.0000,    0.0000,    0.0000) - Na2 (   0.0000,    0.5000,    0.5000) + ( 0,  0,  0): dis =     3.9740 A
  Na1 (   0.0000,    0.0000,    0.0000) - Na3 (   0.5000,    0.0000,    0.5000) + ( 0,  0,  0): dis =     3.9740 A
  Na1 (   0.0000,    0.0000,    0.0000) - Na4 (   0.5000,    0.5000,    0.0000) + ( 0,  0,  0): dis =     3.9740 A
  ... (以下、設定された rmax までの距離が出力されます)

上記の出力例では、Na1 サイトの原子と、現在の単位胞内のCl1 サイトの原子との距離が 2.8100 A であることが示されています。+ ( 0, 0, 0) は、両方の原子が同じ単位胞に属していることを意味します。もし + ( 1, 0, 0) のような値が表示されれば、それは Cl1 サイトの原子がX方向に1単位胞分シフトした隣の単位胞に位置する原子との距離を示します。