crystal_convert_cell.py 技術ドキュメント
プログラムの動作
crystal_convert_cell.py は、指定された結晶の単位格子を、異なる基底ベクトルを持つ新しい単位格子に変換し、その両方の単位格子と含まれる原子サイトを3次元で可視化するPythonプログラムです。
主な機能は以下の通りです。
結晶情報の表示: 初期状態の結晶（格子定数、格子ベクトル、メートルテンソル、体積、原子サイト）を詳細に表示します。
単位格子の変換: ユーザーが指定した変換モード（例: 面心立方格子(FCC)をプリミティブセルに、体心立方格子(BCC)をプリミティブセルに、菱面体晶を六方晶になど）に基づいて、格子ベクトルと原子サイトの座標を変換します。
変換情報の表示: 変換後の単位格子の詳細情報（格子定数、格子ベクトル、メートルテンソル、体積、原子サイト）を表示します。
3D可視化: Matplotlibの3Dプロット機能を使用して、元の単位格子と変換後の単位格子を同時に描画します。これにより、両者の関係や原子サイトの配置を直感的に理解できます。元の格子は黒い線で、変換後の格子は青い線と灰色の原子で表示されます。
結晶学教育・研究支援: 結晶学における単位格子の概念、基底ベクトルの変換、および並進対称性を持つ構造の異なる記述方法を理解するためのツールとして役立ちます。
原理
このプログラムは、結晶学の基本的な概念と線形代数に基づいています。
1. 格子ベクトルと格子パラメータ
結晶構造は、通常、6つの格子パラメータ (\(a, b, c, \alpha, \beta, \gamma\)) で記述されます。 これらは単位格子の3辺の長さと3つの角度を表します。 これらのパラメータから、実空間の基底ベクトル \(\vec{a}_1, \vec{a}_2, \vec{a}_3\) を計算できます。 例えば、直交座標系では以下のように設定できます。
プログラム内部では、tkcrystalbase モジュールの cal_lattice_vectors 関数がこの計算を行います。
2. メートルテンソル
メートルテンソル \(g_{ij}\) は、基底ベクトル間の内積によって定義され、結晶の幾何学的特性を表します。
$\(g_{ij} = \vec{a}_i \cdot \vec{a}_j\)$
プログラム内部では、tkcrystalbase モジュールの cal_metrics 関数がこの計算を行います。
3. 単位格子の体積
単位格子の体積 \(V\) は、基底ベクトルのスカラー三重積によって計算されます。
$\(V = |\vec{a}_1 \cdot (\vec{a}_2 \times \vec{a}_3)|\)$
プログラム内部では、tkcrystalbase モジュールの cal_volume 関数がこの計算を行います。
4. 格子ベクトルの変換
新しい単位格子の基底ベクトル \(\vec{a}'_1, \vec{a}'_2, \vec{a}'_3\) は、元の基底ベクトル \(\vec{a}_1, \vec{a}_2, \vec{a}_3\) と変換行列 \(T\) を用いて記述されます。 変換行列 \(T\) は \(3 \times 3\) の行列 \(T_{ij}\) であり、以下のように定義されます。
すなわち、\(\vec{a}'_i = \sum_{j=1}^3 T_{ij} \vec{a}_j\) です。
プログラム内部では、get_conversion_matrix 関数が T 行列を返し、convert_lattice_vectors 関数がこの変換を行います。
5. 原子サイト（分数座標）の変換
実空間のある点 \(\vec{r}\) が、元の単位格子で分数座標 \((x_1, x_2, x_3)\) を持つ場合、 $\(\vec{r} = x_1\vec{a}_1 + x_2\vec{a}_2 + x_3\vec{a}_3\)\( と書けます。これをベクトル形式で \)\vec{r} = \mathbf{x}^T \mathbf{A}\( とします（ここで \)\mathbf{x}\( は分数座標の列ベクトル、\)\mathbf{A}$ は行ベクトルとして基底ベクトルを並べたもの）。
同様に、変換後の単位格子での分数座標を \((x'_1, x'_2, x'_3)\) とすると、 $\(\vec{r} = x'_1\vec{a}'_1 + x'_2\vec{a}'_2 + x'_3\vec{a}'_3 = \mathbf{x}'^T \mathbf{A}'\)$ となります。
\(\mathbf{A}' = \mathbf{T} \mathbf{A}\) の関係を代入すると、 $\(\mathbf{x}^T \mathbf{A} = \mathbf{x}'^T (\mathbf{T} \mathbf{A}) = (\mathbf{x}'^T \mathbf{T}) \mathbf{A}\)\( したがって、\)\mathbf{x}^T = \mathbf{x}'^T \mathbf{T}\( が成り立ち、両辺を転置すると、 \)\(\mathbf{x} = \mathbf{T}^T \mathbf{x}'\)\( これを \)\mathbf{x}'\( について解くと、 \)\(\mathbf{x}' = (\mathbf{T}^T)^{-1} \mathbf{x}\)$ となります。
プログラムでは、tij_x2xc = np.linalg.inv(tij).transpose() が \((\mathbf{T}^T)^{-1}\) を計算しており、これは \((\mathbf{T}^{-1})^T\) と等価です。
その後、convert_fractional_coordinates_by_tij 関数がこの行列と元の分数座標を用いて新しい分数座標を計算します。
また、結晶学では分数座標を \([0, 1)\) の範囲に収めるため、reduce01 関数が剰余演算を用いてこの処理を行います。
6. 原子サイトの重複判定
原子サイトを追加する際には、rmin (最小距離) を基準に、既存のサイトと重複しないように判定が行われます。
これは、メートルテンソルを用いて2つの分数座標間の実空間距離を計算し、その距離が rmin より小さい場合は同一サイトとみなすことで実現されます。tkcrystalbase モジュールの add_site 関数がこの機能を提供します。
7. 3D可視化
Matplotlibライブラリを用いて、計算された格子ベクトルと原子サイトの座標に基づき、3D空間に単位格子と原子を描画します。
tkcrystalbase モジュールの draw_unitcell 関数が描画の中心的な役割を担います。
必要な非標準ライブラリとインストール方法
このプログラムの実行には、以下の非標準ライブラリが必要です。
NumPy: 数値計算（特に配列操作や線形代数）
Matplotlib: グラフ描画（特に3Dプロット）
tkcrystalbase: プログラムの実行ディレクトリに存在するカスタムモジュール。
NumPy と Matplotlib は pip コマンドでインストールできます。
pip install numpy matplotlib
tkcrystalbase.py は、このプログラム crystal_convert_cell.py と同じディレクトリに配置されている必要があります。このドキュメントの範囲外ですが、tkcrystalbase.py には cal_lattice_vectors, cal_metrics, cal_volume, convert_lattice_vectors, cal_lattice_parameters_from_lattice_vectors, get_conversion_matrix_from_tij, reduce01, add_site, draw_unitcell などの関数が含まれていると推測されます。
必要な入力ファイル
このプログラムはコマンドライン引数を受け取ります。
特定の入力ファイルは必要ありませんが、プログラム自体がcrystal_convert_cell.pyとして存在し、tkcrystalbase.pyが同じディレクトリにある必要があります。
コマンドライン引数として以下の情報を提供できます。
crystal_name: 変換元の結晶の種類を指定します。現在サポートされているのは'FCC','BCC','Hex','Rhomb'です。conversion_mode: 適用する単位格子変換のモードを指定します。現在サポートされているのは'FCCPrim','BCCPrim','HexOrtho','HexRhomb','RhombHex','ACenterPrim','BCenterPrim','CCenterPrim'です。kRatom: 描画される原子半径の表示倍率（係数）を指定する浮動小数点数です。
これらの引数が指定されない場合、プログラムはデフォルト値を使用します。 デフォルト値はコード内で以下のように設定されています。
crystal_name = 'FCC'conversion_mode = 'FCCPrim'kr = 1000.0
生成される出力ファイル
このプログラムはファイル出力を生成しません。
標準出力:
元の結晶の格子パラメータ、格子ベクトル、メートルテンソル、体積、原子サイトの分数座標を出力します。
適用される変換行列を表示します。
変換後の結晶の格子パラメータ、格子ベクトル、メートルテンソル、体積、原子サイトの分数座標を出力します。
描画のために収集された全ての元のサイトと変換後のサイトの詳細を出力します。
グラフィカルウィンドウ:
Matplotlibによって生成された3Dプロットウィンドウが表示されます。このウィンドウには、元の単位格子（黒い線）と、変換された単位格子（青い線と灰色の原子）が同時に描画されます。このウィンドウはユーザーが手動で閉じない限り表示され続けます。
コマンドラインでの使用例 (Usage)
プログラムの基本的な実行コマンドと引数は以下の通りです。
python crystal_convert_cell.py crystal_name conversion_mode kRatom
引数の意味は次の通りです。
crystal_name: 結晶名。'FCC','BCC','Hex','Rhomb'のいずれか。conversion_mode: 変換モード。'FCCPrim','BCCPrim','HexOrtho','HexRhomb','RhombHex','ACenterPrim','BCenterPrim','CCenterPrim'のいずれか。kRatom: 描画時の原子半径の係数（浮動小数点数）。
例:
python crystal_convert_cell.py FCC FCCPrim 1000
コマンドラインでの具体的な使用例
例1: FCC格子をプリミティブセルに変換
面心立方格子 (FCC) をそのプリミティブ（基本）セルに変換し、デフォルトの原子半径係数 (1000.0) で描画します。
python crystal_convert_cell.py FCC FCCPrim
実行結果:
標準出力には、元のFCC格子の情報（格子定数 \(a=5.62\) Å, \(\alpha=\beta=\gamma=90^\circ\)）と、変換後のプリミティブセルの情報が表示されます。プリミティブセルの格子定数は \(a=b=c \approx 3.974\) Å, \(\alpha=\beta=\gamma=60^\circ\) となるはずです。
出力の抜粋:
Crystal name: FCC
Lattice parameters: 5.6200 5.6200 5.6200 90.000 90.000 90.000
Lattice vectors:
ax: ( 5.6200, 0.0000, 0.0000) A
ay: ( 0.0000, 5.6200, 0.0000) A
az: ( 0.0000, 0.0000, 5.6200) A
...
Conversion mode: FCCPrim
(tij) = | 0.500, 0.500, 0.000|
| 0.000, 0.500, 0.500|
| 0.500, 0.000, 0.500|
Converted lattice vectors:
acx: ( 2.8100, 2.8100, 0.0000) A
acy: ( 0.0000, 2.8100, 2.8100) A
acz: ( 2.8100, 0.0000, 2.8100) A
Lattice parameters: 3.9744 3.9744 3.9744 60.000 60.000 60.000
...
グラフィカルウィンドウには、直方体のFCC単位格子（黒）と、その内部に位置する菱面体晶のプリミティブセル（青）が描画されます。
例2: BCC格子をプリミティブセルに変換し、原子サイズを小さく表示
体心立方格子 (BCC) をそのプリミティブセルに変換し、原子の描画サイズを小さく（kRatom=500.0）設定します。
python crystal_convert_cell.py BCC BCCPrim 500.0
実行結果:
標準出力には、元のBCC格子の情報と、変換後のプリミティブセルの情報が表示されます。BCCは \(a=b=c=5.62\) Å, \(\alpha=\beta=\gamma=90^\circ\) から、変換後のプリミティブセルは \(a=b=c \approx 4.867\) Å, \(\alpha=\beta=\gamma \approx 109.47^\circ\) となるはずです。
出力の抜粋:
Crystal name: BCC
Lattice parameters: 5.6200 5.6200 5.6200 90.000 90.000 90.000
...
Conversion mode: BCCPrim
(tij) = | -0.500, 0.500, 0.500|
| 0.500, -0.500, 0.500|
| 0.500, 0.500, -0.500|
Converted lattice vectors:
acx: (-2.8100, 2.8100, 2.8100) A
acy: ( 2.8100, -2.8100, 2.8100) A
acz: ( 2.8100, 2.8100, -2.8100) A
Lattice parameters: 4.8670 4.8670 4.8670 109.471 109.471 109.471
...
グラフィカルウィンドウには、直方体のBCC単位格子（黒）と、その内部に位置するプリミティブセル（体心立方格子の場合は逆平行六面体、青）が描画されます。原子は例1よりも小さく表示されます。
例3: 六方晶を斜方晶に変換
六方晶 (Hex) の単位格子を斜方晶 (Ortho) の単位格子に変換します。
python crystal_convert_cell.py Hex HexOrtho
実行結果:
標準出力には、元の六方晶格子の情報（\(a=b=5.62\) Å, \(c=4.00\) Å, \(\alpha=\beta=90^\circ, \gamma=120^\circ\)）と、変換後の斜方晶格子の情報が表示されます。
出力の抜粋:
Crystal name: Hex
Lattice parameters: 5.6200 5.6200 4.0000 90.000 90.000 120.000
...
Conversion mode: HexOrtho
(tij) = | 1.000, 0.000, 0.000|
| 1.000, 2.000, 0.000|
| 0.000, 0.000, 1.000|
Converted lattice vectors:
acx: ( 5.6200, 0.0000, 0.0000) A
acy: ( -5.6200, 9.7329, 0.0000) A
acz: ( 0.0000, 0.0000, 4.0000) A
Lattice parameters: 5.6200 11.2400 4.0000 90.000 90.000 90.000
...
グラフィカルウィンドウには、六方晶の単位格子（黒）と、それから構築される斜方晶の単位格子（青）が描画されます。