MP_Ewald.py 技術ドキュメント

プログラムの動作

MP_Ewald.py は、結晶構造データ（CIF形式）からマデルングポテンシャルをEwald法によって計算するPythonスクリプトです。このプログラムの主な目的は、結晶中の各原子サイトにおける静電ポテンシャル、および単位格子あたりの全マデルングエネルギーを精度良く算出することです。

主な機能:

  • CIFファイルの読み込みと解析: 入力されたCIFファイルを解析し、結晶の格子パラメータ、原子サイト、原子の電荷情報などを抽出します。

  • Ewald法の適用: 長距離クーロン相互作用の合計を効率的に計算するためにEwald法を使用します。これは、実空間和（短距離相互作用）、逆空間和（長距離相互作用）、および自己エネルギー項の3つの部分に分けて計算されます。

  • Ewaldパラメータの最適化: ユーザーが指定するalpha（Ewaldパラメータ）とprec（計算精度）に基づいて、実空間および逆空間のカットオフ半径を決定し、計算の収束を確保します。

  • マデルングポテンシャルの算出: 各原子サイトにおけるマデルングポテンシャル（単位電荷あたりの電位）を計算し、ジュール(J)および電子ボルト(eV)単位で出力します。

  • 全マデルングエネルギーの計算: 単位格子内の全マデルングエネルギー（eV単位）を計算します。

  • マデルング定数の計算: 代表的な電荷と参照距離（a軸長）を用いて、結晶のマデルング定数を算出します。

  • ログ出力: 計算の詳細、パラメータ、結果を標準出力と指定されたログファイルに出力します。

解決する課題:

結晶中の静電相互作用は長距離にわたるため、単純に全ての原子間のクーロン力を合計すると収束が非常に遅くなるか、発散することがあります。Ewald法は、この長距離クーロン相互作用の合計を数学的に収束する3つの項（実空間和、逆空間和、自己エネルギー項）に分割することで、効率的かつ高精度にマデルングポテンシャルを計算することを可能にします。これにより、イオン結晶の安定性や物性を理解するための基礎データを提供します。

原理

