lattice_plane.py 技術ドキュメント

プログラムの動作

lattice_plane.py は、結晶学における単位格子の指定されたミラー指数 (hkl) 面を3次元で可視化するPythonプログラムです。結晶の格子定数と軸角が与えられた非斜交座標系において、単位格子の辺と(hkl)面の交点を計算し、その交点から形成される多角形として面を描画します。

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

  • 単位格子の定義: ユーザーが指定した格子定数 a, b, c および格子角 alpha, beta, gamma に基づいて、斜交座標系の単位格子をデカルト座標系に変換し、その形状を定義します。

  • ミラー指数 (hkl) 面の計算: コマンドライン引数またはプログラム内のデフォルト設定で与えられたミラー指数 h, k, l を用いて、単位格子を横切る(hkl)面の方程式 \(hx + ky + lz = 1\) (分数座標系) を定義します。

  • 平面と単位格子の交点計算: 単位格子の12本の辺と(hkl)面との交点を計算し、それらの交点からなる多角形を生成します。

  • 3D描画: 計算された単位格子と(hkl)面を matplotlib を使用して3D空間に描画します。描画には格子軸のラベル 'a', 'b', 'c' およびミラーインターセプト \(1/h, 1/k, 1/l\) も含められます。

  • 出力: 生成された3DプロットをPNG画像ファイルとして保存し、同時に画面に表示します。

このプログラムは、結晶構造の理解や教育、あるいは特定の結晶面を視覚的に確認したい場合に有用です。

原理

このプログラムは、結晶学における単位格子とミラー指数面の幾何学的な関係をデカルト座標系で表現する原理に基づいています。

1. 分数座標系からデカルト座標系への変換

結晶学では、単位格子の各辺に沿った軸を基底ベクトル \(\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3\) とし、任意の点 \(\mathbf{r}\) を分数座標 \((x, y, z)\) で次のように表します。

\[\mathbf{r} = x\mathbf{a}_1 + y\mathbf{a}_2 + z\mathbf{a}_3\]

ここで、\(x, y, z\) は各軸に沿った分数的な位置を示します。 このプログラムでは、格子定数 \(a, b, c\) と格子角 \(\alpha, \beta, \gamma\) (度) が与えられています。これらの情報を用いて、デカルト座標系における基底ベクトルを計算します。 慣習的に、\(\mathbf{a}_1\) はX軸上、\(\mathbf{a}_2\) はXY平面上にあるように設定されます。

  • \(\mathbf{a}_1 = (a, 0, 0)\)

  • \(\mathbf{a}_2 = (b \cos \gamma, b \sin \gamma, 0)\)

  • \(\mathbf{a}_3 = (c \cos \beta, c \frac{\cos \alpha - \cos \beta \cos \gamma}{\sin \gamma}, c \sqrt{1 - \cos^2 \beta - \left(\frac{\cos \alpha - \cos \beta \cos \gamma}{\sin \gamma}\right)^2})\)

これらのベクトルを列として持つ変換行列 \(M\) を作成することで、分数座標 \((x, y, z)^T\) をデカルト座標 \((X, Y, Z)^T\) に変換できます。

\[\begin{split} M = \begin{pmatrix} a & b \cos \gamma & c \cos \beta \\ 0 & b \sin \gamma & c \frac{\cos \alpha - \cos \beta \cos \gamma}{\sin \gamma} \\ 0 & 0 & c \sqrt{1 - \cos^2 \beta - \left(\frac{\cos \alpha - \cos \beta \cos \gamma}{\sin \gamma}\right)^2} \end{pmatrix} \end{split}\]

そして、デカルト座標は \(\mathbf{r}_{cart} = M \cdot (x, y, z)^T\) で求められます。

2. ミラー指数 (hkl) 面の方程式

ミラー指数 \((h, k, l)\) で定義される結晶面の方程式は、分数座標系において次のように表されます。

\[hx + ky + lz = 1\]

この方程式は、各軸を \(1/h, 1/k, 1/l\) で切り取る平面を示します（ただし、\(h, k, l\) のいずれかが0の場合、その軸に平行になります）。

3. 平面と単位格子の辺の交点計算

単位格子は、分数座標で頂点 \((0,0,0)\) から \((1,1,1)\) までの8つの頂点と、それらをつなぐ12本の辺で構成されます。 各辺は2つの頂点 \(\mathbf{v}_0\) と \(\mathbf{v}_1\) を持つ線分として表現できます。線分上の任意の点 \(\mathbf{p}(t)\) は、パラメータ \(t\) (\(0 \le t \le 1\)) を用いて次のように表されます。

\[\mathbf{p}(t) = \mathbf{v}_0 + t(\mathbf{v}_1 - \mathbf{v}_0)\]

ここで \(\mathbf{d} = \mathbf{v}_1 - \mathbf{v}_0\) とおきます。この点 \(\mathbf{p}(t) = (x(t), y(t), z(t))\) が(hkl)面上にあるとすると、平面の方程式に代入して \(t\) を解きます。

\[h(v_{0x} + t d_x) + k(v_{0y} + t d_y) + l(v_{0z} + t d_z) = 1\]
\[(h v_{0x} + k v_{0y} + l v_{0z}) + t (h d_x + k d_y + l d_z) = 1\]
\[t = \frac{1 - (h v_{0x} + k v_{0y} + l v_{0z})}{h d_x + k d_y + l d_z}\]

