周期境界条件下における静電相互作用の高速評価：Ewald 法の数理と定式化（理論編）

導入

分子動力学（MD）シミュレーションや第一原理電子状態計算（DFT）において、系内の荷電粒子間に働く静電相互作用（クーロン相互作用）を正確かつ効率的に評価することは、計算コストと精度の双方を左右する決定的な要素である。

結晶や高分子、あるいはバルク溶液の計算では、無限に続くシステムをシミュレートするために周期境界条件（Periodic Boundary Conditions: PBC）が採用される。しかし、長距離相互作用であるクーロン相互作用は距離 $r$ に対して $1/r$ で極めて緩やかに減衰するため、周期境界条件下における実空間の直接和は条件収束（格子ベクトルの和の順序に依存）となり、無秩序な打ち切りは非物理的な結果をもたらす。

本稿では、この周期境界条件下での静電エネルギー評価の難題を数学的・物理的に解決する決定的な手法である Ewald 法（Ewald Summation） の理論的背景と完全な数学的定式化について解説する。実空間と逆格子空間へのクーロン相互作用の厳密な分割、誤差関数を用いた定式化、自己相互作用の補正、さらに電荷中性条件を満たさない系や異方的な周期境界（スラブモデル）に対する双極子補正（Dipole Correction）の数理について、一切の省略を行わずに解説する。

---

1. 周期境界系における静電相互作用と条件収束

1.1 周期境界系における静電エネルギーの定式化

無限に繰り返される3次元周期境界系（基本単位胞の体積を $V$ とする）を考える。単位胞内には $N$ 個の荷電粒子が存在し、各粒子 $i$ は電荷 $q_i$ を持ち、位置 $\mathbf{r}_i$ に配置されている。このとき、基本単位胞とそのすべての周期写像（レプリカ）に含まれる全粒子間の静電相互作用（クーロンエネルギー）の総和を求める。

周期格子ベクトルを $\mathbf{T} = n_1 \mathbf{a}_1 + n_2 \mathbf{a}_2 + n_3 \mathbf{a}_3$ ($n_1, n_2, n_3 \in \mathbb{Z}$) とする。ここで $\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3$ は基本単位胞の基本並進ベクトルである。単位胞1個あたりの全静電エネルギー $E$ は、以下のように定義される：