MP_Ewald.py は、イオン結晶におけるマデルングポテンシャルをEwald法を用いて計算します。マデルングポテンシャルは、結晶中のあるイオンが、その結晶中の他の全てのイオンから受ける静電ポテンシャルの合計です。Ewald法は、無限に広がる周期的な結晶におけるクーロン相互作用の合計を高速に収束させるための強力な手法であり、以下の3つの項に分解して計算します。

  1. 実空間和 (Real-space sum): 短距離的な相互作用を記述する項です。ガウス関数 \(erfc(\alpha r)\) （相補誤差関数）を用いて、急速に収束する部分を実空間で計算します。 各サイト \(i\) における実空間和 \(U_{real}(i)\) は、以下の形式で計算されます。 $\(U_{real}(i) = \sum_{j, \mathbf{n}}' \frac{q_j \cdot erfc(\alpha |\mathbf{r}_{ij} + \mathbf{n}|)}{|\mathbf{r}_{ij} + \mathbf{n}|}\)\( ここで、\)q_j\( はサイト \)j\( の電荷、\)\mathbf{r}_{ij}\( はサイト \)i\( からサイト \)j\( への相対位置ベクトル、\)\mathbf{n}\( は格子ベクトルを表します。\)\sum'\( は \)i=j\( かつ \)\mathbf{n}=\mathbf{0}\( の場合を除くことを意味します。\)\alpha$ はEwaldパラメータと呼ばれる収束パラメータで、実空間と逆空間の寄与のバランスを調整します。 コード中の UC1_list[isite] += qj * erfcar / (rij * 1.0e-10) に対応します。

  2. 逆空間和 (Reciprocal-space sum): 長距離的な相互作用を記述する項で、フーリエ空間（逆格子空間）で計算されます。実空間和と同様に、ガウス関数 \(\exp(-\frac{G^2}{4\alpha^2})\) を用いて急速に収束します。 各サイト \(i\) における逆空間和 \(U_{reciprocal}(i)\) は、以下の形式で計算されます。 $\(U_{reciprocal}(i) = \frac{1}{\pi V} \sum_{\mathbf{G} \neq \mathbf{0}} \frac{\exp\left(-\frac{\pi^2 G^2}{\alpha^2}\right)}{G^2} \left[ \sum_j q_j \cos(2\pi \mathbf{G} \cdot (\mathbf{r}_j - \mathbf{r}_i)) \right]\)\( ここで、\)V\( は単位格子の体積、\)\mathbf{G}\( は逆格子ベクトル、\)G = |\mathbf{G}|\( です。 コード中の `Krec` は \)1/(\pi V)\(、\)Kexp = \pi^2 / \alpha^2\(、\)expg = \exp(-Kexp \cdot G^2) / G^2\( に対応します。`fcal` は \)\sum_j q_j \cos(2\pi \mathbf{G} \cdot (\mathbf{r}_j - \mathbf{r}_i))$ の部分を計算しています。

  3. 自己エネルギー項 (Self-energy term): 各イオンが自分自身と相互作用していると見なされる無限大の寄与を相殺するための補正項です。 各サイト \(i\) における自己エネルギー項 \(U_{self}(i)\) は、以下の形式で計算されます。 $\(U_{self}(i) = - \frac{2 \alpha}{\sqrt{\pi}} q_i\)$ コード中の UC3_list[isite] = 2.0 * qi * (ew_alpha * 1.0e10) / sqrt(pi) に対応します。

これらの3つの項の合計が、各サイトのマデルングポテンシャル \(V_{M,i}\) となります。 $\(V_{M,i} = \frac{1}{4\pi\epsilon_0} (U_{real}(i) + U_{reciprocal}(i) - U_{self}(i))\)\( ただし、プログラム内部では \)q_i\( を原子の電荷数（無次元）とし、最終的な出力で電気素量 \)e\( とクーロン定数 \)K_e = e^2 / (4\pi\epsilon_0)$ を用いて適切な単位（JouleまたはeV）に変換しています。

クーロン定数に関連するコード内の定数:

  • Ke = e * e / 4.0 / pi / e0: \(e^2 / (4\pi\epsilon_0)\) に相当し、単位は \(N \cdot m^2\) です。

  • Ke / e: 物理定数 \(e / (4\pi\epsilon_0)\) に相当し、単位は \(J/C = V\) です。この値が UC_list の各項に乗じられ、ポテンシャル（ボルト）が計算されます。

計算精度とカットオフ: prec パラメータは計算の許容誤差を指定します。この精度に基づいて、実空間の最大距離 rdmax と逆空間の最大逆格子ベクトル G2max が決定され、不要な計算を削減しながら必要な精度を確保します。 距離計算には、格子定数 gij (実空間) および Rgij (逆空間) を用いて、周期境界条件を考慮した正しい距離が使用されます。

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

MP_Ewald.py は、以下の非標準ライブラリに依存しています。

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

  • matplotlib: グラフ描画ライブラリですが、このスクリプトでは直接的なプロット出力には使用されていません。

  • tklib (カスタムライブラリ): 結晶構造の読み込み、解析、ユーティリティ機能を提供します。具体的には tklib.tkapplication, tklib.tkutils, tklib.tkfile, tklib.tkcrystal (特に tklib.tkcrystal.tkcif, tklib.tkcrystal.tkcrystal, tklib.tkcrystal.tkatomtypeobject, tklib.tkcrystal.tkatomtype) がインポートされています。

  • tkcrystalbase (カスタムライブラリ): 結晶学的な基底計算（距離計算など）を提供するライブラリ。

インストール方法:

numpy および matplotlib は、Pythonのパッケージ管理システム pip を用いてインストールできます。

pip install numpy matplotlib

tklib および tkcrystalbase は、汎用的なPyPIリポジトリには登録されていない、特定のプロジェクトで用いられるカスタムライブラリであると推測されます。 したがって、これらのライブラリは別途入手・インストールが必要です。通常、これらはプログラムの配布元から提供されるか、特定のディレクトリに配置されていることを前提とします。

必要な入力ファイル

プログラムは、結晶構造情報が記述された CIF (Crystallographic Information File) 形式のファイルを必要とします。

  • ファイル形式: CIFファイル

  • ファイル名: コマンドライン引数でパスを指定します。

  • データ構造: 結晶の格子パラメータ（a, b, c, alpha, beta, gamma）、原子のサイト座標、原子の種類、および特に 各原子サイトの電荷 (_atom_site_charge) が含まれている必要があります。電荷情報がない場合、プログラムは電荷を既定値（陽イオンは+1、陰イオンは-1など）と仮定するか、エラーを発生させる可能性があります（このコードでは陽イオンを+0.7、陰イオンを-1.4と仮定している部分があるが、atom.Charge() が優先される）。

例 (CIFファイルの一部):

# ... (他のCIF情報) ...
_cell_length_a                   5.62
_cell_length_b                   5.62
_cell_length_c                   5.62
_cell_angle_alpha                60.0
_cell_angle_beta                 60.0
_cell_angle_gamma                60.0
# ...
loop_
 _atom_site_label
 _atom_site_type_symbol
 _atom_site_fract_x
 _atom_site_fract_y
 _atom_site_fract_z
 _atom_site_charge
 Na1 Na 0.0000 0.0000 0.0000 +1.0
 Cl1 Cl 0.5000 0.5000 0.5000 -1.0
# ...

生成される出力ファイル

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

  • ログファイル:

    • ファイル名: 入力CIFファイルと同じディレクトリに、元のファイル名に -out.txt を付け加えた名前で生成されます。 例: 入力ファイルが NaCl.cif の場合、NaCl-out.txt が生成されます。

    • 内容:

      • 入力ファイルパス、Ewaldパラメータ (alpha, prec) の情報。

      • 結晶の格子情報（格子パラメータ、体積など）。

      • Ewald法におけるカットオフパラメータ（RDmax, nrmax, G2max, hgmax）。

      • 実空間和、逆空間和、自己エネルギー項の計算に要した時間。

      • 各原子サイトにおけるマデルングポテンシャル（Joule および eV 単位）。

      • 単位格子あたりの全マデルングエネルギー（eV 単位）。

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

      • その他、プログラムの実行状況やデバッグ情報。

また、上記のすべての情報は標準出力 (コンソール) にもリアルタイムで出力されます。

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

MP_Ewald.py は、コマンドライン引数として入力CIFファイルのパス、Ewaldパラメータ alpha、および計算精度 prec を受け取ります。

python MP_Ewald.py <CIF_path> <alpha> <prec>

  • <CIF_path>: 入力するCIFファイルへのパス（必須）。

  • <alpha>: Ewaldパラメータ。実空間と逆空間の和の収束速度を調整する正の浮動小数点数（オプション、デフォルトは 0.3）。

  • <prec>: 計算精度。計算の打ち切り基準となる正の浮動小数点数（オプション、デフォルトは 1.0e-5）。

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

例えば、カソードとして用いられることが多いLiFePO4結晶のCIFファイル LiFePO4.cif があると仮定し、デフォルトのEwaldパラメータと精度でマデルングポテンシャルを計算する場合、以下のように実行します。

python MP_Ewald.py LiFePO4.cif

または、Ewaldパラメータ alpha を 0.4 に、計算精度 prec を 1.0e-6 に変更して実行する場合:

python MP_Ewald.py LiFePO4.cif 0.4 1.0e-6

実行結果の説明 (例):

上記のコマンドを実行すると、標準出力および LiFePO4-out.txt に以下のような情報が出力されます（具体的な数値はファイルの内容によります）。

Open logfile [./LiFePO4-out.txt]

=============== CIF file read test ============
input: LiFePO4.cif
log file: ./LiFePO4-out.txt
ew_alpha: 0.4
prec: 1e-06

Read [LiFePO4.cif]

Lattice parameters: [ 10.435, 6.096, 4.757, 90.0, 90.0, 90.0]

Ewald parameters
  alpha: 0.4
  precision = 1e-06 = 10^-6.0
  RDmax = 13.9 A, where erfc(alpha*RDmax) = 2.0e-07
  nrmax: 1 2 2
  cal_N(real): 230
  G2max: 13.67
      exp(-pi^2 * G2max^2 / alpha^2) =  2.8e-07
  hgmax: 2 3 4
  cal_N(reciprocal): 735

Time for real space sum     : 0.12345
Time for real reciprocal sum: 0.23456
Total time                  : 0.35801

Madelung potential (electrostatic potential):
Li  : Li1   : q=    1.0: MP = -1.000000e-18 J  (= -2.000000e-18 + 1.000000e-18 - 0.000000e+00)
                          MP = -6.241509e+00 eV (= -1.248302e+01 + 6.241509e+00 - 0.000000e+00)
Fe  : Fe1   : q=    2.0: MP = -2.000000e-18 J  (= -4.000000e-18 + 2.000000e-18 - 0.000000e+00)
                          MP = -1.248302e+01 eV (= -2.496604e+01 + 1.248302e+01 - 0.000000e+00)
P   : P1    : q=    5.0: MP = -3.000000e-18 J  (= -6.000000e-18 + 3.000000e-18 - 0.000000e+00)
                          MP = -1.872453e+01 eV (= -3.744907e+01 + 1.872453e+01 - 0.000000e+00)
O   : O1    : q=   -2.0: MP =  1.000000e-18 J  (=  2.000000e-18 - 1.000000e-18 + 0.000000e+00)
                          MP =  6.241509e+00 eV (=  1.248302e+01 - 6.241509e+00 + 0.000000e+00)
... (他のサイトについても同様に出力) ...

Total Madelung energy in unit cell: -1.234567e+02 eV

Madelung constant
NOTE: The a-axis length is taken as the representative atomic distance: a = 10.435
      The charge of the 0-th ion is taken as the representative ion charge: q = 1.0
NOTE: This value is in the unit cell chemical formula
       The following value must be devided by Z to get the Madeluing constant in the standard definition
Madelung constant in unit cell:   2.34567890