前言

一、本文重科普轻推导，主要是介绍Ewald Summation 与 Partical Mesh Ewald为什么会 出现，它们做了什么以及为什么这么做。

二、book内容是对以下两篇文章的简述，欢迎大家阅读。

1. Z. Hu, J. Chem. Theory Comput. 2014, 10, 5254–5264.
2. Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577−8593.

三、有疑问可联系zhaoyihao@protonmail.com

若读者只对分子模拟有基本了解，那本文可以让你对Ewald Summation方法本身以及发展脉络有一个清晰的理解，若读者推导过Ewald Summmation方法并写过相应的code，那么恭喜你将是本文最大的受益者，深刻理解Ewald Summation背后的物理意义是比数学推导更加困难的事情。

Background

静电相互作用势为： $$u^{\mathrm{ele}}(\mathbf{r})=\frac{1}{r}\text{\hspace{1em}(1-1)}$$ 对于N个带电粒子的体系而言,体系的静电相互作用能为： $$E^{\mathrm{ele}}=\frac{1}{2}\sum _i^N\sum _{j\ne i}^N{q}_i{q}_j{u}^{\mathrm{ele}}({\mathbf{r}}_{ij})\text{\hspace{1em}(1-2)}$$ 受算力的限制我们无法模拟一个宏观体系，因此我们通常会施加周期性边界条件，使得能从微观体系中获得和宏观体系一般的性质。 此时对于周期性体系而言，其静电相互作用能为： E^{\mathrm{ele}} = \dfrac{1}{2}\sum_{\mathbf{n}} \sum_{i}^{N} \sum_{j}^{N} !' q_iq_j u^{\mathrm{ele}}(\mathbf{r}_{ij}+\mathbf{n}\mathbf{L}) \tag{1-3} 其中 $\mathbf{nL}=(n_x{L}_x,n_y{L}_y,n_z{L}_z)$ , $L_x,L_y,L_z$ 为盒子边长，$n_x,n_y,n_z\in \mathbb{Z}$ ，“$\prime$”代表当 $\mathbf{n}=0$ 时 $i\ne j$。 显然我们不能真的对所有的格矢 $\mathbf{n}$ 进行求和，需要像其它作用势那样加截断距离。然而静电势是$r^{-1}$衰减,在一维空间及以上时其从0积分至 $\infty$ 时是发散的, 不能直接加截断。 因此出现Ewald Summation方法。

$$\begin{array}{c}{\int }_0^{\infty }\frac{1}{r^a}drconvergesifa>1;itdivergesifa\le 1\end{array}\text{\hspace{1em}(1-4)}$$

Ewald Summation

1921年Paul Peter Ewald提出, 当静电势 $r^{-1}$ 分解为两部分: $$\frac{1}{r}\equiv \frac{\mathrm{erfc}(\alpha r)}{r}+\frac{\mathrm{erf}(\alpha r)}{r}\text{\hspace{1em}(2-1)}$$

$\alpha$ 是调控两部分占比的系数。由图1我们可以看到两部分中 $\frac{\mathrm{erfc}(\alpha r)}{r}$ 部分衰减迅速, $\frac{\mathrm{erf}(\alpha r)}{r}$ 部分衰减缓慢。分解后对 于衰减迅速的 $\frac{\mathrm{erfc}(\alpha r)}{r}$ 部分就可以像其它作用势一般加截断在实空间进行求和，对于 $\frac{\mathrm{erf}(\alpha r)}{r}$ 部分做特殊处理。

当作用势被分解后，相应地静电相互作用能也被分解为多部分：

$$E^{\mathrm{ele}}=E^R+(E^K+E^{\mathrm{IB}})-E^{\mathrm{self}}\text{\hspace{1em}(2-2)}$$

E^{R} = \dfrac{1}{2}\sum_{\mathbf{n}} \sum_{i}^{N} \sum_{j}^{N} !'q_iq_j \dfrac{\mathrm{erfc}(\alpha|\mathbf{r}{ij}+\mathbf{n}\mathbf{L}|)}{|\mathbf{r}{ij}+\mathbf{a}|} \tag{2-3} $$(E^K+E^{\mathrm{IB}})=\frac{1}{2}\sum _{\mathbf{n}}\sum _i^N\sum _j^N{q}_i{q}_j\frac{\mathrm{erf}(\alpha |{\mathbf{r}}_{ij}+\mathbf{nL}|)}{|{\mathbf{r}}_{ij}+\mathbf{nL}|}\text{\hspace{1em}(2-4)}$$