\[E = \frac{1}{2} \sum_{\mathbf{T}} \sum_{i, j}' \frac{q_i q_j}{\vert \mathbf{r}_{ij} + \mathbf{T} \vert}\]

ここで、$\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$ であり、ダブルサメーション（二重和）のプライム記号（$’$）は、$\mathbf{T} = \mathbf{0}$ のときに自己相互作用項（$i = j$）を除外することを表す。係数 $1/2$ は、ペア相互作用の重複カウントを防ぐためのものである。

1.2 静電ポテンシャルの遠距離減衰と条件収束の数学的起源

実空間におけるクーロン相互作用は距離 $r$ に対して $1/r$ で減衰する。この長距離相互作用の和を無限格子にわたって直接評価しようとすると、数学的な困難に直面する。この困難の本質は、この無限級数が絶対収束（absolute convergence）せず、条件収束（conditional convergence）あるいは発散することにある。

この数学的性質を理解するために、離散的な格子和を連続的な空間積分で近似して評価する。半径 $R$ の球内での静電ポテンシャルまたはエネルギーの積分を考えると、3次元空間の微小体積要素は $d^3\mathbf{r} = 4\pi r^2 dr$ であるため、動径方向の積分は以下の形をとる：

\[\int_0^R \frac{1}{r} d^3\mathbf{r} = \int_0^R \frac{1}{r} \cdot 4\pi r^2 dr = 4\pi \int_0^R r dr = 2\pi R^2\]

この積分結果は、カットオフ半径 $R$ を無限大に飛ばしたとき（$R \to \infty$）、二次式（$\propto R^2$）で発散する。各格子点からの寄与が減衰する速さよりも、シェル（球殻）に含まれる格子点の数が距離の2乗（$\propto r^2$）で増加する速さの方が勝るためである。したがって、電荷の具体的な打ち消し合い（charge cancellation）がない限り、この系全体の静電エネルギーは容易に発散する。

さらに重要なのは、電荷中性系（$\sum_i q_i = 0$）であっても、この級数和が絶対収束しない（すなわち、各項の絶対値の和が有限値に収束しない）という点である。リーマンの再配列定理（Riemann rearrangement theorem）によれば、条件収束する無限級数は、その和をとる順序（極限操作の順序）を変更することによって、任意の実数値に収束させたり、あるいは発散させたりすることができる。

1.3 和の順序（境界形状）依存性：球状境界 vs 立方体状境界

この和の順序依存性を物理的に説明する。無限格子和を計算する際、単位胞の「集まり」をどのような形状で無限大に拡張していくかによって、得られる静電エネルギーの値が変化する。代表的な例として、球状に領域を拡張していく場合（球状境界）と、直方体・立方体状に拡張していく場合（立方体状境界）を比較する。

1. 球状の拡大シェル（Spherical expanding shell）: 境界を半径 $R$ の球（すなわち $\vert \mathbf{T} \vert \le R$ を満たす格子ベクトル $\mathbf{T}$ の集合）とし、$R \to \infty$ とする。単位胞が実効的な双極子モーメント $\mathbf{M} = \sum_i q_i \mathbf{r}i$ を持つ場合、この球状境界の表面に分極電荷（表面電荷）が誘起される。この表面電荷は球の内部に一様な脱分極場（depolarizing field） $\mathbf{E}{\text{dep}} = -\frac{4\pi}{3V} \mathbf{M}$ を生じさせる。この脱分極場と単位胞内の双極子モーメントとの相互作用により、以下の補正項（表面エネルギー）がエネルギーに加わる：

   \[E_{\text{surf}}^{\text{sphere}} = \frac{2\pi}{3V} \vert \mathbf{M} \vert^2\]

2. 立方体状の拡大シェル（Cubic expanding shell）: 境界を各辺の長さが $2L$ の立方体 $C_L = { \mathbf{T} \ \text{s.t.} \ -L \le T_x, T_y, T_z \le L }$ とし、$L \to \infty$ とする。この場合、立方体の表面に誘起される電荷分布は球状の場合と異なり、球対称な脱分極場とは異なる電場が発生する。結果として、得られる静電エネルギーの値は球状境界の場合からシフトする。

このように、無限和の順序（マクロな系境界の幾何学的形状）は、系のマクロな静電分極状態を通じてエネルギー値に物理的な差異をもたらす。したがって、周期境界条件下での静電エネルギーを一意に定義するためには、背景の境界条件（周囲の媒質の誘電率など）を数学的に指定する必要がある。通常、Ewald法では、球状の拡大シェルを仮定し、さらに周囲を導体（誘電率 $\epsilon_{\infty} = \infty$、錫箔境界条件 / tinfoil boundary condition）で取り囲むことで脱分極場を完全に相殺し、表面エネルギー項をゼロとする基準を採用することが一般的である。

1.4 電荷中性条件の絶対的必要性

上述の通り、周期系において静電エネルギーが有限値として確定するための前提として、電荷中性条件（charge neutrality condition）：

\[\sum_{i=1}^N q_i = 0\]

が満たされていることが不可欠である。

もし単位胞内の全電荷の総和がゼロでない（$\sum_i q_i = Q \neq 0$、すなわち系全体が正または負に帯電している）場合、エネルギーの定義式におけるモノポール-モノポール相互作用（最も長距離に及ぶ寄与）は以下のように振る舞う：

\[E_{\text{mono}} \approx \frac{1}{2} \sum_{\mathbf{T} \neq \mathbf{0}} \frac{Q^2}{\vert \mathbf{T} \vert}\]

この和を積分で近似すると、その積分は上述したように $\propto R^2$ であり、無限大への極限で確実に二次発散する。物理的には、無限に並んだ同種の電荷が互いに反発し合い、系全体の静電電位が無限大に発散することを意味する。

実際の第一原理計算や古典分子動力学において、電子の脱離や不純物の導入などによって過剰電荷を持つ系（チャージドシステム）を扱う場合でも、バックグラウンドに均一な反対電荷（uniform neutralizing background charge）を仮想的に配置することで、数学的な電荷中性を強制的に満たし、この発散を回避する手法が取られる。

---

2. クーロン相互作用の分割（Ewald 分割）と物理的イメージ

2.1 クーロンポテンシャルの恒等式

Ewald 法の核心は、長距離に及ぶ緩やかな減衰（実空間で条件収束）と短距離での特異性（$r \to 0$ での発散）という二面性を持つクーロンポテンシャル $1/r$ を、誤差関数（error function, $\text{erf}$）とその相補誤差関数（complementary error function, $\text{erfc}$）を用いて、数学的に厳密に2つの寄与に分割することである。

まず、以下の誤差関数の定義から出発する：

\[\text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt\] \[\text{erfc}(x) = 1 - \text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_x^{\infty} e^{-t^2} dt\]

定義より明らかに、任意の正の実数 $\eta > 0$ に対して次の恒等式が成立する：

\[\text{erf}(\eta r) + \text{erfc}(\eta r) = 1\]

この両辺を $r$ で除すことで、クーロンポテンシャル $1/r$ に関する極めて重要な恒等式（Ewald 分割の基礎）が得られる：

\[\frac{1}{r} = \frac{\text{erfc}(\eta r)}{r} + \frac{\text{erf}(\eta r)}{r}\]

ここで、パラメータ $\eta$ はEwaldパラメータ（または遮蔽パラメータ）と呼ばれ、実空間と逆格子空間のどちらに計算の重きを置くかを調整する正の定数である。この恒等式は数学的に厳密であり、どのような $\eta > 0$ に対しても成立する。

2.2 遮蔽電荷による物理的描像

この数学的分割は、極めて明確な物理的描像と対応している。

電荷 $q_i$ を持つ各点電荷（位置 $\mathbf{r}_i$）について、その周囲に点電荷とは逆符号の電荷 $-q_i$ を持つ、広がりを持った電荷雲（遮蔽電荷雲, screening cloud）を仮想的に重ね合わせることを考える。Ewald法では、この遮蔽電荷雲の形状として以下の3次元ガウス分布を採用する：

\[\rho_{\text{screen}, i}(\mathbf{r}) = -q_i \left( \frac{\eta^2}{\pi} \right)^{3/2} \exp\left( -\eta^2 \vert \mathbf{r} - \mathbf{r}_i \vert^2 \right)\]

このガウス分布による総電荷量を空間全体で積分すると、次のように正確に $-q_i$ となる：

\[\int \rho_{\text{screen}, i}(\mathbf{r}) d^3\mathbf{r} = -q_i \left( \frac{\eta^2}{\pi} \right)^{3/2} \int_0^{\infty} \exp\left( -\eta^2 r^2 \right) \cdot 4\pi r^2 dr = -q_i\]

この逆符号の電荷雲の存在により、もとの点電荷 $q_i$ の作るポテンシャルは近傍では点電荷そのものの寄与が残るものの、距離が離れるにつれてガウス電荷雲によって電気的に遮蔽され、急速（指数関数的）にゼロへと減衰するようになる。この「遮蔽された点電荷」が作る有効的なポテンシャルが、恒等式の第1項である短距離項（実空間項）に対応する：

\[\phi_{\text{short}, i}(\mathbf{r}) = \frac{q_i \text{erfc}(\eta \vert \mathbf{r} - \mathbf{r}_i \vert)}{\vert \mathbf{r} - \mathbf{r}_i \vert}\]

しかし、実際の系にはこのような逆符号の遮蔽電荷雲は存在しない。したがって、数学的・物理的整合性を保つためには、挿入した遮蔽電荷雲を完全に打ち消す（相殺する）ための、点電荷と同符号のガウス電荷雲（相殺電荷雲, neutralizing/compensating cloud）を同時に配置しなければならない。この相殺電荷雲の電荷密度は以下で与えられる：

\[\rho_{\text{comp}, i}(\mathbf{r}) = +q_i \left( \frac{\eta^2}{\pi} \right)^{3/2} \exp\left( -\eta^2 \vert \mathbf{r} - \mathbf{r}_i \vert^2 \right)\]

この同符号のガウス電荷雲が作るスムーズな（特異性のない）ポテンシャルが、恒等式の第2項である長距離項（逆格子空間項）に対応する：

\[\phi_{\text{long}, i}(\mathbf{r}) = \frac{q_i \text{erf}(\eta \vert \mathbf{r} - \mathbf{r}_i \vert)}{\vert \mathbf{r} - \mathbf{r}_i \vert}\]

2.3 実空間・逆格子空間への役割分担

以上の物理的描像に基づき、全静電エネルギーの計算は以下の2つの異なる空間における和へと役割が分担される：

1. 実空間和（Real-space sum）: 各点電荷 $q_i$ と、他の「遮蔽された」電荷 $q_j$ との間の相互作用を計算する。この相互作用ポテンシャル $\phi_{\text{short}}$ は距離に対して指数関数的に急激に減衰するため、実空間での打ち切り（カットオフ）が完全に正当化され、絶対収束する短い範囲の和として高速かつ高精度に計算できる。

2. 逆格子空間和（Reciprocal-space sum / Fourier-space sum）: 相殺電荷雲 $\rho_{\text{comp}, i}(\mathbf{r})$ の全空間における重ね合わせが作るポテンシャルと、系全体の点電荷との相互作用を計算する。この電荷分布は周期境界条件下で連続的かつ滑らか（smooth）な関数であるため、フーリエ変換と非常に相性が良い。滑らかな周期関数のフーリエ級数は逆格子空間において非常に速く収束するため、少数の波数（逆格子ベクトル）の和として効率的かつ厳密に評価可能である。

この実空間と逆格子空間の協調により、直接評価では条件収束かつ $O(N^2)$ であった静電エネルギーの計算が、厳密性を保ったまま高効率に計算可能となる。

---

3. 実空間和の定式化と収束性解析

3.1 実空間エネルギー項の定式化

前節の Ewald 分割に基づき、周期境界系における全静電エネルギー $E$ のうち、実空間の寄与である $E_{\text{real}}$ は、各点電荷と遮蔽された電荷との相互作用の総和として以下のように定式化される：

\[E_{\text{real}} = \frac{1}{2} \sum_{\mathbf{T}} \sum_{i, j}' q_i q_j \frac{\text{erfc}(\eta \vert \mathbf{r}_{ij} + \mathbf{T} \vert)}{\vert \mathbf{r}_{ij} + \mathbf{T} \vert}\]

ここで、$\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$ であり、$\mathbf{T}$ は周期格子ベクトルである。プライム記号（$’$）は、同一単位胞内（$\mathbf{T} = \mathbf{0}$）における自己相互作用項（$i = j$）を除外することを意味する。この実空間エネルギーは、相補誤差関数 $\text{erfc}(x)$ の性質により、格子和の極限をとる順番に依存せず一意に絶対収束する。

3.2 相補誤差関数の数学的性質と漸近展開

実空間和が急速に収束する数学的根拠は、相補誤差関数 $\text{erfc}(x)$ の漸近的な振る舞いにある。$\text{erfc}(x)$ は $x \ge 0$ において単調減少する関数であり、$x \to 0$ で $\text{erfc}(0) = 1$、$x \to \infty$ で $\text{erfc}(x) \to 0$ となる。

十分大きな引数 $x \gg 1$ に対して、相補誤差関数は以下の漸近展開（asymptotic expansion）を持つ：

\[\text{erfc}(x) = \frac{e^{-x^2}}{x\sqrt{\pi}} \left( 1 - \frac{1}{2x^2} + \frac{3}{4x^4} - \frac{15}{8x^6} + \mathcal{O}(x^{-8}) \right)\]

この展開式が示すように、十分に離れた距離 $r = \vert \mathbf{r}_{ij} + \mathbf{T} \vert$ において、実空間ポテンシャル項は以下のように近似される：

\[\frac{\text{erfc}(\eta r)}{r} \approx \frac{e^{-\eta^2 r^2}}{\eta \sqrt{\pi} r^2}\]

すなわち、通常のクーロンポテンシャル $1/r$ が距離の逆数という非常に緩やかな減衰を示すのに対し、Ewald 分割後の実空間ポテンシャルは、ガウス分布の指数減衰因子 $e^{-\eta^2 r^2}$ によって支配され、きわめて急速にゼロに漸近する。これが、実空間和における直接的なカットオフ（打ち切り）の数学的正当性の根拠である。

3.3 カットオフ基準と Ewald パラメータ $\eta$ の選定

実空間での高速な評価を実現するため、ある有限のカットオフ半径 $R_c$ を導入し、$\vert \mathbf{r}_{ij} + \mathbf{T} \vert \ge R_c$ となる格子ベクトルおよび粒子ペアからの寄与を無視する近似を行う。

このとき生じるエネルギーの切り捨て誤差（truncation error）を許容精度 $\delta$（たとえば $\delta = 10^{-5} \sim 10^{-6}$）以下に抑えたい。エネルギー寄与の絶対値が最大となる上限を見積もるため、最も厳しい基準として以下の条件を課す：

\[\frac{\text{erfc}(\eta R_c)}{R_c} \le \delta\]

あるいは、簡略化のために分子の相補誤差関数の減衰のみに着目して、以下の基準を用いることもある：

\[\text{erfc}(\eta R_c) \le \delta\]

ここで、十分大きな引数に対する近似 $\text{erfc}(x) \approx e^{-x^2}$ を用いると、上式は以下のように書き換えられる：

\[e^{-\eta^2 R_c^2} \le \delta \implies \eta^2 R_c^2 \ge -\ln \delta \implies \eta \ge \frac{\sqrt{- \ln \delta}}{R_c}\]

この関係式は、実空間カットオフ半径 $R_c$ と Ewald パラメータ $\eta$ の間にトレードオフが存在することを示している：

• $\eta$ を大きく設定すると、実空間ポテンシャルはより急激に減衰するため、より小さなカットオフ半径 $R_c$ で十分な精度が得られ、実空間の計算コストは削減される。しかしこの場合、相殺電荷雲の分布が急峻（デルタ関数に近く）になり、逆格子空間での収束が遅くなるため、逆格子空間の計算コスト（必要な波数の数）が増大する。
• 逆に $\eta$ を小さく設定すると、逆格子空間での収束は早まるが、実空間ポテンシャルの減衰が緩やかになり、大きな $R_c$ が必要となって実空間の計算コストが増大する。

通常、古典 MD や DFT 計算パッケージにおいては、単位胞の最小の辺の長さ $L_{\text{min}}$ に対して、実空間計算において隣接するセル（レプリカ）の二重カウントを避けるために最小基本並進写像基準（minimum image convention）：

\[R_c < \frac{L_{\text{min}}}{2}\]

を課すことが多い。この制約下で、要求精度 $\delta$ から Ewald パラメータ $\eta$ の値が決定され、それに伴って逆格子空間でのカットオフ波数 $K_c$ が決定される。実空間和の計算量 $O(N^2)$ と逆格子空間和の計算量（Ewald 法では $O(N^{3/2})$、粒子メッシュ Ewald（PME）法では $O(N \log N)$）のバランスが最適となるよう、$\eta$ が自動的にチューニングされる。

---

4. 逆格子空間和の定式化と自己相互作用補正

4.1 ガウス型電荷分布に対するポアソン方程式のフーリエ空間での解法

Ewald 分割における逆格子空間（フーリエ空間）の計算対象は、点電荷と同符号の広がりを持ったガウス型電荷分布（相殺電荷雲） $\rho_{\text{comp}}(\mathbf{r})$ の集合が作る静電ポテンシャル $\phi_{\text{comp}}(\mathbf{r})$ と、点電荷群との相互作用エネルギーである。

まず、単一のガウス電荷雲による電荷密度分布を考える。位置 $\mathbf{r}_i$ に中心を持つ電荷 $q_i$ のガウス電荷雲およびその周期境界条件下の全レプリカによる電荷密度分布 $\rho_i(\mathbf{r})$ は、次のように定義される：

\[\rho_i(\mathbf{r}) = \sum_{\mathbf{T}} q_i \left( \frac{\eta^2}{\pi} \right)^{3/2} \exp\left( -\eta^2 \vert \mathbf{r} - \mathbf{r}_i - \mathbf{T} \vert^2 \right)\]

この電荷密度は格子ベクトル $\mathbf{T}$ に対して周期的な関数であるため、逆格子ベクトル $\mathbf{G}$ を用いて次のようにフーリエ級数展開することができる：

\[\rho_i(\mathbf{r}) = \sum_{\mathbf{G}} \tilde{\rho}_i(\mathbf{G}) \exp(i \mathbf{G} \cdot \mathbf{r})\]

ここで、フーリエ成分 $\tilde{\rho}_i(\mathbf{G})$ は、基本単位胞の体積を $V$ とすると、単位胞内での積分として以下のように与えられる：

\[\tilde{\rho}_i(\mathbf{G}) = \frac{1}{V} \int_V \rho_i(\mathbf{r}) \exp(-i \mathbf{G} \cdot \mathbf{r}) d^3\mathbf{r}\]

上式に $\rho_i(\mathbf{r})$ の定義式を代入する：

\[\tilde{\rho}_i(\mathbf{G}) = \frac{1}{V} \int_V \left[ \sum_{\mathbf{T}} q_i \left( \frac{\eta^2}{\pi} \right)^{3/2} \exp\left( -\eta^2 \vert \mathbf{r} - \mathbf{r}_i - \mathbf{T} \vert^2 \right) \right] \exp(-i \mathbf{G} \cdot \mathbf{r}) d^3\mathbf{r}\]

ここで、単位胞 $V$ 内での積分とすべての周期格子ベクトル $\mathbf{T}$ についての総和の順序を交換し、$\mathbf{r}’ = \mathbf{r} - \mathbf{T}$ という変数変換を行う。この変換により、単位胞 $V$ 内の積分を全空間 $\mathbb{R}^3$ における積分へと統合することができる：

\[\tilde{\rho}_i(\mathbf{G}) = \frac{q_i}{V} \left( \frac{\eta^2}{\pi} \right)^{3/2} \int_{\mathbb{R}^3} \exp\left( -\eta^2 \vert \mathbf{r}' - \mathbf{r}_i \vert^2 \right) \exp(-i \mathbf{G} \cdot (\mathbf{r}' + \mathbf{T})) d^3\mathbf{r}'\]

ここで、任意の逆格子ベクトル $\mathbf{G}$ と任意の周期格子ベクトル $\mathbf{T}$ に対して $\exp(-i \mathbf{G} \cdot \mathbf{T}) = 1$ が成立するため、上式は次のように簡略化される：

\[\tilde{\rho}_i(\mathbf{G}) = \frac{q_i}{V} \left( \frac{\eta^2}{\pi} \right)^{3/2} \int_{\mathbb{R}^3} \exp\left( -\eta^2 \vert \mathbf{r}' - \mathbf{r}_i \vert^2 \right) \exp(-i \mathbf{G} \cdot \mathbf{r}') d^3\mathbf{r}'\]

さらに、$\mathbf{s} = \mathbf{r}’ - \mathbf{r}_i$ と置くと、$\mathbf{r}’ = \mathbf{s} + \mathbf{r}_i$ であり、$d^3\mathbf{r}’ = d^3\mathbf{s}$ となる。これらを代入すると：

\[\tilde{\rho}_i(\mathbf{G}) = \frac{q_i}{V} \left( \frac{\eta^2}{\pi} \right)^{3/2} \exp(-i \mathbf{G} \cdot \mathbf{r}_i) \int_{\mathbb{R}^3} \exp\left( -\eta^2 s^2 \right) \exp(-i \mathbf{G} \cdot \mathbf{s}) d^3\mathbf{s}\]

被積分関数はデカルト座標の各軸（$x, y, z$）について独立に変数分離できる。すなわち、3次元の積分は以下の3つの1次元ガウス積分の積として表される：

\[\int_{\mathbb{R}^3} \exp\left( -\eta^2 s^2 \right) \exp(-i \mathbf{G} \cdot \mathbf{s}) d^3\mathbf{s} = \prod_{\alpha=x,y,z} \int_{-\infty}^{\infty} \exp\left( -\eta^2 s_{\alpha}^2 \right) \exp(-i G_{\alpha} s_{\alpha}) ds_{\alpha}\]

各デカルト成分 $\alpha$ について平方完成を行うと、次のように変形できる：

\[-\eta^2 s_{\alpha}^2 - i G_{\alpha} s_{\alpha} = -\eta^2 \left( s_{\alpha} + \frac{i G_{\alpha}}{2\eta^2} \right)^2 - \frac{G_{\alpha}^2}{4\eta^2}\]

これより、1次元のガウス積分は以下のようになる：

\[\int_{-\infty}^{\infty} \exp\left( -\eta^2 s_{\alpha}^2 - i G_{\alpha} s_{\alpha} \right) ds_{\alpha} = \exp\left( -\frac{G_{\alpha}^2}{4\eta^2} \right) \int_{-\infty}^{\infty} \exp\left( -\eta^2 \left( s_{\alpha} + \frac{i G_{\alpha}}{2\eta^2} \right)^2 \right) ds_{\alpha}\]

複素平面における積分経路のシフトを考慮すると、右辺の積分値は実軸上のガウス積分 $\int_{-\infty}^{\infty} \exp(-\eta^2 x^2) dx = \frac{\sqrt{\pi}}{\eta}$ に等しい。したがって、

\[\int_{-\infty}^{\infty} \exp\left( -\eta^2 s_{\alpha}^2 - i G_{\alpha} s_{\alpha} \right) ds_{\alpha} = \frac{\sqrt{\pi}}{\eta} \exp\left( -\frac{G_{\alpha}^2}{4\eta^2} \right)\]

この3成分の積をとることで、3次元ガウス積分の結果が得られる：

\[\int_{\mathbb{R}^3} \exp\left( -\eta^2 s^2 \right) \exp(-i \mathbf{G} \cdot \mathbf{s}) d^3\mathbf{s} = \left( \frac{\sqrt{\pi}}{\eta} \right)^3 \exp\left( -\frac{G^2}{4\eta^2} \right) = \left( \frac{\pi}{\eta^2} \right)^{3/2} \exp\left( -\frac{G^2}{4\eta^2} \right)\]

ここで $G^2 = G_x^2 + G_y^2 + G_z^2$ である。この結果を $\tilde{\rho}_i(\mathbf{G})$ の式に代入すると、係数部分が相殺され、次の簡潔な表現が得られる：

\[\tilde{\rho}_i(\mathbf{G}) = \frac{q_i}{V} \exp(-i \mathbf{G} \cdot \mathbf{r}_i) \exp\left( -\frac{G^2}{4\eta^2} \right)\]

もしガウス電荷雲が原点に中心を持つ場合（すなわち $\mathbf{r}_i = \mathbf{0}$）、そのフーリエ成分は次のようになる：

\[\tilde{\rho}(\mathbf{G}) = \frac{q_i}{V} \exp\left( -\frac{G^2}{4\eta^2} \right)\]

4.2 構造因子を用いた逆格子空間エネルギーの導出

系全体のガウス電荷分布 $\rho_{\text{comp}}(\mathbf{r})$ は、すべての粒子 $i$ によるガウス電荷雲の重ね合わせである：

\[\rho_{\text{comp}}(\mathbf{r}) = \sum_{i=1}^N \rho_i(\mathbf{r})\]

この全電荷密度のフーリエ成分 $\tilde{\rho}_{\text{comp}}(\mathbf{G})$ は、各成分の線形結合となる：

\[\tilde{\rho}_{\text{comp}}(\mathbf{G}) = \sum_{i=1}^N \tilde{\rho}_i(\mathbf{G}) = \frac{1}{V} \exp\left( -\frac{G^2}{4\eta^2} \right) \sum_{i=1}^N q_i \exp(-i \mathbf{G} \cdot \mathbf{r}_i)\]

ここで、逆格子空間における電荷分布の位相的な重ね合わせを表す物理量として、構造因子（structure factor） $S(\mathbf{G})$ を以下のように定義する：

\[S(\mathbf{G}) = \sum_{i=1}^N q_i \exp(i \mathbf{G} \cdot \mathbf{r}_i)\]

構造因子の複素共役 $S^(\mathbf{G})$ は $S^(\mathbf{G}) = \sum_i q_i \exp(-i \mathbf{G} \cdot \mathbf{r}_i)$ であるため、全電荷密度のフーリエ成分は次のように書ける：

\[\tilde{\rho}_{\text{comp}}(\mathbf{G}) = \frac{1}{V} \exp\left( -\frac{G^2}{4\eta^2} \right) S^*(\mathbf{G})\]

次に、この電荷分布が形成する静電ポテンシャル $\phi_{\text{comp}}(\mathbf{r})$ を求めるため、ポアソン方程式を解く：

\[\nabla^2 \phi_{\text{comp}}(\mathbf{r}) = -4\pi \rho_{\text{comp}}(\mathbf{r})\]

静電ポテンシャルも電荷密度と同様にフーリエ級数展開する：

\[\phi_{\text{comp}}(\mathbf{r}) = \sum_{\mathbf{G}} \tilde{\phi}_{\text{comp}}(\mathbf{G}) \exp(i \mathbf{G} \cdot \mathbf{r})\]

これらをポアソン方程式に代入すると、ラプラシアンの作用 $\nabla^2 \exp(i \mathbf{G} \cdot \mathbf{r}) = -G^2 \exp(i \mathbf{G} \cdot \mathbf{r})$ により、各 $\mathbf{G} \neq \mathbf{0}$ に対して以下の代数方程式が得られる：

\[-G^2 \tilde{\phi}_{\text{comp}}(\mathbf{G}) = -4\pi \tilde{\rho}_{\text{comp}}(\mathbf{G})\]

したがって、$\mathbf{G} \neq \mathbf{0}$ におけるポテンシャルのフーリエ成分は次のようになる：

\[\tilde{\phi}_{\text{comp}}(\mathbf{G}) = \frac{4\pi}{G^2} \tilde{\rho}_{\text{comp}}(\mathbf{G}) = \frac{4\pi}{V G^2} \exp\left( -\frac{G^2}{4\eta^2} \right) S^*(\mathbf{G})\]

逆格子空間における静電エネルギー $E_{\text{recip}}$ は、系内の全点電荷群とポテンシャル $\phi_{\text{comp}}(\mathbf{r})$ とのペア相互作用エネルギーとして定義される：

\[E_{\text{recip}} = \frac{1}{2} \sum_{j=1}^N q_j \phi_{\text{comp}}(\mathbf{r}_j)\]

ポテンシャルのフーリエ級数展開を代入すると：

\[E_{\text{recip}} = \frac{1}{2} \sum_{j=1}^N q_j \sum_{\mathbf{G} \neq \mathbf{0}} \tilde{\phi}_{\text{comp}}(\mathbf{G}) \exp(i \mathbf{G} \cdot \mathbf{r}_j)\]

ここで $\mathbf{G} = \mathbf{0}$ 項は系の全体的な中性バックグラウンドに対応し、周期境界系における平均電位の基準点（無限遠での基準）を考慮して除外する。$\tilde{\phi}_{\text{comp}}(\mathbf{G})$ の具体的な表式を代入して整理する：

\[E_{\text{recip}} = \frac{1}{2} \sum_{j=1}^N q_j \sum_{\mathbf{G} \neq \mathbf{0}} \left[ \frac{4\pi}{V G^2} \exp\left( -\frac{G^2}{4\eta^2} \right) S^*(\mathbf{G}) \right] \exp(i \mathbf{G} \cdot \mathbf{r}_j)\]

和の順序を入れ替えると：

\[E_{\text{recip}} = \frac{2\pi}{V} \sum_{\mathbf{G} \neq \mathbf{0}} \frac{1}{G^2} \exp\left( -\frac{G^2}{4\eta^2} \right) S^*(\mathbf{G}) \left[ \sum_{j=1}^N q_j \exp(i \mathbf{G} \cdot \mathbf{r}_j) \right]\]

括弧内の和はまさに構造因子 $S(\mathbf{G})$ そのものである。したがって、$S^*(\mathbf{G}) S(\mathbf{G}) = \vert S(\mathbf{G}) \vert^2$ を用いると、逆格子空間エネルギーは次のように書き下される：

\[E_{\text{recip}} = \frac{2\pi}{V} \sum_{\mathbf{G} \neq \mathbf{0}} \frac{1}{G^2} \exp\left( -\frac{G^2}{4\eta^2} \right) \vert S(\mathbf{G}) \vert^2\]

なお、数式表現の便宜上、前定数を $\frac{4\pi}{V}$ に統一し、総和の内部に因子 $1/2$ を含める形で定式化されることも多い：

\[E_{\text{recip}} = \frac{4\pi}{V} \sum_{\mathbf{G} \ne 0} \frac{1}{2G^2} \exp\left( -\frac{G^2}{4\eta^2} \right) \vert S(\mathbf{G}) \vert^2\]

なお、一部の文献やシミュレーションパッケージのソースコードでは、前定数として $\frac{4\pi}{V}$ を用いて記述されている場合があるが、これは逆格子空間の総和を半空間（位相的な重複を除いた一意なベクトル群）に制限して重複カウントを排除しているか、あるいは一粒子あたりの静電電位を定義しているという前提があるため、数理モデルを実装する際には総和の範囲の定義に細心の注意を払う必要がある。

4.3 自己相互作用補正の厳密な導出

逆格子空間でのポテンシャル $\phi_{\text{comp}}(\mathbf{r})$ は、すべての粒子が作るガウス電荷雲を重ね合わせた連続的な電荷密度分布から求められている。このため、逆格子エネルギー $E_{\text{recip}} = \frac{1}{2} \sum_j q_j \phi_{\text{comp}}(\mathbf{r}j)$ には、各点電荷 $q_i$ が自分自身に付随する相殺ガウス電荷雲 $\rho{\text{comp}, i}$ が作る電位と相互作用する「自己相互作用エネルギー」（spurious self-interaction energy）が含まれてしまっている。

物理的な静電エネルギーには、粒子自身と自身の電荷雲との自己相互作用は存在しないため、この寄与を厳密に計算して全エネルギーから差し引かなければならない。

位置 $\mathbf{r}i$ に置かれた点電荷 $q_i$ と、それに付随するガウス電荷雲（電荷総量 $+q_i$、幅パラメータ $\eta$）との自己相互作用エネルギー $E{\text{self}, i}$ を求める。 ガウス電荷雲（原点中心）が作る静電ポテンシャル $\phi_{\text{Gaussian}}(\mathbf{r})$ は、ポアソン方程式を球対称条件のもとで実空間で解く（あるいは逆フーリエ変換を行う）ことで、以下のように相補誤差関数を用いて表される：

\[\phi_{\text{Gaussian}}(\mathbf{r}) = \frac{q_i \text{erf}(\eta r)}{r}\]

このガウス電荷雲がその中心（$r \to 0$）に作る電位 $\phi_{\text{Gaussian}}(\mathbf{0})$ は、誤差関数のテイラー展開 $\text{erf}(x) = \frac{2}{\sqrt{\pi}} \left( x - \frac{x^3}{3} + \dots \right)$ を用いて極限をとることで求められる：

\[\phi_{\text{Gaussian}}(\mathbf{0}) = \lim_{r \to 0} q_i \frac{\text{erf}(\eta r)}{r} = q_i \lim_{r \to 0} \frac{1}{r} \left[ \frac{2}{\sqrt{\pi}} \left( \eta r - \frac{(\eta r)^3}{3} + \mathcal{O}(r^5) \right) \right] = q_i \frac{2\eta}{\sqrt{\pi}}\]

したがって、この自己電位と粒子 $i$ の点電荷 $q_i$ 自体との相互作用エネルギー $E_{\text{self}, i}$ は、係数 $1/2$ を伴って次のようになる：

\[E_{\text{self}, i} = \frac{1}{2} q_i \phi_{\text{Gaussian}}(\mathbf{0}) = \frac{\eta}{\sqrt{\pi}} q_i^2\]

これをすべての粒子について足し合わせることで、差し引くべき全自己相互作用エネルギー（自己相互作用補正項） $E_{\text{self}}$ が得られる：

\[E_{\text{self}} = \frac{\eta}{\sqrt{\pi}} \sum_{i=1}^N q_i^2\]

極限操作による定式化との対応

この自己相互作用補正の物理的意味を、点電荷が自身の実効的な静電電位から無限大の発散（裸の Coulomb 自己エネルギー $q_i^2/r$）を差し引いて正則化する極限操作としても解釈できる。すなわち、点電荷の位置 $\mathbf{r} \to \mathbf{0}$ において、ガウス電荷雲によるポテンシャルから点電荷自身の特異的なクーロンポテンシャル $q_i/r$ を引いた差分の極限として、次のように表される：

\[E_{\text{self}, i} = \frac{1}{2} \lim_{r \to 0} \left[ q_i \phi_{\text{Gaussian}}(\mathbf{r}) - \frac{q_i^2}{r} \right]\]

実際に括弧内の数式を整理すると、次のようになる：

\[q_i \phi_{\text{Gaussian}}(\mathbf{r}) - \frac{q_i^2}{r} = q_i \left( \frac{q_i \text{erf}(\eta r)}{r} \right) - \frac{q_i^2}{r} = \frac{q_i^2}{r} \left( \text{erf}(\eta r) - 1 \right) = -q_i^2 \frac{\text{erfc}(\eta r)}{r}\]

この $r \to 0$ における極限は、ロピタルの定理または相補誤差関数の展開式を用いて評価できる：

\[\lim_{r \to 0} \left[ -q_i^2 \frac{\text{erfc}(\eta r)}{r} \right] = -q_i^2 \lim_{r \to 0} \left( \frac{- \frac{2\eta}{\sqrt{\pi}} e^{-\eta^2 r^2}}{1} \right) = \frac{2\eta}{\sqrt{\pi}} q_i^2\]

この正則化された値に、ペア相互作用の重複を防ぐための係数 $1/2$ を乗じる（あるいは、符号を含めた補正エネルギー項としての差し引き分を定義する）ことで、差し引くべき正のエネルギー項として次の表現が厳密に得られる：

\[E_{\text{self}} = \frac{\eta}{\sqrt{\pi}} \sum_{i=1}^N q_i^2\]

このようにして、Ewald 法における逆格子空間和から生じる unphysical な自己相互作用は完全に相殺され、有限で物理的な静電エネルギーが保証される。

---

5. バックグラウンド電荷と双極子補正

5.1 $\mathbf{G}=\mathbf{0}$ 項の除外と一様中性背景電荷（Jellium）の物理的意味

逆格子空間における静電エネルギー $E_{\text{recip}}$ の定式化において、波数ベクトル $\mathbf{G}=\mathbf{0}$ の項を除外した。この数学的処理の物理的意味と、それが系全体に与える影響について詳細に議論する。

単位胞内の全電荷の総和（正味の電荷）を $Q = \sum_{i=1}^N q_i$ とする。系が電荷中性でない場合（$Q \neq 0$）、あるいは Ewald 分割において各点電荷に付随する相殺ガウス電荷雲の総和を考える場合、全空間における平均電荷密度はゼロにならない。このとき、ポアソン方程式を逆格子空間で解く際、$\mathbf{G}=\mathbf{0}$ における電荷密度のフーリエ成分 $\tilde{\rho}(\mathbf{0})$ は以下のようになる：

\[\tilde{\rho}(\mathbf{0}) = \frac{1}{V} \int_V \rho(\mathbf{r}) d^3\mathbf{r} = \frac{Q}{V}\]

もし $\mathbf{G}=\mathbf{0}$ の寄与をそのままポアソン方程式 $-G^2 \tilde{\phi}(\mathbf{G}) = -4\pi \tilde{\rho}(\mathbf{G})$ に代入しようとすると、左辺がゼロになるのに対し、右辺は $Q \neq 0$ の場合に非ゼロの定数（$-4\pi Q/V$）となり、数学的矛盾（発散）が生じる。

物理的には、これは電荷密度 $\rho(\mathbf{r})$ が空間全体で平均値 $Q/V$ を持つため、無限に広い系全体の静電ポテンシャルが二次関数的に発散することに対応する。

この発散を回避するために、Ewald 法では $\mathbf{G}=\mathbf{0}$ の項を明示的に除外して静電ポテンシャルを計算する。この「$\mathbf{G}=\mathbf{0}$ の除外」という数学的操作は、物理的には系全体に一様で平坦な中性化背景電荷（uniform neutralizing background charge、いわゆる Jellium 背景）を導入することと等価である。

具体的には、元の電荷分布 $\rho(\mathbf{r})$ に対し、密度 $\rho_{\text{jellium}} = -Q/V$ の一様な背景電荷を重ね合わせた新しい周期電荷分布 $\rho_{\text{neutral}}(\mathbf{r})$ を定義する：

\[\rho_{\text{neutral}}(\mathbf{r}) = \rho(\mathbf{r}) - \frac{Q}{V}\]

この中性化された電荷分布のフーリエ成分 $\tilde{\rho}_{\text{neutral}}(\mathbf{G})$ は、$\mathbf{G} \neq \mathbf{0}$ では元の $\tilde{\rho}(\mathbf{G})$ と完全に一致し、$\mathbf{G} = \mathbf{0}$ では以下のようになる：

\[\tilde{\rho}_{\text{neutral}}(\mathbf{0}) = \frac{1}{V} \int_V \left( \rho(\mathbf{r}) - \frac{Q}{V} \right) d^3\mathbf{r} = \frac{Q}{V} - \frac{Q}{V} = 0\]

したがって、$\mathbf{G} = \mathbf{0}$ における電荷密度が厳密にゼロとなり、ポアソン方程式の特異性が解消される。結果として、$\mathbf{G}=\mathbf{0}$ 項を除外して計算された逆格子ポテンシャルは、この一様な背景電荷が存在するもとでのポテンシャルと正確に一致する。

第一原理計算（DFT）などで孤立分子やスラブモデルを周期境界セルに閉じ込めて荷電状態（チャージドシステム）の計算を行う際、この背景電荷の導入によってエネルギーが有限値に収まる。ただし、背景電荷自体と点電荷との相互作用、および背景電荷同士の自己相互作用が含まれるため、セルのサイズ $L$ に依存する人工的な補正項（Makov-Payne 補正など）が必要となる。

---

5.2 真空境界における表面・双極子補正の導出

Ewald法の実空間和と逆格子空間和への分割は、マクロな境界条件を一意に定めることで初めて完全なものとなる。1.3節で述べたように、無限格子和は条件収束であるため、系の外側の境界条件（周囲を取り囲む媒質）の性質によってエネルギーが変化する。

一般に、Ewald法の標準的な定式化（脱分極場を無視した形式）は、系の周囲が無限遠において導体（誘電率 $\epsilon_r = \infty$、tin-foil 境界条件）で囲まれていることを前提としている。しかし、現実の物理系や、スラブモデルのように一方向に真空領域を持つ異方的な周期境界系（真空境界、$\epsilon_r = 1$）をシミュレートする場合には、系が持つ正味の双極子モーメントによって誘起される脱分極場（depolarizing field）および表面電荷の寄与を適切に補正する必要がある。

周囲の誘電率を $\epsilon_r$ とし、体積 $V$ の基本単位胞内の全双極子モーメントを以下のように定義する：

\[\mathbf{M} = \sum_{i=1}^N q_i \mathbf{r}_i\]

このとき、周囲の媒質が任意の誘電率 $\epsilon_r$ を持つ場合の双極子エネルギー補正項 $E_{\text{dipole}}$ は、以下のように導出される：

\[E_{\text{dipole}} = \frac{2\pi}{(2\epsilon_r + 1) V} \left\lvert \sum_{i=1}^N q_i \mathbf{r}_i \right\rvert^2\]

この補正項が周囲の媒質の誘電率 $\epsilon_r$ の変化によってどのようにシフトするか、物理的極限を交えて解析する。

1. 錫箔境界（Tin-foil boundary: $\epsilon_r \to \infty$）: 周囲が完全導体（無限の誘電率）で満たされている場合、境界における静電遮蔽効果により、表面に誘起される分極電荷が完全に相殺される。このとき、分極電場はゼロとなるため：

   \[E_{\text{dipole}} = \lim_{\epsilon_r \to \infty} \frac{2\pi}{(2\epsilon_r + 1) V} \vert \mathbf{M} \vert^2 = 0\]

   これは、通常の Ewald 法で追加の補正を必要としない標準的な設定に対応する。

2. 真空境界（Vacuum boundary: $\epsilon_r = 1$）: 周囲が真空（あるいは気相）である場合、単位胞内の双極子モーメント $\mathbf{M}$ によって分極が生じ、境界の表面電荷によって内部に一様な脱分極場（depolarizing field）が誘起される。この場合、補正項は以下のようになる：

   \[E_{\text{dipole}} = \frac{2\pi}{3 V} \left\lvert \sum_{i=1}^N q_i \mathbf{r}_i \right\rvert^2\]

   このエネルギーは、双極子モーメントが作る巨視的な電場が真空空間に蓄える静電エネルギーに対応しており、tin-foil 境界条件を前提とした通常の Ewald エネルギーにこの項を加算することで、真空境界条件における正しい静電エネルギーが得られる。

3. 一般的な誘電媒質（General dielectric medium: $1 < \epsilon_r < \infty$）: 周囲が水（$\epsilon_r \approx 80$）や極性溶媒などの有限の誘電率を持つ媒質である場合、遮蔽効果は部分的になり、補正エネルギーは 0 から $\frac{2\pi}{3V}\vert\mathbf{M}\vert^2$ の間の値をとる。

特に、2次元の周期性（スラブモデル）を持つシステムで $z$ 軸方向にのみ周期性を切断し真空領域を設ける場合、双極子モーメントの $z$ 成分 $M_z$ によって生じる電場が問題となる。この場合、1次元的な脱分極場の相殺（Yeh-Berkowitz 補正）として、以下の形の補正項が導入される：

\[E_{\text{dipole}}^{1\text{D}} = \frac{2\pi}{V} M_z^2\]

これにより、3次元周期境界条件用の Ewald コードを用いながら、2次元スラブ境界を極めて厳密に再現することが可能となる。

---

結論

本稿では、周期境界条件下における長距離静電相互作用の評価における条件収束の問題から出発し、Ewald 法の厳密な数学的定式化について解説した。

Ewald 法の本質は、クーロンポテンシャル $1/r$ を誤差関数を用いて実空間の短距離成分と逆格子空間の長距離成分に厳密に分割し、それぞれ絶対収束する急速な和として処理する点にある。さらに、逆格子空間の計算で生じる人工的な自己相互作用の補正（$E_{\text{self}}$）、$\mathbf{G}=\mathbf{0}$ の除外による一様背景電荷（Jellium）の導入、そしてマクロな境界の幾何学的形状から生じる双極子補正（$E_{\text{dipole}}$）の数理を示すことで、この手法が備える数学的一貫性と物理的妥当性を明らかにした。

現代の分子シミュレーションや電子状態計算においては、Ewald 法の基本思想を拡張した粒子メッシュ Ewald（PME）法や高速多重極展開法（FMM）が標準的に用いられており、これらは巨大な生体高分子や複雑な凝縮系物理の解析において不可欠なインフラとなっている。本稿で扱った定式化の理解は、これらの高度なアルゴリズムの挙動や、境界条件の選択が物理量に与える影響を深く洞察するための確固たる基礎となるだろう。