\(t\) の値が \(0 \le t \le 1\) の範囲内であれば、その交点は線分上に存在します。計算されたすべての交点は、(hkl)面を構成する多角形の頂点となります。

4. 多角形頂点のソート

計算された交点群は、ランダムな順序で得られます。これらを正しく多角形として描画するためには、平面上で連続する頂点順序に並べ替える必要があります。 このプログラムでは、以下の手順でソートを行います。

  1. 交点群の重心 cen を計算します。

  2. 平面の法線ベクトル \(\mathbf{n} = (h, k, l)\) を正規化します。

  3. 法線ベクトル \(\mathbf{n}\) に直交する2つのベクトル \(\mathbf{u}\) と \(\mathbf{v}\) を生成し、平面内の局所的な座標系を定義します。

    • \(\mathbf{u}\) は \(\mathbf{n}\) とX軸（またはY軸）の外積から生成されます。

    • \(\mathbf{v}\) は \(\mathbf{n}\) と \(\mathbf{u}\) の外積から生成されます。

  4. 各交点 \(\mathbf{p}\) について、重心 cen からの相対ベクトル \(\mathbf{p} - \text{cen}\) を計算し、\(\mathbf{u}\) と \(\mathbf{v}\) を基底とする平面上での角度 \(\text{arctan2}(\text{dot}(\mathbf{p}-\text{cen}, \mathbf{v}), \text{dot}(\mathbf{p}-\text{cen}, \mathbf{u}))\) を求めます。

  5. この角度に基づいて交点群をソートすることで、多角形の辺を連続的に描画するための正しい頂点順序が得られます。

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

このプログラムの実行には、以下の非標準ライブラリが必要です。

  • numpy: 数値計算を効率的に行うためのライブラリ。

  • matplotlib: 2D/3Dグラフ描画のためのライブラリ。mpl_toolkits.mplot3d は3D描画機能を提供します。

これらのライブラリは、Pythonのパッケージマネージャー pip を使ってインストールできます。コマンドプロンプトやターミナルで以下のコマンドを実行してください。

pip install numpy matplotlib

必要な入力ファイル

lattice_plane.py は、実行するために特定の入力ファイルを必要としません。 ミラー指数 h, k, l は、プログラム内のデフォルト値（h=1, k=2, l=3）を使用するか、コマンドライン引数として直接渡すことができます。

生成される出力ファイル

プログラムを実行すると、以下の画像ファイルが生成されます。

  • ファイル名: lattice_planeHKL.png

    • H, K, L は、実行時に指定されたミラー指数（h, k, l）の実際の値に置き換わります。

    • 例: h=1, k=2, l=3 の場合は lattice_plane123.png が生成されます。

  • 内容: 指定された格子定数と軸角を持つ単位格子、およびその単位格子を横切る(hkl)面が3Dで描画された画像。格子軸のラベル 'a', 'b', 'c' とミラーインターセプト \(1/h, 1/k, 1/l\) も表示されます。

  • 形式: PNG画像 (300 dpi)

  • その他: 画像ファイル保存後、MatplotlibのGUIウィンドウに描画結果がポップアップ表示されます。

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

lattice_plane.py は、ミラー指数 h, k, l をコマンドライン引数として受け取ります。引数が与えられない場合、プログラム内のデフォルト値が使用されます。

python lattice_plane.py [h] [k] [l]

  • [h]: ミラー指数 h を指定します（整数）。省略可能。

  • [k]: ミラー指数 k を指定します（整数）。省略可能。

  • [l]: ミラー指数 l を指定します（整数）。省略可能。

引数はスペースで区切って指定します。指定しない引数はデフォルト値が使われますが、指定する場合は左から順に h, k, l を指定する必要があります。例えば、l だけを指定することはできません。

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

例1: デフォルトのミラー指数を使用する

ミラー指数をコマンドラインで指定しない場合、プログラム内のデフォルト値 h=1, k=2, l=3 が使用されます。

python lattice_plane.py

実行結果: lattice_plane123.png というファイルが生成され、単位格子と(123)面が描画されます。Matplotlibウィンドウにも同様の図が表示されます。 図には、格子軸 a, b, c とミラーインターセプト \(1/1, 1/2, 1/3\) の位置が示されます。

例2: (110)面を描画する

h=1, k=1, l=0 のミラー指数を指定して実行します。

python lattice_plane.py 1 1 0

実行結果: lattice_plane110.png というファイルが生成され、単位格子と(110)面が描画されます。Matplotlibウィンドウにも同様の図が表示されます。 図には、格子軸 a, b, c とミラーインターセプト \(1/1, 1/1\) の位置が示されます（\(l=0\) のため \(1/0\) は表示されず、c軸に平行な面として描画されます）。

例3: (211)面を描画する

h=2, k=1, l=1 のミラー指数を指定して実行します。

python lattice_plane.py 2 1 1

実行結果: lattice_plane211.png というファイルが生成され、単位格子と(211)面が描画されます。Matplotlibウィンドウにも同様の図が表示されます。 図には、格子軸 a, b, c とミラーインターセプト \(1/2, 1/1, 1/1\) の位置が示されます。