$$E^{\mathrm{self}}=\frac{1}{2}\sum _i^N{q}_i^2\underset{r\to 0}{\lim }\frac{\mathrm{erf}(\alpha r)}{r}=\frac{\alpha }{\sqrt{\pi }}\sum _i^N{q}_i^2\text{\hspace{1em}(2-5)}$$ $E^R$ 与 $(E^K+E^{\mathrm{IB}})$ 两项作用能分别对应分解出的两部分作用势，需注意在 $(E^K+E^{\mathrm{IB}})$ 中没有去掉 $\mathbf{n}=0$ 时 $i=j$ 项，即自相互作用项，因此需要额外减去 $E^{\mathrm{self}}$ 。

对于衰减缓慢的 $\frac{\mathrm{erf}(\alpha r)}{r}$ 部分所作的特殊处理是，不再将其在实空间对格矢 $\mathbf{nL}$ 进行求和, 而是变换到 Fourier空间(也称倒易空间)对波矢进行求和。

这一步是通过Poisson's summation formula完成的: $$\sum _{\mathbf{n}}f(\mathbf{r}+\mathbf{nL})=\frac{1}{V}\sum _{\mathbf{k}}\mathcal{F}[f(\mathbf{r})]e^{i\mathbf{k}\cdot \mathbf{r}}\text{\hspace{1em}(2-6)}$$

$$\begin{array}{rl}(E^K+E^{\mathrm{IB}})=& \frac{1}{2}\sum _{\mathbf{n}}\sum _i^N\sum _j^N{q}_i{q}_j\frac{erf(\alpha |{\mathbf{r}}_{ij}+\mathbf{nL}|)}{|{\mathbf{r}}_{ij}+\mathbf{nL}|}\\ =& \frac{1}{2}\sum _i^N\sum _j^N{q}_i{q}_j\frac{1}{V}\sum _{\mathbf{k}}\frac{4\pi }{k^2}\exp (-\frac{k^2}{4{\alpha }^2})\exp (i\mathbf{k}\cdot {\mathbf{r}}_{ij})\end{array}\text{\hspace{1em}(2-7)}$$

其中 $\mathbf{k}=(\frac{2\pi m_1}{L_x},\frac{2\pi m_2}{L_y},\frac{2\pi m_3}{L_z})m_1,m_2,m_3\in \mathbb{Z}$ 。由于 $k$ 处在分母位置上，因此将 $\mathbf{k}=0$ 项单列为 $E^{\mathrm{IB}}$ 称为无穷边界项，该项是条件收敛的且与体系的介电响应等性质息息相关(1921年Ewald提出Ewald Summation方法时并没有注意到这一项，直到1981年人们才认识到这一项的重要性，无穷边界项是Ewald Summation插图的最后一角，人们至今为止依旧没能完全理解这一项)。

$$\begin{array}{rl}{E}^K=& \frac{1}{2V}\sum _i^N\sum _j^N{q}_i{q}_j\sum _{\mathbf{k}\ne 0}\frac{4\pi }{k^2}\exp (-\frac{k^2}{4{\alpha }^2})\exp (i\mathbf{k}\cdot {\mathbf{r}}_{ij})\\ E^{\mathrm{IB}}=& \underset{\mathbf{k}\to 0}{\lim }\frac{1}{2V}\sum _i\sum _j{q}_i{q}_j\exp (-\frac{k^2}{4{\alpha }^2})\exp (i\mathbf{k}\cdot {\mathbf{r}}_{ij})\end{array}\text{\hspace{1em}(2-8)}$$

此时，对于 $E^K$ 项也是通过在倒空间加截断求和的(如图2中的虚线)。

