crystal_distance.py 技術ドキュメント
プログラムの動作
crystal_distance.py は、結晶構造の格子パラメーターと原子サイト情報に基づいて、原子間の距離を計算し表示するPythonプログラムです。
主な機能:
格子情報の計算と表示: 与えられた格子パラメーター（格子定数 \(a, b, c\) と角度 \(\alpha, \beta, \gamma\)）から、実空間および逆空間の基本格子ベクトル、メトリックテンソル、および単位胞体積を計算し、標準出力に表示します。
原子間距離の計算: プログラム内に定義された原子サイト情報と周期境界条件を考慮し、指定された最小距離 (
rmin) と最大距離 (rmax) の範囲内にある全ての原子ペア間の距離を計算します。結果のソートと表示: 計算された原子間距離のリストは、距離の昇順、続いて原子サイトのラベル、周期境界条件によるシフト量の順にソートされ、詳細な情報とともに標準出力に表示されます。
解決する課題:
結晶構造において、特定の原子とその周囲に存在する他の原子との間の距離を正確に特定することは、材料科学や物性物理学における構造解析や物性予測の基礎となります。このプログラムは、単位胞内に定義された原子だけでなく、周期的に繰り返される隣接単位胞内の原子との距離も自動的に考慮することで、この課題を解決します。
原理
このプログラムは、結晶学における基本的な幾何学的関係と距離計算の原則に基づいています。
1. 格子ベクトルとメトリックテンソル
結晶格子は、3つの基本格子ベクトル \(\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3\) によって定義されます。これらのベクトルは、格子定数 \(a, b, c\) とその間の角度 \(\alpha, \beta, \gamma\) から計算されます。
実空間におけるメトリックテンソル \(\mathbf{G}\) (コード内の gij) は、格子ベクトルの内積によって定義されます。
つまり、\(g_{ij} = \mathbf{a}_i \cdot \mathbf{a}_j\) です。
2. 原子間距離の計算
分数座標で表現された2つの原子 \(\mathbf{p}_0 = (x_0, y_0, z_0)\) と \(\mathbf{p}_1 = (x_1, y_1, z_1)\) の間の距離 \(d\) は、分数座標の差 \(\Delta \mathbf{p} = (x_1-x_0, y_1-y_0, z_1-z_0)\) とメトリックテンソル \(\mathbf{G}\) を用いて以下の式で計算されます。
具体的には、
$\(
d^2 = g_{11}(\Delta x)^2 + g_{22}(\Delta y)^2 + g_{33}(\Delta z)^2 + 2g_{12}(\Delta x)(\Delta y) + 2g_{23}(\Delta y)(\Delta z) + 2g_{31}(\Delta z)(\Delta x)
\)$
この計算は、外部モジュール tkcrystalbase.py の distance 関数によって実行されます。
3. 周期境界条件の考慮
結晶は無限に繰り返される構造を持つため、単一の単位胞内に存在する原子だけでなく、隣接する単位胞内の原子との距離も考慮する必要があります。これは、原子 \(\mathbf{p}_1\) の分数座標に整数シフト \((ix, iy, iz)\) を加えることで実現されます。すなわち、\(\mathbf{p}_1' = \mathbf{p}_1 + (ix, iy, iz)\) とし、この \(\mathbf{p}_1'\) を用いて \(\mathbf{p}_0\) との距離を計算します。
整数シフトの範囲 \((ix, iy, iz)\) は、計算対象となる最大距離 rmax と格子パラメーターに基づいて動的に決定されます。これにより、rmax の範囲内にある全ての原子ペアが探索されるようになっています。
4. 逆格子
プログラムは、実空間の格子情報と同様に、逆格子ベクトル、逆格子パラメーター、逆格子メトリックテンソル、および逆格子単位胞体積も計算し表示します。逆格子は回折現象の解析などで重要となります。これらの計算も tkcrystalbase.py の関数群によって行われます。
必要な非標準ライブラリとインストール方法
このプログラムの実行には、以下の非標準ライブラリが必要です。
numpy: 数値計算を効率的に行うためのライブラリ。ベクトルや行列の操作に使用されます。
matplotlib: グラフ描画ライブラリ。プログラム中では
mpl_toolkits.mplot3dとmatplotlib.pyplotがインポートされていますが、明示的なプロット機能は使用されていません。tkcrystalbase: このプログラムと共に提供される、結晶学計算用のカスタムPythonモジュール。
インストール方法:
numpy および matplotlib は、pip を使用してインストールできます。
pip install numpy matplotlib
tkcrystalbase.py は、crystal_distance.py と同じディレクトリに配置するか、Pythonのモジュール検索パスが通っているディレクトリに配置する必要があります。
必要な入力ファイル
このプログラムは、外部の入力ファイルを読み込みません。すべての必要なデータは、プログラムのソースコード内に直接記述されています。
lattice_parameters: 格子定数 \(a, b, c\) (Å) および格子角 \(\alpha, \beta, \gamma\) (度) のリスト。デフォルトでは立方晶のパラメーターが設定されています。 例:[ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]sites: 単位胞内の原子サイト情報を含むリストのリスト。各原子サイトは以下の情報で構成されます。原子名 (str)
サイトラベル (str)
原子番号 (int)
原子質量 (float)
電荷 (float)
原子半径 (float)
色 (str)
単位胞内分数座標 (numpy.array) 例:
['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])]
rmin: 原子間距離計算の最小閾値 (Å)。これより短い距離は無視されます。また、同一サイトの判定にも使用されます。rmax: 原子間距離計算の最大閾値 (Å)。これより長い距離は無視されます。
これらの変数は、プログラムの冒頭部分で直接編集することで変更できます。
生成される出力ファイル
このプログラムは、いかなるファイルも生成しません。すべての計算結果は、標準出力 (stdout) に直接表示されます。
コマンドラインでの使用例 (Usage)
このプログラムは、引数を取らずに直接実行されます。
python crystal_distance.py
コマンドラインでの具体的な使用例
crystal_distance.py をデフォルト設定のまま実行した場合の例を以下に示します。
python crystal_distance.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.584, 0, 0) A
( 0, 31.584, 0) A
( 0, 0, 31.584) A
Volume: 177.5 A^3
Unit cell volume: 177.5 A^3
Reciprocal lattice parameters: [1.1187, 1.1187, 1.1187, 90.0, 90.0, 90.0]
Reciprocal lattice vectors:
Rax: ( 1.1187, 0, 0) A^-1
Ray: ( 0, 1.1187, 0) A^-1
Raz: ( 0, 0, 1.1187) A^-1
Reciprocal lattice metric tensor:
Rgij: ( 1.2514, 0, 0) A^-1
( 0, 1.2514, 0) A^-1
( 0, 0, 1.2514) A^-1
Reciprocal unit cell volume: 0.005634 A^-3
nmax: 0 0 0
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
Na1 ( 0.0000, 0.0000, 0.0000) - Cl1 ( 0.5000, 0.0000, 0.0000) + (-1, 0, 0): dis = 2.8100 A
Na1 ( 0.0000, 0.0000, 0.0000) - Cl3 ( 0.0000, 0.0000, 0.5000) + ( 0, 0, -1): dis = 2.8100 A
Na1 ( 0.0000, 0.0000, 0.0000) - Cl4 ( 0.0000, 0.5000, 0.0000) + ( 0, -1, 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
Na1 ( 0.0000, 0.0000, 0.0000) - Na2 ( 0.0000, 0.5000, 0.5000) + ( 0, -1, -1): dis = 3.9740 A
Na1 ( 0.0000, 0.0000, 0.0000) - Na3 ( 0.5000, 0.0000, 0.5000) + (-1, 0, -1): dis = 3.9740 A
Na1 ( 0.0000, 0.0000, 0.0000) - Na4 ( 0.5000, 0.5000, 0.0000) + (-1, -1, 0): dis = 3.9740 A
Na1 ( 0.0000, 0.0000, 0.0000) - Cl2 ( 0.5000, 0.5000, 0.5000) + ( 0, 0, 0): dis = 4.8606 A
Cl1 ( 0.5000, 0.0000, 0.0000) - Na2 ( 0.0000, 0.5000, 0.5000) + ( 0, 0, 0): dis = 4.8606 A
...
実行結果の説明:
格子パラメーターと格子ベクトルの情報: プログラムの冒頭で定義された格子パラメーターと、それから計算された実空間の格子ベクトル、メトリックテンソル、および単位胞体積が表示されます。デフォルトでは立方晶のパラメーターが設定されており、それに伴い格子ベクトルやメトリックテンソルは対角成分のみを持ちます。
逆格子情報の表示: 同様に、逆格子パラメーター、逆格子ベクトル、逆格子メトリックテンソル、および逆格子単位胞体積が計算され表示されます。
探索範囲 (
nmax): 周期境界条件を適用する際の単位胞の繰り返し数 (nxmax,nymax,nzmax) が表示されます。ここではrmaxが 4.5 Å でlattice_parameters[0]が 5.62 Å のため、int(4.5 / 5.62) + 1が0 + 1 = 1となりnmax: 0 0 0と表示されています。これは、\(ix, iy, iz\) のループが-0から+0まで、つまり0のみで実行されることを意味します。そのため、実質的に単位胞内の原子と、その原子の近傍に位置する隣接単位胞内の原子との距離が計算されます。原子間距離のリスト: 最後に、「Interatomic distances:」の見出しの下に、計算された全ての原子間距離がソートされて表示されます。各行は以下の情報を含みます。
サイト0のラベルと分数座標
サイト1のラベルと分数座標
周期境界条件によるシフト量
(ix, iy, iz): これは、サイト1の原子が基準となる単位胞からどの単位胞に位置するかを示します。(0, 0, 0) は同じ単位胞内、例えば (-1, 0, 0) は \(x\) 方向に1つ隣の単位胞に位置するサイト1の原子を指します。計算された距離
dis(Å)
この出力により、結晶内の原子配置と、各原子間の結合長や相互作用距離を詳細に把握することができます。