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

プログラムの動作

interatomic_distance.py は、指定された格子定数と原子サイト情報に基づいて、結晶構造における原子間距離を計算するPythonスクリプトです。このプログラムの主な目的は、単位胞内の原子だけでなく、周期境界条件を考慮して指定された距離範囲内の全ての原子ペア間の距離を列挙し、距離順にソートして出力することです。

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

結晶パラメータの計算: 実空間および逆空間の格子定数から、格子ベクトル、計量テンソル、単位胞体積を計算します。
周期境界条件の適用: 指定された最大距離 (rmax) に基づいて、周囲の単位胞を考慮し、原子間の距離を計算します。
原子間距離の計算: 結晶学における計量テンソルを用いて、分数座標で指定された原子間の距離を正確に計算します。
結果のソートと表示: 計算された全ての原子間距離を昇順にソートし、どの原子ペアがどの並進ベクトルによって関連付けられているかを含め、標準出力に詳細情報を表示します。

このプログラムは、特定の結晶構造（デフォルトではNaCl構造）における原子間相互作用の分析や、構造内の最短距離結合の特定などの課題を解決するのに役立ちます。

原理

このプログラムは、結晶学の基本的な概念と計算手法に基づいています。

格子ベクトルと計量テンソル: 結晶構造は3つの基本格子ベクトル \(\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3\) によって定義されます。これらのベクトルは、格子定数 \((a, b, c, \alpha, \beta, \gamma)\) から tkcrystalbase モジュールの cal_lattice_vectors 関数によって計算されます。計量テンソル \(\mathbf{G}\) は、格子ベクトル間の内積から構成される対称行列であり、結晶空間における距離の計算に不可欠です。その要素 \(g_{ij}\) は以下のように定義されます。

\[g_{ij} = \mathbf{a}_i \cdot \mathbf{a}_j\]

行列形式では、

\[\begin{split} \mathbf{G} = \begin{pmatrix} \mathbf{a}_1 \cdot \mathbf{a}_1 & \mathbf{a}_1 \cdot \mathbf{a}_2 & \mathbf{a}_1 \cdot \mathbf{a}_3 \\ \mathbf{a}_2 \cdot \mathbf{a}_1 & \mathbf{a}_2 \cdot \mathbf{a}_2 & \mathbf{a}_2 \cdot \mathbf{a}_3 \\ \mathbf{a}_3 \cdot \mathbf{a}_1 & \mathbf{a}_3 \cdot \mathbf{a}_2 & \mathbf{a}_3 \cdot \mathbf{a}_3 \end{pmatrix} \end{split}\]

これは cal_metrics 関数によって計算されます。
単位胞の体積: 単位胞の体積 \(V\) は、3つの基本格子ベクトルのスカラー三重積の絶対値として与えられます。

\[V = |\mathbf{a}_1 \cdot (\mathbf{a}_2 \times \mathbf{a}_3)|\]

これは cal_volume 関数によって計算されます。
原子間距離の計算: 結晶学において、原子の位置は分数座標 \(\mathbf{r} = (x, y, z)\) で表されます。2つの原子の分数座標がそれぞれ \(\mathbf{r}_0\) と \(\mathbf{r}_1\) である場合、これらの間の距離 \(d\) は計量テンソル \(\mathbf{G}\) を用いて次のように計算されます。

\[\Delta\mathbf{r} = \mathbf{r}_0 - \mathbf{r}_1\]

\[d = \sqrt{\Delta\mathbf{r}^T \mathbf{G} \Delta\mathbf{r}}\]

ここで、\(\Delta\mathbf{r}^T\) はベクトル \(\Delta\mathbf{r}\) の転置を表します。プログラムでは、distance 関数がこの計算を実行します。
周期境界条件の適用: 結晶は無限に繰り返される構造であるため、特定の原子から見た隣接原子は、現在の単位胞だけでなく、周囲の並進した単位胞にも存在します。プログラムは、最大距離 rmax に基づいて、考慮すべき並進ベクトルの範囲 \((ix, iy, iz)\) を決定します。これにより、原子 site0 と、並進ベクトル \(\mathbf{T} = ix\mathbf{a}_1 + iy\mathbf{a}_2 + iz\mathbf{a}_3\) でシフトされた位置にある site1 との間の距離が計算されます。分数座標では、これは pos1 + np.array([ix, iy, iz]) として表現されます。

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

このプログラムは以下の非標準ライブラリを使用します。

NumPy: 数値計算、特にベクトルや行列の操作に使用されます。
- インストール方法:
```
pip install numpy
```
Matplotlib: グラフ描画ライブラリですが、このプログラムではインポートされているものの、直接的な描画処理には使用されていません。
- インストール方法:
```
pip install matplotlib
```
tkcrystalbase: 結晶学的な計算（格子ベクトル、計量テンソル、距離計算など）を提供するカスタムモジュールです。このプログラムの実行には、interatomic_distance.py と同じディレクトリに tkcrystalbase.py ファイルが存在する必要があります。

必要な入力ファイル

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

