技術ドキュメント: crystal_MP_Evjen.py

プログラムの動作

crystal_MP_Evjen.py は、イオン結晶におけるマデルングポテンシャルをEvjen (エフイェン) の方法を用いて計算するPythonスクリプトです。このプログラムは、指定された格子パラメーターと単位セル内のサイト情報（原子種、電荷、位置など）に基づいて、結晶の構造情報を分析し、特定のイオン周囲のマデルングポテンシャルとマデルング定数を計算します。

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

  • 結晶構造の解析: 入力された格子パラメーターから、格子ベクトル、メトリックテンソル、単位セル体積、およびそれらの逆格子情報を計算し出力します。

  • サイト情報の処理: 単位セル内の全サイトをリストアップし、それぞれのマルチプリシティ（多重度）を考慮した重み付けを行い、総電荷のバランスを確認します。

  • マデルングポテンシャルの計算: 単位セル内の特定のイオン（デフォルトではリストの最初のイオン）を中心に据え、Evjenの方法に基づき周囲のイオンからの静電ポテンシャルを合計します。計算範囲はコマンドライン引数で指定可能です。

  • 結果の出力: 計算されたマデルングポテンシャル（ジュールおよび電子ボルト単位）とマデルング定数を標準出力に表示します。

このプログラムは、主にイオン結晶のエネルギー特性を理解するための基礎的な静電相互作用の計算ツールとして機能します。

原理

このプログラムは、イオン結晶中の特定のイオンが受ける静電ポテンシャルの総和であるマデルングポテンシャルを計算します。計算には、条件収束の問題を回避するためのEvjen法が用いられます。

  1. 静電ポテンシャル: クーロンの法則に基づき、電荷 \(q_j\) を持つイオンが距離 \(r_{ij}\) 離れた位置にあるイオン \(i\) に及ぼすポテンシャル \(V_{ij}\) は次のように表されます。 $\(V_{ij} = \frac{k_e q_j}{r_{ij}}\)\( ここで \)k_e\( はクーロン定数であり、\)k_e = \frac{1}{4 \pi \epsilon_0}\( です。プログラム内では電荷 \)e\( を用いて \)k_e = \frac{e^2}{4 \pi \epsilon_0}$ と表現されています。

  2. マデルングポテンシャル: イオン \(i\) が結晶全体から受けるマデルングポテンシャル \(V_i\) は、他のすべてのイオンからのポテンシャルの総和です。 $\(V_i = \sum_{j \neq i} \frac{k_e q_j}{r_{ij}}\)$ この無限級数の直接的な合計は条件収束であるため、適切な方法で計算する必要があります。

  3. Evjen の方法: Evjenの方法は、電荷中性の単位ブロックを定義し、それを広げていくことでポテンシャルを計算する手法の一つです。プログラムでは、単位セル内に存在する各サイトの電荷 \(q_j\) に、そのサイトの単位セルにおけるマルチプリシティ \(m_j\) の逆数 \(1/m_j\) を重みとして乗じることで、実効的に電荷中性の単位セルを構成しています。これにより、各単位セルが全体として電荷中性であるかのように扱われ、級数の収束が改善されます。計算範囲は nmax で指定される格子定数の倍数で区切られた直方体領域です。

  4. マデルング定数: マデルング定数 \(M\) は、マデルングポテンシャルを無次元化したものです。プログラムでは以下の式を用いて計算されます。 $\(M = \frac{1}{2} \cdot \frac{MP}{|q_0|} \cdot (a \cdot 10^{-10})\)\( ここで \)MP\( はクーロン定数 \)k_e\( を含まないポテンシャルの和 \)\sum_{j \neq i} \frac{q_j}{r_{ij}}\(、 \)q_0\( は基準となるイオンの電荷、\)a\( は結晶の代表的な格子定数（プログラムでは `lattice_parameters[0]` を使用）です。係数 \)1/2$ は、マデルング定数の定義やEvjen法における特定の補正に由来する可能性があります。

  5. 結晶学の基礎: プログラムでは、tkcrystalbase モジュール内の関数を使用して、格子パラメーターから直交座標系での格子ベクトル、メトリックテンソル、単位セル体積、およびそれらの逆格子表現を計算しています。

    • 格子ベクトル: 単位セルを定義する3つのベクトル \(a_x, a_y, a_z\)。

    • メトリックテンソル: \(g_{ij} = \vec{a_i} \cdot \vec{a_j}\) で定義され、結晶座標系における距離計算に使用されます。

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

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

  • NumPy: 数値計算（配列操作、線形代数など）に広く使われるライブラリ。

  • Matplotlib: グラフ描画ライブラリ（本プログラムでは直接的な描画は行われていませんが、インポートされています）。

これらのライブラリは、pip を用いてインストールできます。

pip install numpy matplotlib

