crystal_MP_simple.py 技術ドキュメント

プログラムの動作

crystal_MP_simple.py は、指定された結晶格子におけるマデルングポテンシャルの距離依存性を計算し、可視化するPythonプログラムです。主にイオン結晶において、ある中心イオンからの距離に応じて他のすべてのイオンが寄与する静電ポテンシャルを計算することを目的としています。

主な機能:

  1. 格子情報の計算と表示: 定義された格子定数とサイト情報から、実空間および逆空間の格子ベクトル、メトリックテンソル、単位胞体積などを計算し、標準出力に表示します。

  2. マデルングポテンシャル計算: 結晶内の基準となるイオン（デフォルトでは sites リストの最初のイオン）を中心とし、そこからの距離に応じた静電ポテンシャル（マデルングポテンシャル）を、指定された距離範囲とステップ数で計算します。この計算は、周期的に配置された多数の単位胞内のイオンからの寄与を直接合計することで行われます。

  3. 結果の可視化: 計算されたマデルングポテンシャルの距離依存性をグラフとしてプロットし、matplotlibウィンドウに表示します。

解決する課題:

イオン結晶におけるマデルングポテンシャルは、長距離にわたるクーロン相互作用の合計であり、その収束にはEwald和のような特殊な手法が用いられることが一般的です。本プログラムは、比較的シンプルな直接和の方法でマデルングポテンシャルを計算し、その距離による変化を視覚的に理解する手段を提供します。これにより、マデルングポテンシャルの基本的な概念とその計算アプローチを学ぶことができます。

原理

結晶の幾何学

プログラムは、格子定数 \(a, b, c, \alpha, \beta, \gamma\) を入力として、以下の結晶学的な情報を計算します。

  1. 実格子ベクトル: 直交座標系における格子ベクトル \(\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3\) を計算します。 例えば、 $\( \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{V_{cell}}{abc \sin\gamma}) \)\( ただし、\)V_{cell}$ は単位胞の体積です。

  2. メトリックテンソル: 格子ベクトル間の内積から構成されるメトリックテンソル \(g_{ij}\) を計算します。これは、結晶座標系における2点間の距離を計算する際に使用されます。 $\( g_{ij} = \mathbf{a}_i \cdot \mathbf{a}_j \)\( 2つの分率座標 \)\mathbf{u} = (u_1, u_2, u_3)\( と \)\mathbf{v} = (v_1, v_2, v_3)\( を持つ点の間の距離 \)d\( は、以下のように計算されます。 \)\( d = \sqrt{\sum_{i=1}^3 \sum_{j=1}^3 (u_i - v_i) (u_j - v_j) g_{ij}} \)$

  3. 逆格子ベクトル: 実格子ベクトルから、逆格子ベクトル \(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3\) を計算します。 $\( \mathbf{b}_1 = 2\pi \frac{\mathbf{a}_2 \times \mathbf{a}_3}{\mathbf{a}_1 \cdot (\mathbf{a}_2 \times \mathbf{a}_3)} \)\( 同様に \)\mathbf{b}_2, \mathbf{b}_3$ も計算されます。

マデルングポテンシャルの計算

プログラムは、ある中心イオン（sites[0]で定義される）の分率座標 \(\mathbf{r}_0\) を基準とし、結晶内の他のすべてのイオンからの静電ポテンシャルを合計することでマデルングポテンシャルを計算します。

各イオン \(j\) の電荷を \(q_j\) とし、その位置を \(\mathbf{r}_j\) とすると、中心イオン \(\mathbf{r}_0\) におけるポテンシャル \(V\) は以下のクーロンの法則に基づく合計として与えられます。

\[ V(\mathbf{r}_0) = \sum_{j \neq \text{center}} \frac{K_e q_j}{|\mathbf{r}_j - \mathbf{r}_0|} \]

ここで、\(K_e = \frac{e}{4 \pi \epsilon_0}\) はクーロン定数であり、\(e\) は電気素量、\(\epsilon_0\) は真空の誘電率です。 プログラム内では $K_e = e \cdot e / (4.0 \cdot \text{pi} \cdot e0)$ と定義されています。