ただし、結晶学的な計算を行うための補助モジュールである tkcrystalbase.py が、interatomic_distance.py と同じディレクトリに配置されている必要があります。

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

このプログラムは、いかなるファイルも生成しません。全ての計算結果は標準出力（コンソール）に表示されます。出力には、格子定数、格子ベクトル、計量テンソル、単位胞体積、およびソートされた原子間距離のリストが含まれます。

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

interatomic_distance.py はコマンドライン引数を取らずに実行されます。

python interatomic_distance.py

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

プログラムを実行すると、以下のような情報が標準出力に表示されます。（出力は、プログラム内で設定されているデフォルトのNaCl構造の格子定数 [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0] とサイト情報に基づいています。）

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: [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.005639 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
  Cl1 (  0.5000,   0.0000,   0.0000) - Na2 (  0.0000,   0.5000,   0.5000) + ( 1, -1, -1): dis =     2.8100 A
  Cl1 (  0.5000,   0.0000,   0.0000) - Na3 (  0.5000,   0.0000,   0.5000) + ( 0,  0, -1): dis =     2.8100 A
  Cl1 (  0.5000,   0.0000,   0.0000) - Na4 (  0.5000,   0.5000,   0.0000) + ( 0, -1,  0): dis =     2.8100 A
  Na2 (  0.0000,   0.5000,   0.5000) - Cl1 (  0.5000,   0.0000,   0.0000) + ( 0,  1,  1): dis =     2.8100 A
  Na2 (  0.0000,   0.5000,   0.5000) - Cl2 (  0.5000,   0.5000,   0.5000) + ( 0,  0,  0): dis =     2.8100 A
  Na2 (  0.0000,   0.5000,   0.5000) - Cl3 (  0.0000,   0.0000,   0.5000) + ( 0,  1,  0): dis =     2.8100 A
  Na2 (  0.0000,   0.5000,   0.5000) - Cl4 (  0.0000,   0.5000,   0.0000) + ( 0,  0,  1): dis =     2.8100 A
  Na3 (  0.5000,   0.0000,   0.5000) - Cl1 (  0.5000,   0.0000,   0.0000) + ( 0,  0,  0): dis =     2.8100 A
  Na3 (  0.5000,   0.0000,   0.5000) - Cl2 (  0.5000,   0.5000,   0.5000) + ( 0, -1,  0): dis =     2.8100 A
  Na3 (  0.5000,   0.0000,   0.5000) - Cl3 (  0.0000,   0.0000,   0.5000) + ( 1,  0,  0): dis =     2.8100 A
  Na3 (  0.5000,   0.0000,   0.5000) - Cl4 (  0.0000,   0.5000,   0.0000) + ( 1, -1,  1): dis =     2.8100 A
  Na4 (  0.5000,   0.5000,   0.0000) - Cl1 (  0.5000,   0.0000,   0.0000) + ( 0,  1,  0): dis =     2.8100 A
  Na4 (  0.5000,   0.5000,   0.0000) - Cl2 (  0.5000,   0.5000,   0.5000) + ( 0,  0, -1): dis =     2.8100 A
  Na4 (  0.5000,   0.5000,   0.0000) - Cl3 (  0.0000,   0.0000,   0.5000) + ( 1,  1, -1): dis =     2.8100 A
  Na4 (  0.5000,   0.5000,   0.0000) - Cl4 (  0.0000,   0.5000,   0.0000) + ( 1,  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
  Cl1 (  0.5000,   0.0000,   0.0000) - Cl2 (  0.5000,   0.5000,   0.5000) + ( 0,  0,  0): dis =     3.9740 A
  Cl1 (  0.5000,   0.0000,   0.0000) - Cl3 (  0.0000,   0.0000,   0.5000) + ( 1,  0,  0): dis =     3.9740 A
  Cl1 (  0.5000,   0.0000,   0.0000) - Cl4 (  0.0000,   0.5000,   0.0000) + ( 1,  0,  0): dis =     3.9740 A
  Na2 (  0.0000,   0.5000,   0.5000) - Na3 (  0.5000,   0.0000,   0.5000) + ( 0,  0,  0): dis =     3.9740 A
  Na2 (  0.0000,   0.5000,   0.5000) - Na4 (  0.5000,   0.5000,   0.0000) + ( 0,  0,  0): dis =     3.9740 A
  Na3 (  0.5000,   0.0000,   0.5000) - Na4 (  0.5000,   0.5000,   0.0000) + ( 0,  0,  0): dis =     3.9740 A
  Cl2 (  0.5000,   0.5000,   0.5000) - Cl3 (  0.0000,   0.0000,   0.5000) + ( 1,  1,  0): dis =     3.9740 A
  Cl2 (  0.5000,   0.5000,   0.5000) - Cl4 (  0.0000,   0.5000,   0.0000) + ( 1,  0,  1): dis =     3.9740 A
  Cl3 (  0.0000,   0.0000,   0.5000) - Cl4 (  0.0000,   0.5000,   0.0000) + ( 0,  0,  0): dis =     3.9740 A