また、このプログラムは tkcrystalbase というモジュールに依存しています。このモジュールは標準のPythonパッケージインデックス (PyPI) には登録されていません。したがって、このモジュールは crystal_MP_Evjen.py と同じディレクトリに配置されているか、Pythonのモジュール探索パスが通っている必要があります。別途、このモジュールのソースコードを入手し、適切に配置してください。

必要な入力ファイル

crystal_MP_Evjen.py は、外部からの入力ファイルを必要としません。格子パラメーター、単位セル内のサイト情報、計算範囲のデフォルト値は、すべてプログラムのソースコード内にハードコードされています。

  • lattice_parameters: 格子定数 (a, b, c) と格子角 (alpha, beta, gamma) をリストとして定義。

  • sites: 単位セル内の各原子サイトの詳細（原子名、ラベル、原子番号、質量、電荷、半径、色、分率座標）をリストのリストとして定義。

  • nmax: マデルングポテンシャル計算の範囲を決定する整数。デフォルト値は 1 ですが、コマンドライン引数で上書き可能です。

  • rmin: 同一サイトと見なすための最小距離。

将来的にこれらのデータを外部ファイル（例: CIF, POSCAR, YAML, JSONなど）から読み込むように拡張することも可能です。

生成される出力ファイル

crystal_MP_Evjen.py は、計算結果をファイルとして保存しません。すべての出力は標準出力 (stdout) に直接表示されます。

出力される情報には、以下の内容が含まれます。

  • 入力された格子パラメーター

  • 計算された格子ベクトル、メトリックテンソル、単位セル体積

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

  • 単位セル内の拡張されたサイトリストとそれらのマルチプリシティ、および総電荷

  • 中心イオン周囲のマデルングポテンシャル（ジュールおよび電子ボルト単位）

  • 計算されたマデルング定数

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

crystal_MP_Evjen.py を実行する際の基本的なコマンドライン構文は以下の通りです。

python crystal_MP_Evjen.py nmax

  • nmax: Evjen法の計算で考慮する単位セルの範囲を指定する整数値です。例えば 1 と指定すると、中心セルとその隣接する1層のセル（\([-1, 1]\) の範囲のセル）が考慮されます。値が大きいほど計算範囲が広がり、より正確な結果が得られますが、計算時間も増加します。

例:

python crystal_MP_Evjen.py 1

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

ここでは、nmax に 1 を指定して crystal_MP_Evjen.py を実行する具体的な例と、その実行結果の一部を示します。

実行コマンド:

python crystal_MP_Evjen.py 1

実行結果の例:

（実際の出力はシステムや計算精度によって若干異なる場合があります。）

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.117436067711466, 1.117436067711466, 1.117436067711466, 90.0, 90.0, 90.0]
Reciprocal lattice vectors:
  Rax: (   1.117,        0,        0) A^-1
  Ray: (       0,    1.117,        0) A^-1
  Raz: (       0,        0,    1.117) A^-1
Reciprocal lattice metric tensor:
  Rgij: (   1.249,        0,        0) A^-1
        (       0,    1.249,        0) A^-1
        (       0,        0,    1.249) A^-1
Reciprocal unit cell volume:   0.005639 A^-3

Site information (all sites in the unit cell with the range: [[-0.1, 1.1], [-0.1, 1.1], [-0.1, 1.1]])
  Na1 (   0.0,    0.0,    0.0) q=  1.0 mult= 1 weight=   1.0000
  Na2 (   0.0,    0.5,    0.5) q=  1.0 mult= 1 weight=   1.0000
  Na3 (   0.5,    0.0,    0.5) q=  1.0 mult= 1 weight=   1.0000
  Na4 (   0.5,    0.5,    0.0) q=  1.0 mult= 1 weight=   1.0000
  Cl1 (   0.5,    0.0,    0.0) q= -1.0 mult= 1 weight=   1.0000
  Cl2 (   0.5,    0.5,    0.5) q= -1.0 mult= 1 weight=   1.0000
  Cl3 (   0.0,    0.0,    0.5) q= -1.0 mult= 1 weight=   1.0000
  Cl4 (   0.0,    0.5,    0.0) q= -1.0 mult= 1 weight=   1.0000
qtot= 0.0

Calculate Madelung potential around the zero-th ion by Evjen method
  nmax: 1
  Origin: Na1 (0.0, 0.0, 0.0)

  Madelung potential: -4.03714e-19 J
  Madelung potential: -2.52042 eV
  Madelung constant:  -1.74756314

この例では、NaCl構造の格子パラメーターが与えられ、nmax=1 で計算された結果が表示されています。中心の Na1 イオンにおけるマデルングポテンシャルは約 \(-2.52\) eV、マデルング定数は約 \(-1.74756\) と算出されています。このマデルング定数の値は、NaCl構造の理論値（約 \(1.74756\)) と一致しており、Evjen法の有効性を示しています。