計算アルゴリズム:

  1. 距離範囲とステップの設定: コマンドライン引数またはデフォルト値に基づいて、最小距離 rmin から最大距離 rmax までを nr 個のステップで分割し、距離のリスト rlist を作成します。

  2. 単位胞の範囲決定: rmax と格子定数に基づいて、計算に含める単位胞の最大範囲 (nxmax, nymax, nzmax) を決定します。これにより、中心イオンの周りの多くの隣接する単位胞を考慮します。

  3. 差分ポテンシャルの計算:

    • MPdiff という配列を初期化します。これは、各距離ステップ ir において、その距離範囲にあるすべてのイオンからのポテンシャルの寄与を合計するためのものです。

    • -nzmax から nzmax までの \(z\) 方向、-nymax から nymax までの \(y\) 方向、-nxmax から nxmax までの \(x\) 方向のすべての単位胞についてループします。

    • 各単位胞内のすべてのサイト (sites リストの全イオン) についてループし、そのイオン \(j\) と中心イオン \(0\) との間の距離 \(r = |\mathbf{r}_j - \mathbf{r}_0|\) を計算します。距離計算には、結晶のメトリックテンソル gij を考慮した distance 関数が使用されます。

    • 計算された距離 \(r\) が rmin と rmax の間にある場合、そのイオンからのポテンシャル寄与 \(K_e \frac{q_j}{r}\) を計算し、対応する MPdiff[ir] に加算します。プログラムでは、結果を電子ボルト (eV) 単位にするため、クーロン定数 \(K_e\) に電気素量 \(e\) を乗算し、\(r\) をメートル単位に変換（オングストロームから \(1.0e-10\) を乗算）した後、\(e\) で除算しています。

  4. 積算マデルングポテンシャルの計算: MP という配列を初期化し、MP[0] = MPdiff[0] とします。その後、ループで MP[i] = MP[i-1] + MPdiff[i] と計算することで、rmin から \(r_i\) までの距離範囲にあるすべてのイオンからのポテンシャルを累積的に合計し、マデルングポテンシャルを得ます。

  5. 結果の表示とプロット: 計算された距離 rlist と積算ポテンシャル MP を標準出力に表示し、matplotlib を使用してグラフを生成します。

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

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

  • numpy: 数値計算（配列操作、数学関数など）

  • matplotlib: データのプロットと可視化

  • tkcrystalbase: 結晶学的な基本計算を行うカスタムライブラリ。このファイルは、crystal_MP_simple.py と同じディレクトリ、またはPythonのパスが通っている場所に配置する必要があります。このライブラリのソースコードは本ドキュメントには含まれていませんが、結晶学の一般的な計算（格子ベクトル、メトリックテンソル、距離計算など）を提供します。

これらのライブラリは、Pythonのパッケージ管理システム pip を使用してインストールできます。

pip install numpy matplotlib

tkcrystalbase はカスタムライブラリのため、別途入手または作成し、プログラムが参照できる場所に配置してください。

必要な入力ファイル

このプログラムは、外部の入力ファイルを必要としません。 格子定数 (lattice_parameters) やサイト情報 (sites) は、プログラムのソースコード内に直接ハードコードされています。 もしこれらのデータを変更したい場合は、ソースコードを直接編集する必要があります。

生成される出力ファイル

このプログラムは、いかなるファイルも生成しません。

計算結果は以下の形式で標準出力に表示されます。

  • 格子定数、格子ベクトル、メトリックテンソル、単位胞体積

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

  • 計算に用いる単位胞の最大範囲

  • 距離 (r (A)) と対応するマデルングポテンシャル (Madelung potential (eV)) の表

さらに、計算されたマデルングポテンシャルの距離依存性を示すグラフが matplotlib を用いてGUIウィンドウに表示されます。

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

crystal_MP_simple.py は、以下の形式でコマンドラインから実行できます。

python crystal_MP_simple.py [rmax] [nr]

  • rmax: マデルングポテンシャルを計算する最大距離をオングストローム (A) 単位で指定します。これは必須の引数ではありません。指定しない場合、プログラム内のデフォルト値 (100.0) が使用されます。

  • nr: rmin (プログラム内のデフォルト値は 0.1 A) から rmax までの距離を分割するステップ数を指定します。これは必須の引数ではありません。指定しない場合、プログラム内のデフォルト値 (101) が使用されます。

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

例1: デフォルト値での実行

python crystal_MP_simple.py

実行結果の説明: このコマンドを実行すると、rmax=100.0 A、nr=101 (デフォルト値) でマデルングポテンシャルが計算されます。標準出力には、格子情報、マデルングポテンシャルの数値リストが表示され、距離に対するポテンシャルのグラフが新しいウィンドウに表示されます。グラフは、距離が増加するにつれてマデルングポテンシャルがある値に収束していく様子を示します。

例2: 最大距離とステップ数を指定して実行

python crystal_MP_simple.py 20.0 51

実行結果の説明: このコマンドは、マデルングポテンシャルの計算範囲を最大距離 20.0 A までに制限し、その範囲を 51 個のステップで分割して計算を実行します。 標準出力には、計算された格子情報と、0.1 A から 20.0 A までのマデルングポテンシャルの値が51行にわたって表示されます。 matplotlib ウィンドウには、0.1 A から 20.0 A までのマデルングポテンシャルのグラフが表示され、より短い距離範囲でのポテンシャルの収束挙動を詳細に観察することができます。