crystal_draw_cell.py 技術ドキュメント

プログラムの動作

crystal_draw_cell.py は、指定された格子定数と原子サイト情報に基づき、結晶の実空間単位格子と逆空間単位格子を3Dで可視化するPythonスクリプトです。このプログラムは、物理的な格子構造を視覚的に理解することを目的としており、原子サイトの配置とそれに対応する逆格子空間の構造を同時に表示します。

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

  1. 格子情報の計算: 入力された格子定数から、実空間の格子ベクトル、単位格子体積、メトリックテンソルを計算します。

  2. 逆格子情報の計算: 計算された実空間格子ベクトルから、逆空間の格子ベクトル、逆格子定数、逆格子体積、逆格子メトリックテンソルを計算します。

  3. 原子サイトの生成と重複判定: 定義された単位格子内の原子サイト情報に基づき、指定された範囲内で周期的に原子サイトを生成します。この際、rmin で指定された距離閾値内にあるサイトは同一とみなし、重複を排除します。

  4. 3D可視化: Matplotlibライブラリを使用して、計算された実空間単位格子（原子サイトを含む）と、適切なスケールで表示された逆空間単位格子を3Dグラフィックで描画します。

このプログラムは、結晶学の教育や研究において、実空間と逆空間の関係を直感的に把握するためのツールとして役立ちます。プログラムの冒頭にある3重引用符で囲まれたコメント（docstring）には、「単位格子と逆格子単位格子を描画する」という目的と「tkcrystalbase.py が必要であること」が簡潔に記述されています。

原理

crystal_draw_cell.py は、結晶学の基本的な数理に基づいています。主要な計算原理は以下の通りです。

1. 格子ベクトル、体積、メトリックテンソル

結晶の格子定数 \(a, b, c\) (長さ) および \(\alpha, \beta, \gamma\) (角度) から、実空間の基本格子ベクトル \(\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3\) がデカルト座標系で計算されます。tkcrystalbase モジュール内の cal_lattice_vectors 関数がこの計算を行います。一般的に、これらは以下のように表現できます。

\[\begin{split}\begin{aligned} \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 \frac{\sqrt{1 - \cos^2\alpha - \cos^2\beta - \cos^2\gamma + 2\cos\alpha\cos\beta\cos\gamma}}{\sin \gamma}) \end{aligned}\end{split}\]

単位格子の体積 \(V\) は、これら3つの格子ベクトルの混合積として計算されます。

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

または、格子定数から直接計算する式もあります。

\[V = abc \sqrt{1 - \cos^2\alpha - \cos^2\beta - \cos^2\gamma + 2\cos\alpha\cos\beta\cos\gamma}\]

メトリックテンソル \(g_{ij}\) は、格子ベクトルの内積によって定義され、結晶内の距離計算に不可欠です。

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

行列形式では以下のようになります。

\[\begin{split} G = \begin{pmatrix} a^2 & ab \cos\gamma & ac \cos\beta \\ ab \cos\gamma & b^2 & bc \cos\alpha \\ ac \cos\beta & bc \cos\alpha & c^2 \end{pmatrix} \end{split}\]

tkcrystalbase の cal_metrics 関数は、このメトリックテンソルなどを計算します。

2. 逆格子ベクトルと逆格子体積

実空間の格子ベクトル \(\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3\) から、逆空間の基本格子ベクトル \(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3\) が計算されます。結晶学の分野では、通常 \(2\pi\) の因子を省略して定義されることが多く、このプログラムもその定義を採用していると推測されます。

\[\begin{split}\begin{aligned} \mathbf{b}_1 &= \frac{\mathbf{a}_2 \times \mathbf{a}_3}{V} \\ \mathbf{b}_2 &= \frac{\mathbf{a}_3 \times \mathbf{a}_1}{V} \\ \mathbf{b}_3 &= \frac{\mathbf{a}_1 \times \mathbf{a}_2}{V} \end{aligned}\end{split}\]

ここで \(V\) は実空間の単位格子体積です。逆格子単位格子の体積 \(V^*\) は、実空間単位格子体積の逆数に等しくなります。

\[V^* = |\mathbf{b}_1 \cdot (\mathbf{b}_2 \times \mathbf{b}_3)| = \frac{1}{V}\]

tkcrystalbase の cal_reciprocal_lattice_vectors および cal_reciprocal_lattice_parameters 関数がこれらの計算を行います。