将 $(E^K+E^{\mathrm{IB}})$ 项变换到倒空间对倒格矢进行求和的必要性是: 尽管我们依旧无法对无限的格矢求和, 但是对有限的倒格矢求和还是使得我们能够抓住两个相距无穷远的带电粒子之间的correlation。 而这是直接对静电势 $r^{-1}$ 在实空间加截断无法做到的，因为只有截断距离是有限的，那么有限+1尺度上的correlation还是无法抓住的。

那么自然而然地就存在一个问题，为什么对$(E^K+E^{\mathrm{IB}})$ 在倒易空间中对波矢加截断不会切断实空间中两个相距无穷远的电荷之间的关联？ 对比笔者认为最易于理解的角度是点阵的角度出发：由于我们模拟的体系是具有周期性边界条件可以将整个体系视作晶格，当把周期性体系抽象为点阵后，eq（1-3）中的对无穷的格矢求和其实就是对无限的点阵点求和，同样也等价于对无限的晶面族求和，而倒易点阵的点阵点与正点阵中晶面族一一对应（见附录）,对倒易点阵点加截断即为只对有限的晶面族求和， 而任一晶面族中都不仅包括包含离正点阵原点近的正点阵点也包括相距无穷远的正点阵点。 这就是为什么Ewald Summation能够抓住相距无穷远的电荷之间的关联。

Partical Mesh Ewald

从eq(2-3)与eq(2-8)中我们不难发现，不论是实空间能量 $E^R$ 项还是倒空间能量 $E^K$ 项，它们的复杂度均为 $\mathcal{O}(N^2)$ 。 对于实空间中的 $E^R$ 项我们可以使用neighbour list方法来降低复杂度。而对于倒空间的 $E^K$ 项则需要通过Partical Mesh Ewald方法来降低复杂度 至 $\mathcal{O}(NlogN)$

Partical Mesh Ewald的核心是Fast Fourier Transform (FFT), 而为了能够FFT需要先离散插值。 为保持与参考文献的一致性，将eq(2-8)中做 $\mathbf{k}=2\pi \mathbf{m}$ 替换： $$\begin{array}{rl}{E}^K=& \frac{1}{2\pi V}\sum _{\mathbf{m}\ne 0}\frac{1}{m^2}\exp (-\frac{{\pi }^2{m}^2}{{\alpha }^2})\sum _i^N\sum _j^N{q}_i{q}_j\exp (i2\pi \mathbf{m}\cdot {\mathbf{r}}_{ij})\\ =& \frac{1}{2\pi V}\sum _{\mathbf{m}\ne 0}\frac{1}{m^2}\exp (-\frac{{\pi }^2{m}^2}{{\alpha }^2})S(\mathbf{m})S(-\mathbf{m})\end{array}\text{\hspace{1em}(3-1)}$$

其中：

$$S(\mathbf{m})=\sum _i^N{q}_i{e}^{i2\pi \mathbf{m}\cdot {\mathbf{r}}_i}=\sum _i^N{q}_i\exp (i2\pi \frac{m_1}{K_1}{u}_1)\exp (i2\pi \frac{m_2}{K_2}{u}_2)\exp (i2\pi \frac{m_3}{K_3}{u}_3)\text{\hspace{1em}(3-2)}$$

u是带电粒子的分数坐标, 例 $u_1=K_1{x}_i/L_x$ 。

首先将 $S(\mathbf{m})$ 离散到$K_1\times K_2\times K_3$的网格上： $$\begin{array}{rl}S(\mathbf{m})\approx & \sum _i^N{q}_i{b}_1(m_1)\sum _{k_1=-\infty }^{\infty }{M}_n(u_{1i}-k_1-n_1{K}_1)\cdot \exp (i2\pi \frac{m_1}{K_1}{k}_1)\\ & \cdot b_2(m_2)\sum _{k_2=-\infty }^{\infty }{M}_n(u_{2i}-k_2-n_2{K}_2)\cdot \exp (i2\pi \frac{m_2}{K_2}{k}_2)\\ & \cdot b_3(m_3)\sum _{k_3=-\infty }^{\infty }{M}_n(u_{3i}-k_3-n_3{K}_3)\cdot \exp (i2\pi \frac{m_3}{K_3}{k}_3)\\ =& b_1(m_1)b_1(m_1)b_1(m_1)\mathcal{F}[Q](m_1,m_2,m_3)\end{array}\text{\hspace{1em}(3-3)}$$

