crystal_MP_Ewald.py 技術ドキュメント
プログラムの動作
crystal_MP_Ewald.py は、イオン結晶におけるMadelungポテンシャルおよびMadelung定数をEwald和の手法を用いて計算するPythonプログラムです。
このプログラムは、指定された格子定数と単位胞内の原子サイト情報（原子種、電荷、位置）に基づいて、結晶中の特定のサイトにおける静電ポテンシャルを決定します。
主な機能は以下の通りです。
結晶の格子ベクトル、計量テンソル、単位胞の体積、およびそれらの逆格子空間における対応する量を計算します。
Ewaldパラメータ (\(\alpha\)) と計算精度 (
prec) に基づいて、実空間および逆空間の和の計算範囲を自動的に決定します。Ewald和を構成する3つの主要な項（実空間項、逆空間項、自己項）を計算し、それらを合計してMadelungポテンシャルを導出します。
計算されたMadelungポテンシャルをジュール（J）および電子ボルト（eV）単位で出力します。
Madelung定数を計算し、出力します。
各計算ステップにかかる時間を測定し、出力します。
このプログラムは、イオン結晶のクーロン相互作用の合計が無限大に発散する問題を、Ewald和という強力なアルゴリズムを用いて収束させ、物理的に意味のある結晶内部のポテンシャルを算出することを目的としています。
原理
crystal_MP_Ewald.py は、無限に広がるイオン結晶のMadelungポテンシャルを計算するために、Ewald和のアルゴリズムを利用しています。Ewald和は、長距離にわたるクーロン相互作用の合計がゆっくりと収束するため、これを高速かつ正確に計算するために開発された手法です。
Ewald和では、結晶内の点電荷 \(q_j\) を、短距離で急激に減衰するガウス関数で広がった電荷分布 \(\rho_G(\mathbf{r})\) と、それを打ち消すような逆符号のガウス分布を持つ補償電荷分布 \(\rho_{comp}(\mathbf{r})\) に分割します。これにより、静電ポテンシャルは以下の3つの項の和として計算されます。
実空間項 (\(V_{real}\)) 広がったガウス電荷間の相互作用の合計です。これは短距離で収束が速く、実空間で直接計算されます。プログラムでは
UC1に対応します。\[V_{real}(\mathbf{r}_i) = \sum_{j, \mathbf{n}}' \frac{q_j \cdot \text{erfc}(\alpha |\mathbf{r}_{ij} + \mathbf{n}|)}{|\mathbf{r}_{ij} + \mathbf{n}|}\]ここで、\(q_j\) はサイト \(j\) の電荷、\(erfc\) は相補誤差関数、\(\alpha\) はEwaldパラメータ、\(\mathbf{r}_{ij}\) はサイト \(i\) から \(j\) への距離ベクトル、\(\mathbf{n}\) は格子ベクトルです。プライム記号 \('\) は、\(i=j\) かつ \(\mathbf{n}=0\) の項を除外することを示します。
逆空間項 (\(V_{reciprocal}\)) 点電荷と、その点電荷がガウス分布で広がったものとの差（補償電荷）の相互作用のフーリエ変換の合計です。これは長距離で収束が速く、逆格子空間で計算されます。プログラムでは
UC2に対応します。\[V_{reciprocal}(\mathbf{r}_i) = \frac{1}{\pi V} \sum_{\mathbf{G} \neq 0} \frac{e^{-\mathbf{G}^2 / (4\alpha^2)}}{\mathbf{G}^2} \left[ \sum_j q_j \cos(\mathbf{G} \cdot (\mathbf{r}_j - \mathbf{r}_i)) \right]\]ここで、\(V\) は単位胞の体積、\(\mathbf{G}\) は逆格子ベクトルです。プログラムでは、\(\sum_j q_j \cos(\mathbf{G} \cdot (\mathbf{r}_j - \mathbf{r}_i))\) の部分を
fcal変数で計算しています。自己項 (\(V_{self}\)) 各点電荷に重ねられたガウス関数が、その点電荷自身に与えるポテンシャルを補償するための項です。プログラムでは
UC3に対応します。\[V_{self}(\mathbf{r}_i) = -\frac{2 q_i \cdot \alpha}{\sqrt{\pi}}\]ここで、\(q_i\) はポテンシャルを計算するサイト \(i\) の電荷です。
Madelungポテンシャル \(V_M\) は、これら3つの項の合計 MP = UC1 + UC2 - UC3 として計算されます。
最終的なMadelung定数 \(\alpha_M\) は、計算されたMadelungポテンシャル（eV単位）から以下のように導出されます。
ここで、\(V_M^{\text{eV}}\) はeV単位のMadelungポテンシャル、\(|q_i|\) は対象サイトの電荷の絶対値、\(a\) は結晶の格子定数（プログラムでは lattice_parameters[0]）です。この定義は文献によって異なる場合がありますが、プログラムはこの式を使用しています。
プログラム中で使われる結晶学的な計算（格子ベクトル、体積、距離など）は、カスタムライブラリである tkcrystalbase に含まれる関数によって行われます。
必要な非標準ライブラリとインストール方法
crystal_MP_Ewald.py を実行するために、以下の非標準ライブラリが必要です。
numpy: 数値計算（配列操作、数学関数など）
matplotlib: グラフ描画（ただし、このプログラムでは直接使用されていませんが、インポートされています）
tkcrystalbase: 結晶学計算のためのカスタムライブラリ（格子ベクトル計算、距離計算など）
numpy と matplotlib は pip コマンドでインストールできます。
pip install numpy matplotlib
tkcrystalbase はこのソースコードに付属していないカスタムライブラリであるため、別途入手してPythonの検索パスに配置するか、crystal_MP_Ewald.py と同じディレクトリに配置する必要があります。通常、tkcrystalbase.py というファイル名で提供されます。
必要な入力ファイル
このプログラムは、外部の入力ファイルを必要としません。 結晶の格子パラメータと単位胞内の原子サイトに関する情報（原子名、電荷、位置など）は、プログラムのソースコード内にPythonのリストとNumPy配列としてハードコードされています。
具体的には、以下の変数が初期データとして設定されています。
lattice_parameters: 格子定数（a, b, c, alpha, beta, gamma）のリスト。sites: 各サイトの原子情報（原子名、サイトラベル、原子番号、原子質量、電荷、半径、色、分数座標）のリスト。
生成される出力ファイル
crystal_MP_Ewald.py は、いかなるファイルも生成または保存しません。
全ての計算結果は標準出力（コンソール）に直接表示されます。
出力内容には、格子情報、Ewaldパラメータの詳細、Madelungポテンシャル（JおよびeV単位）、Madelung定数、そして各計算フェーズにかかった時間が含まれます。
コマンドラインでの使用例 (Usage)
プログラムは、Ewaldパラメータ \(\alpha\) (alpha) と計算精度 (prec) をコマンドライン引数として受け取ることができます。引数が提供されない場合、デフォルト値が使用されます。
python crystal_MP_Ewald.py [alpha] [prec]
alpha: Ewaldパラメータ。実数。デフォルト値は0.3。この値は実空間と逆空間の和のバランスを決定します。prec: 計算精度。実数。デフォルト値は1.0e-5。和の打ち切り条件に影響を与えます。
例
デフォルト値で実行する場合：
python crystal_MP_Ewald.py
特定のEwaldパラメータと精度を指定して実行する場合：
python crystal_MP_Ewald.py 0.4 1.0e-6
コマンドラインでの具体的な使用例
ここでは、crystal_MP_Ewald.py をデフォルトのEwaldパラメータ 0.3 と精度 1.0e-5 で実行した場合の出力例を示します。
この例は、NaCl（塩化ナトリウム）結晶（立方晶、格子定数 5.62 Å）の単位胞における \(\mathrm{Na}^{+}\) イオンのMadelungポテンシャルと定数を計算したものです。
実行コマンド:
python crystal_MP_Ewald.py
実行結果の例:
Lattice parameters: [5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
Lattice vectors:
ax: ( 5.62, 0.0, 0.0) A
ay: ( 0.0, 5.62, 0.0) A
az: ( 0.0, 0.0, 5.62) A
Metric tensor:
gij: ( 5.62, 0.0, 0.0) A
( 0.0, 5.62, 0.0) A
( 0.0, 0.0, 5.62) A
Volume: 177.3069 A^3
Unit cell volume: 177.3069 A^3
Reciprocal lattice parameters: [0.1779360853472064, 0.1779360853472064, 0.1779360853472064, 90.0, 90.0, 90.0]
Reciprocal lattice vectors:
Rax: ( 0.1779, 0.0, 0.0) A^-1
Ray: ( 0.0, 0.1779, 0.0) A^-1
Raz: ( 0.0, 0.0, 0.1779) A^-1
Reciprocal lattice metric tensor:
Rgij: ( 0.1779, 0.0, 0.0) A^-1
( 0.0, 0.1779, 0.0) A^-1
( 0.0, 0.0, 0.1779) A^-1
Reciprocal unit cell volume: 0.0056 A^-3
Ewald parameters
alpha: 0.3
precision = 1e-05 = 10^-5
RDmax = 17.6 A, where erfc(alpha*RDmax) = 1.05e-05
nrmax: 4 4 4
cal_N(real): 1801
G2max: 1.834451634842183
exp(-pi^2 * G2max^2 / alpha^2) = 4.998495147871638e-06
hgmax: 5 5 5
cal_N(reciprocal): 13
Time for real space sum : 0.005998134613037109
Time for real reciprocal sum: 0.004999876022338867
Total time : 0.011998891830444336
Madelung potential: -1.24647e-18 J (= -5.35246e-19 + -8.84661e-19 - 1.73431e-19)
Madelung potential: -7.77977 eV (= -3.34079 + -5.52222 - 1.08276)
Madelung constant: -0.87105060
実行結果の説明:
Lattice parameters: 入力された格子定数と、それらから計算された直交座標系での格子ベクトル、計量テンソル、および単位胞の体積が表示されます。
Reciprocal lattice parameters: 逆格子空間における対応する量が計算・表示されます。
Ewald parameters: 使用されたEwaldパラメータ \(\alpha\) と精度
prec、およびそれらに基づいて決定された実空間 (RDmax,nrmax,cal_N(real)) と逆空間 (G2max,hgmax,cal_N(reciprocal)) の計算範囲と項目数が示されます。Time for ... sum: 実空間和と逆空間和、および全体の計算にかかった時間が秒単位で表示されます。
Madelung potential: 計算されたMadelungポテンシャルがジュール（J）と電子ボルト（eV）単位で表示されます。カッコ内は、実空間項、逆空間項、自己項の寄与です。
Madelung constant: 計算されたMadelung定数が表示されます。この値は、プログラム内で定義された式に基づいており、NaClの典型的なMadelung定数 (約 1.74756) に対応する値が出力されます。符号が負であるのは、イオン結晶のMadelungポテンシャルの安定化エネルギーに起因します。この特定のコードの定義では、単位胞の辺長 \(a\) を用いて計算されるため、一般的に最近接距離 \(a_{NN}\) を用いて定義されるMadelung定数の値とは異なることがあります（この例では、格子定数 \(a=5.62\) Å のNaCl立方晶の場合、最近接距離 \(a_{NN} = a/2 = 2.81\) Å。計算結果は \(-0.87105060\) であり、標準的な \(1.74756\) の約半分に符号を反転させた値となっています）。