3. 原子サイトの正規化と重複判定

原子サイトの座標は、単位格子内での相対座標 \((u, v, w)\) で与えられます。reduce01 関数は、任意の浮動小数点数を \(0 \le x < 1\) の範囲に正規化します。これはPythonの % 演算子を使用して実現されます（例: x % 1.0）。

add_site 関数では、nrange で指定された範囲で生成された新しいサイトが、既存のサイトと重複していないかどうかの判定が行われます。2つのサイト \(\mathbf{r}_1, \mathbf{r}_2\) 間の距離は、実空間のメトリックテンソル \(G\) を用いた以下の式で計算されます。

\[d(\mathbf{r}_1, \mathbf{r}_2) = \sqrt{(\mathbf{r}_1 - \mathbf{r}_2)^T G (\mathbf{r}_1 - \mathbf{r}_2)}\]

もしこの距離が rmin (距離の閾値) 以下であれば、それらのサイトは同一であると判断され、リストには追加されません。

4. 3D描画

Matplotlibの mplot3d ツールキットを用いて、実空間の単位格子（および原子）と逆空間の単位格子が描画されます。逆空間単位格子は kRUC 係数によって実空間単位格子に対する相対的な大きさが調整されます。draw_unitcell 関数がこれらの描画処理を担当します。

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

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

  1. NumPy: 数値計算（ベクトル、行列演算、三角関数など）のために必要です。

    • インストールコマンド: pip install numpy

  2. Matplotlib: 3Dグラフの描画のために必要です。プログラム内で直接 import matplotlib.pyplot as plt は行われていませんが、plt.figure などの関数が使用されているため、tkcrystalbase.py モジュール内でインポートされていると推測されます。

    • インストールコマンド: pip install matplotlib

  3. tkcrystalbase.py: このプログラムの主要な計算ロジック（格子定数から格子ベクトルへの変換、メトリックテンソルの計算、逆格子ベクトルの計算、サイトの重複判定、描画ヘルパー関数など）を提供するカスタムモジュールです。

    • インストールは不要ですが、crystal_draw_cell.py と同じディレクトリにこのファイルが存在する必要があります。

これらのライブラリをインストールするには、コマンドプロンプトまたはターミナルで以下のコマンドを実行してください。

pip install numpy matplotlib

tkcrystalbase.py は、このプログラムが動作するために必須です。このファイルが見つからない場合、プログラムは ModuleNotFoundError を発生させます。

必要な入力ファイル

crystal_draw_cell.py は、外部の入力ファイルを必要としません。すべての設定データはプログラムのソースコード内に直接ハードコードされています。

具体的には、以下の変数を変更することで、描画する結晶構造をカスタマイズできます。

  • lattice_parameters (リスト):

    • 格子定数 \(a, b, c\) (単位: オングストローム, A) と、角度 \(\alpha, \beta, \gamma\) (単位: 度数法) を指定します。

    • 例: [a, b, c, alpha, beta, gamma]

  • sites (リストのリスト):

    • 単位格子内の原子サイト情報を定義します。各サイトは以下の要素を持つリストです。

      • atom name (文字列): 原子名（例: 'Na', 'Cl'）

      • site label (文字列): サイトの識別子（例: 'Na1', 'Cl1'）

      • atomic number (整数): 原子番号

      • atomic mass (浮動小数点数): 原子質量

      • charge (浮動小数点数): 形式電荷

      • radius (浮動小数点数): 原子半径（描画時のサイズに影響）

      • color (文字列): 描画時の原子の色（例: 'red', 'blue'）

      • position (NumPy配列): 単位格子内での相対座標 [u, v, w]

これらの変数は、main 関数の実行前にプログラムコードの冒頭部分で定義されています。

その他の設定変数（kr, rmin, kRUC, nrange, figsize）も同様にコード内で直接変更可能です。

生成される出力ファイル

crystal_draw_cell.py は、実行時に直接的な出力ファイルを生成しません。