其中$M_n$为b样条系数,可由 Cox-de Boor formula递归得到: $${\mathrm{M}}_{\mathrm{n}}(\mathrm{x})=\frac{\mathrm{x}}{\mathrm{n}-1}{\mathrm{M}}_{\mathrm{n}-1}(\mathrm{x})+\frac{\mathrm{n}-\mathrm{x}}{\mathrm{n}-1}{\mathrm{M}}_{\mathrm{n}-1}(\mathrm{u}-1)\text{\hspace{1em}(3-4)}$$

$${\mathrm{M}}_2(\mathrm{x})=\{\begin{array}{ll}1-|x-1|& \text{for }0\ge x\le 2\\ 0& \text{for }x<0\text{ or }x>2\end{array}\text{\hspace{1em}(3-5)}$$

需要说明的是样条插值的系数由线性方程组的解得到,但会出现不满秩的情况(出现基函数为0)。 对于复指数的样条插值有特殊的解决方法称为欧拉指数样条，在阶数$n$为偶数时：

b_{\alpha}(m_{\alpha}) = \exp(i2\pi(n-1)\frac{m_{\alpha}}{K_{\alpha}})\big{[} \sum_{k=0}^{n-2}M_n(k+1)\exp(i2\pi \frac{m_{\alpha}}{K_1}k)\big{]}^{-1} \tag{3-6}

eq(3-3)中被Fouier变换的函数 $Q(k_1,k_2,k_3)$ 的表达式如下, 其物理意义是将电荷离散到网格的格点上(如图3所示), $n_1,n_2,n_3\in \mathbb{Z}$ 的存在是保持电荷离散时格点在同一个周期内。

$$\begin{array}{rl}Q(k_1,k_2,k_3)=& \sum _i^N{q}_i\sum _{n_1,n_2,n_3}{M}_n(u_{1i}-k_1-n_1{K}_1)\\ & \times M_n(u_{2i}-k_2-n_2{K}_2)\times M_n(u_{3i}-k_3-n_3{K}_3)\end{array}\text{\hspace{1em}(3-7)}$$

此时将离散后的 $S(\mathbf{m})$ 与 $S(-\mathbf{m})$ 带入到 $E^K$ 中： $$\begin{array}{rl}{E}^K=& \frac{1}{2\pi V}\sum _{\mathbf{m}\ne 0}\frac{1}{m^2}{e}^{-\frac{{\pi }^2{m}^2}{{\alpha }^2}}B(m_1,m_2,m_3)\\ & \cdot \mathcal{F}[Q](m_1,m_2,m_3)\mathcal{F}[Q](-m_1,-m_2,-m_3)\end{array}\text{\hspace{1em}(3-8)}$$

其中 $B(m_1,m_2,m_3)=|b_1(m_1){|}^2|b_2(m_2){|}^2|b_3(m_3){|}^2$ 。此时对于 $\mathcal{F}[Q]$ 使用Cooley-Tukey FFT算法已经降低复杂度至 $\mathcal{O}(N\log N)$ 了，但通常会更进一步简化，即通过 ${\sum }_{\mathbf{m}}\mathcal{F}[A]\cdot B={\sum }_{\mathbf{m}}A\cdot \mathcal{F}[B]$ 以及卷积的 Fourier 变换等于 Fourier 的乘积得到最终表达式: $$\begin{array}{rl}{E}^K=& \frac{1}{2\pi V}\sum _{\mathbf{m}\ne 0}Q(m_1,m_2,m_3)(\theta \star Q)((m_1,m_2,m_3)\end{array}\text{\hspace{1em}(3-9)}$$ 其中： $$\theta \star Q=\mathcal{F}[\frac{1}{m^2}{e}^{-\frac{{\pi }^2{m}^2}{{\alpha }^2}}\cdot B(m_1,m_2,m_3)\cdot \mathcal{F}[Q](-m_1,-m_2,-m_3)]\text{\hspace{1em}(3-10)}$$

PME code 实现

示例代码力求可读性，不具备工作效率。

求解倒易空间的能量项时，所需的的物理量如下：

1. N ；粒子的总数目;
2. $\{{\mathbf{r}}_1,{\mathbf{r}}_2,...,{\mathbf{r}}_N\}$ : 粒子坐标的集合；
3. $\{{\mathbf{q}}_1,{\mathbf{q}}_2,...,{\mathbf{q}}_N\}$: 粒子电荷的集合；
4. $K_1,K_2,K_3$ : 三个维度上离散格点的数目(三者必须均为2的整数次幂,Cooley-Tukey FFT算法结构决定);
5. $\mathrm{n}$ : B样条阶数,(必须为偶数，且小于等于 $\min (K_1,K_2,K_3)$，由离散电荷时不超过一个周期决定)

• 第一步：计算欧拉指数样条系数 $B(m_1,m_2,m_3)$

1. 通过Cox-de Boor formula eq(3-4)和初始条件eq(3-5)递归计算B样条系数 ${\mathrm{M}}_n(i)i=1,\dots ,n-1$(递归结构如图4所示);
2. 由eq(3-6)计算 ${\mathrm{b}}_1,{\mathrm{b}}_2,{\mathrm{b}}_3$;
3. $\mathrm{B}({\mathrm{m}}_1,{\mathrm{m}}_2,{\mathrm{m}}_3)=|{\mathrm{b}}_1({\mathrm{m}}_1){|}^2|{\mathrm{b}}_2({\mathrm{m}}_2){|}^2|{\mathrm{b}}_3({\mathrm{m}}_3){|}^2$。

• 第二步计算实空间中电荷网格 $Q(k_1,k_2,k_3)$

1. 通过Cox-de Boor formula eq(3-4)和初始条件eq(3-5)递归计算B样条系数 ${\mathrm{M}}_n(x+i)x=u_{\alpha }(j)-\mathrm{int}(u_{\alpha }(j)),i=1,\dots ,n-1$ (递归结构如图5所示);
2. 由eq(3-7)离散电荷到网格上得到 $Q(k_1,k_2,k_3)$ 。

• 第三步对电荷网格做FFT变换$\mathcal{F}[Q](-m_1,-m_2,-m_3)$

$$\begin{array}{rl}\mathcal{F}[Q](-m_1,-m_2,-m_3)=& \sum _{k_1=0}^{K_1-1}\sum _{k_2=0}^{K_2-1}\sum _{k_3=0}^{K_3-1}Q(k_1,k_2,k_3)\\ & \cdot \exp [i2\pi (\frac{-m_1{k}_1}{K_1}+\frac{-m_2{k}_2}{K_2}+\frac{-m_3{k}_3}{K_3})]\end{array}\text{\hspace{1em}(4-1)}$$

• 第四步计算 coefficient $\frac{\exp (-{\pi }^2{\mathbf{m}}^2/{\alpha }^2)}{{\mathbf{m}}^2}$, $B$，$\mathcal{F}[Q](-m_1,-m_2,-m_3)$ 三者乘积CBF

1. 注意倒易空间具有截断.

• 第五步计算第四步三者乘积CBF的Fouier变换 $\theta \star Q(k_1,k_2,k_3)$

$$\begin{array}{rl}\theta \star Q(k_1,k_2,k_3)=& \sum _{m_1=0}^{K_1-1}\sum _{m_2=0}^{K_2-1}\sum _{m_3=0}^{K_3-1}\mathrm{CBF}(m_1,m_2,m_3)\\ & \cdot \exp [i2\pi (\frac{m_1{k}_1}{K_1}+\frac{m_2{k}_2}{K_2}+\frac{m_3{k}_3}{K_3})]\end{array}\text{\hspace{1em}(4-2)}$$

• 第六步计算倒易空间能量项 $E^K$

1. 计算$Q(k_1,k_2,k_3)$与$(\theta \star Q)(k_1,k_2,k_3)$的乘积;
2. 由eq(3-9)计算乘积得到的数组之和，乘上系数与被reduce的$\frac{1}{4\pi {\epsilon }_0}$得到 $E^K$

附录

倒易空间

• 点阵

当在分子模拟中使用周期性边界条件时可视为晶体, 而晶体则可以抽象为点阵。额外需要说明的是由于分子模拟中Box里的粒子通常并不是直观有序的, 因此粒子并不需要一定处在晶格点上, 但粒子的排列必须满足晶格对称性(对于直观有序的如NaCl晶体,粒子可以不在晶格点上)。

• 密勒指数

密勒指数就是取晶面族在坐标轴上的截距，截距的单位就是基向量大小，取倒数。

如图6所示，该族晶面在${\mathbf{a}}_1$轴上的截距是$\frac{1}{2}{a}_1$，在${\mathbf{a}}_2$轴上的截距是$\frac{1}{3}{a}_2$，同时与${\mathbf{a}}_3$轴平行截距为$\infty$， 因此是(230)晶面族。正点阵中晶面族的数目是无穷的。

• 倒基矢

倒格基矢 $({\mathbf{b}}_1,{\mathbf{b}}_2,{\mathbf{b}}_3)$ 的方向分别为(100)、(010)、(001)三个晶面族的法向量方向，大小是三个晶面族的间距的倒数乘以 $2\pi$ 。下面我们以 ${\mathbf{b}}_1$ 为例来推导倒基矢与正基矢的关系，其它方向上的倒基矢推导同理。

$$|{\mathbf{b}}_1|=\frac{2\pi }{d}d=\frac{{\mathbf{a}}_1\cdot ({\mathbf{a}}_2\times {\mathbf{a}}_3)}{|{\mathbf{a}}_2\times {\mathbf{a}}_3|}\text{\hspace{1em}(A1)}$$

$d$ 就是(100)晶面族的间距，它等与${\mathbf{a}}_1$在该族晶面法向量上的投影，而${\mathbf{a}}_2$叉乘${\mathbf{a}}_3$就是该族晶面的法向量。因此${\mathbf{b}}_1$就等于$|{\mathbf{b}}_1|$乘上单位法向量： $${\mathbf{b}}_1=\frac{2\pi }{\frac{{\mathbf{a}}_1\cdot ({\mathbf{a}}_2\times {\mathbf{a}}_3)}{|{\mathbf{a}}_2\times {\mathbf{a}}_3|}}\frac{{\mathbf{a}}_2\times {\mathbf{a}}_3}{|{\mathbf{a}}_2\times {\mathbf{a}}_3|}=2\pi \frac{{\mathbf{a}}_2\times {\mathbf{a}}_3}{{\mathbf{a}}_1\cdot ({\mathbf{a}}_2\times {\mathbf{a}}_3)}\text{\hspace{1em}(A2)}$$ 由上述表达式我们可得： $${\mathbf{a}}_1\cdot {\mathbf{b}}_1=2\pi {\mathbf{a}}_2\cdot {\mathbf{b}}_1={\mathbf{a}}_3\cdot {\mathbf{b}}_1=0\text{\hspace{1em}(A3)}$$ 因此我们可以看到一些教材上是这样写的: $${\mathbf{a}}_i\cdot {\mathbf{b}}_j=2\pi {\delta }_{ij}\text{\hspace{1em}(A4)}$$

• 倒易点阵

倒易点阵的一个格矢是和正点阵中的一族晶面相对应的, 又因为一个格点对应一个格矢, 因此可以说倒易点阵中一个点阵点对应正点阵中的一族晶面。

Fast Fourier Transform

按笔者的理解, FFT的核心就是利用复指数的特性eq(B1)再加上分而治之的思想(图7)。

$$\begin{array}{rl}{e}^{i\pi \frac{2k}{2n}}=& e^{i\pi \frac{k}{n}}\\ e^{i\pi \frac{k+n/2}{n}}=& -e^{i\pi \frac{k}{n}}\end{array}\text{\hspace{1em}(B1a)(B1b)}$$

根据分而治之的思想中可以获得两点知识, 一是说明为什么要离散为2的整数次幂, 二是FFT的复杂度$\mathcal{O}(N\log N)$中的 $\log N$ 是以2为底, 这同样也是分而治之的层数。 受篇幅所限不再对FFT做过多介绍。