プログラムの出力は以下の2種類です。

  1. 標準出力へのテキスト情報:

    • 計算された格子パラメータ、格子ベクトル、メトリックテンソル、単位格子体積。

    • 逆格子パラメータ、逆格子ベクトル、逆格子メトリックテンソル、逆格子体積。

    • 描画される全原子サイトのリスト（原子名、サイトラベル、座標、電荷）。

    • これらの情報は、プログラムが正常に動作しているか、計算結果が期待通りかを確認するために利用されます。

  2. Matplotlibによる3Dプロットウィンドウ:

    • 実空間の単位格子（黒い線）と、その中に配置された原子サイト（色付きの球）が描画されます。

    • 逆空間の単位格子（赤い線）も描画され、実空間単位格子との相対的なスケールが kRUC で調整されます。

    • このウィンドウはインタラクティブであり、ユーザーは視点を変更したり、画像を保存したりすることができます（Matplotlibの標準機能）。

プログラムの実行が完了すると、プロットウィンドウが閉じられ、標準出力に何も問題がなければプログラムは終了します。

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

crystal_draw_cell.py は、コマンドライン引数を必要としません。プログラムの実行は非常にシンプルです。

プログラムが保存されているディレクトリに移動し、以下のコマンドを実行してください。

python crystal_draw_cell.py

このコマンドを実行すると、プログラムが起動し、標準出力に計算情報が表示された後、Matplotlibによる3Dプロットウィンドウが表示されます。

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

crystal_draw_cell.py を実行する具体的な例を以下に示します。

まず、crystal_draw_cell.py とその依存モジュール tkcrystalbase.py が同じディレクトリにあることを確認します。

# プログラムを実行する
python crystal_draw_cell.py

実行結果

上記のコマンドを実行すると、ターミナルまたはコマンドプロンプトに以下のような情報が出力されます（表示される数値はコード内の lattice_parameters や sites の設定に依存します）。

Lattice parameters: [5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
Lattice vectors:
  ax: (    5.62,     0.00,     0.00) A
  ay: (   2.810,    4.866,     0.00) A
  az: (   2.810,    1.622,    4.621) A
Metric tensor:
  gij: (    31.58,     15.80,     15.80) A
       (    15.80,     31.58,     15.80) A
       (    15.80,     15.80,     31.58) A
Volume:      102.7 A^3

Unit cell volume:      102.7 A^3
Reciprocal lattice parameters: [ 0.2312, 0.2312, 0.2312, 60.0, 60.0, 60.0]
Reciprocal lattice vectors:
  Rax: (  0.2312,   0.0000,   0.0000) A^-1
  Ray: (  0.1156,   0.2002,   0.0000) A^-1
  Raz: (  0.1156,   0.0667,   0.1905) A^-1
Reciprocal lattice metric tensor:
  Rgij: (  0.0535,   0.0268,   0.0268) A^-1
        (  0.0268,   0.0535,   0.0268) A^-1
        (  0.0268,   0.0268,   0.0535) A^-1
Reciprocal unit cell volume:  0.009406 A^-3

All sites to draw:
  Na  : Na1 : (   0.000,    0.000,    0.000) Z=   1.000
  Na  : Na2 : (   0.000,    0.500,    0.500) Z=   1.000
  Na  : Na3 : (   0.500,    0.000,    0.500) Z=   1.000
  Na  : Na4 : (   0.500,    0.500,    0.000) Z=   1.000
  Cl  : Cl1 : (   0.500,    0.000,    0.000) Z=  -1.000
  Cl  : Cl2 : (   0.500,    0.500,    0.500) Z=  -1.000
  Cl  : Cl3 : (   0.000,    0.000,    0.500) Z=  -1.000
  Cl  : Cl4 : (   0.000,    0.500,    0.000) Z=  -1.000
# ... (nrangeに応じてさらに多くのサイト情報が出力されます)

同時に、Matplotlibによって新しいウィンドウがポップアップ表示されます。このウィンドウには、以下のような3Dグラフィックが表示されます。

  • 実空間単位格子: 黒い線で描画された格子フレーム。この例では、NaCl構造に対応する原子サイトが配置されています。

  • 原子サイト: sites リストで定義されたNa (赤色の球) とCl (青色の球) が、それぞれの半径と色で描画されます。nrange の設定によっては、複数の単位格子にまたがる原子が描画されることがあります。

  • 逆空間単位格子: 赤い線で描画された格子フレーム。実空間単位格子とは異なる向きとスケールで描画され、実空間と逆空間の関係を視覚的に示します。

ユーザーはこのプロットウィンドウをマウスで操作し、視点を回転させたりズームイン/アウトしたりして、結晶構造を多角的に観察することができます。プロットウィンドウを閉じると、プログラムは終了します。