本NB初衷：

在 AISI 电池团队的一个项目中，需要求周期性体系的静电相互作用能，最初是采用原始的 Ewald 进行实现，后期为了提高效率，决定采用 PME 替代 Ewald，因此有了本NB。

前言

在分子动力学（molecular dynamics, MD）模拟中，在为模拟系统选取边界条件时经常会采用周期性边界条件（periodic boundary conditions, PBC）以降低边界效应对系统中粒子的影响。考虑真空中单元系统内的 $N$ 个点电荷（point charge）$q_1,q_2,...,q_N$，坐标分别为 ${\mathbf{r}}_1,{\mathbf{r}}_2,...,{\mathbf{r}}_N$，且满足条件 ${\sum }_{i=1}^N{q}_i=0$。使用周期性边界条件的无穷大系综的静电相互作用能（electrostatic energy）的总和为 $$E=\frac{1}{4\pi {\epsilon }_0}\frac{1}{2}\sum _{\mathbf{n}}^{\ast }\sum _{i,j=1}^N\frac{q_i{q}_j}{\left|{\mathbf{r}}_{ij}+\mathbf{n}\right|}\text{\hspace{1em}(1)}$$ 其中 ${\epsilon }_0=8.854\times 10^{-12}{\mathrm{C}}^2{\mathrm{N}}^{-1}{\mathrm{m}}^{-2}$ 为真空中介电常数（vacuum permittivity），${\mathbf{r}}_{ij}={\mathbf{r}}_j-{\mathbf{r}}_i$，$\mathbf{n}=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}}_{\alpha }$（$\alpha =1,2,3$）为单元系统的边界矢量。引入 $\ast$ 记号表示 $i\ne j$，当且仅当 $\mathbf{n}=0$，即排除单元系统中点电荷与自身的静电相互作用。方程 $1$ 指出，在周期性的无穷大系综中，直接计算带电粒子的静电相互作用能是困难的，这是因为静电相互作用能随 $r^{-1}$ 衰减，并且能量的无穷项求和是条件收敛的。

德国物理学家 Ewald 为了计算晶体内静电相互作用能，于1921年提出 Ewald 求和算法1。包含无穷多重复单元的周期性系综和包含无穷多相同晶胞的晶体在计算性质上是高度相似的，因此在分子动力学模拟中计算静电相互作用也可以采用相同的算法。Ewald 求和的时间复杂度为 $\mathcal{O}(N^2)$，而为了处理大体量系统中的静电相互作用，Darden、York 和 Pedersen 于1993年基于经典 Ewald 求和提出了 particle mesh Ewald（PME）求和，将 Ewald 求和中最耗时的一项的时间复杂度缩减为 $\mathcal{O}(N{\log }_{2}N)$2。1995年，Essmann、Perera 和 Berkowitz 等人通过改良 PME 求和的算法提高了计算精度，提出 smooth particle mesh Ewald（SPME）求和算法3。

本文将简要介绍 Ewald 求和和 PME 求和的基本原理，并运用 C 语言实现 PME 求和对3D周期性系综的静电相互作用能的计算。

Ewald 求和

基本思想

Ewald 求和的基本思想是引入一个快速衰减函数 $f(r)$，将 $r^{-1}$ 阶的长程库仑相互作用重写为两项之和： $$\frac{1}{r}=\frac{f(r)}{r}+\frac{1-f(r)}{r}\text{\hspace{1em}(2)}$$ 其中，第一项 $\frac{f(r)}{r}$ 是快速衰减的，且其求和是收敛的。第二项 $\frac{1-f(r)}{r}$ 是慢变函数，可以通过 Fourier 变换将其变换到倒易空间（reciprocal space）后对波矢（wave vector）$\mathbf{k}$ 求和得到。Ewald 求和选取了余误差函数（complementary error function）作为衰减函数，其表达式为： $$\mathrm{erfc}(x)=1-\mathrm{erf}(x)=1-\frac{1}{\sqrt{\pi }}{\int }_{-x}^x{e}^{-t^2}\mathrm{d}t=1-\frac{2}{\sqrt{\pi }}{\int }_0^x{e}^{-t^2}\mathrm{d}t\text{\hspace{1em}(3)}$$

图 1. 通过选择余误差函数作为衰减函数，将长程静电相互作用（黑色实线）分离为短程项（蓝色实线）和慢衰减项（红色虚线）。

物理图像

下面将探讨方程 $2$ 中对点电荷间长程静电相互作用进行的分离的物理意义。点电荷在空间中的电荷密度分布可以通过 Dirac delta 函数的集合来描述。特别地，点电荷 $q_i$ 在位置 $\mathbf{r}$ 处的电荷密度为 $${\rho }_i^{\mathrm{P}}(\mathbf{r})=q_i\delta (\mathbf{r}-{\mathbf{r}}_i)\text{\hspace{1em}(4)}$$ 其中 Dirac delta 函数满足 $\delta (\mathbf{r})=0(\mathbf{r}\ne 0)$，且 $\int \delta (\mathbf{r})\mathrm{d}\mathbf{r}=1$。 假想以点电荷的位置为中心在空间中引入一个总电荷量相同、电性相反、球对称（spherical symmetry）且平滑变化的电荷分布。这种云状的电荷分布有效地屏蔽了点电荷的长程静电相互作用（可以设想在无穷远处，电性相反的云状电荷分布近似等同于点电荷，完全抵消点电荷的静电作用），仅保留短程静电相互作用。同时为了补偿引入的电荷分布，还需引入与云状电荷分布形状相同，且与点电荷电性相同另一种电荷分布。在 Ewald 求和中选取 Gauss 分布描述云状电荷分布： $${\rho }_i^{\mathrm{G}}(\mathbf{r})=q_i{\mathrm{G}}_{\sigma }(\mathbf{r}-{\mathbf{r}}_i)\text{\hspace{1em}(5)}$$ $${\mathrm{G}}_{\sigma }(\mathbf{r})=\frac{1}{(2\pi {\sigma }^2{)}^{3/2}}\exp \left(-\frac{{\left|\mathbf{r}\right|}^2}{2{\sigma }^2}\right)\text{\hspace{1em}(6)}$$ 其中，$\sigma$ 是 Gauss 分布的标准差，指示其分布宽度。值得注意的是，delta 函数也可以视为 Gauss 分布对 $\sigma \to 0$ 的极限，即 $$\underset{\sigma \to 0}{\lim }{\mathrm{G}}_{\sigma }(\mathbf{r})=\delta (\mathbf{r})\text{\hspace{1em}(7)}$$ 规定与对应点电荷电性相同且服从 Gauss 分布的电荷称为 Gauss 电荷，与 Gauss 电荷形状相同且电性相反的电荷称为反 Gauss 电荷，点电荷与反 Gauss 电荷合并称为屏蔽电荷（screened charge）。

图2. 一组点电荷可以被认为是由 (a) 一组屏蔽电荷与 (b) 一组 Gauss 电荷组合而成。颜色用于指示电荷的符号。

引入 Poisson 方程： $${\nabla }^2\varphi (\mathbf{r})=-\frac{\rho (\mathbf{r})}{{\epsilon }_0}\text{\hspace{1em}(8)}$$ 其中 $\varphi (\mathbf{r})$ 为 $\mathbf{r}$ 处的静电势（electrostatic potential）。可以解得点电荷 $q_i$ 在 $\mathbf{r}$ 处产生的静电势为 $${\varphi }_i^{\mathrm{P}}(\mathbf{r})=\frac{1}{4\pi {\epsilon }_0}\frac{q_i}{\left|\mathbf{r}-{\mathbf{r}}_i\right|}\text{\hspace{1em}(9)}$$ Gauss 电荷在 $\mathbf{r}$ 处产生的静电势为 \phi_{i}^{\mathrm{G}}(\mathbf{r}) = \frac{1}{4\pi\varepsilon_{0}} \frac{q_i}{\left|\mathbf{r}-\mathbf{r}_i\right|} \operatorname{erf}\left(\frac{\left|\mathbf{r}-\mathbf{r}_i\right|}{\sqrt{2}\sigma}\right) \tag{10}\ 下面将分别在实空间（real space）中对屏蔽电荷、在倒易空间中对 Gauss 电荷的静电相互作用能进行求和，最后将两项相加即可得到系综的静电相互作用能总和。

实空间

点电荷 $q_i$ 与其对应的反 Gauss 电荷组成的屏蔽电荷在 $\mathbf{r}$ 处产生的静电势为 $$\begin{array}{rl}{\varphi }_i^{\mathrm{S}}(\mathbf{r})& ={\varphi }_i^{\mathrm{P}}(\mathbf{r})-{\varphi }_i^{\mathrm{G}}(\mathbf{r})\\ & =\frac{1}{4\pi {\epsilon }_0}\frac{q_i}{\left|\mathbf{r}-{\mathbf{r}}_i\right|}-\frac{1}{4\pi {\epsilon }_0}\frac{q_i}{\left|\mathbf{r}-{\mathbf{r}}_i\right|}\mathrm{erf}\left(\frac{\left|\mathbf{r}-{\mathbf{r}}_i\right|}{\sqrt{2}\sigma }\right)\\ & =\frac{1}{4\pi {\epsilon }_0}\frac{q_i}{\left|\mathbf{r}-{\mathbf{r}}_i\right|}\mathrm{erfc}\left(\frac{\left|\mathbf{r}-{\mathbf{r}}_i\right|}{\sqrt{2}\sigma }\right)\end{array}\text{\hspace{1em}(11)}$$ 点电荷 $q_j$ 与系综中除自身的屏蔽电荷外所有屏蔽电荷的静电相互作用能为 $$\begin{array}{rl}{E}_j^{\mathrm{real}}& =\frac{1}{2}\sum _{\mathbf{n}}^{\ast }\sum _{i=1}^N{q}_j{\varphi }_i^{\mathrm{S}}({\mathbf{r}}_j+\mathbf{n})\\ & =\frac{1}{2}\sum _{\mathbf{n}}^{\ast }\sum _{i=1}^N\frac{1}{4\pi {\epsilon }_0}\frac{q_i{q}_j}{\left|{\mathbf{r}}_j-{\mathbf{r}}_i+\mathbf{n}\right|}\mathrm{erfc}\left(\frac{\left|{\mathbf{r}}_j-{\mathbf{r}}_i+\mathbf{n}\right|}{\sqrt{2}\sigma }\right)\\ & =\frac{1}{4\pi {\epsilon }_0}\frac{1}{2}\sum _{\mathbf{n}}^{\ast }\sum _{i=1}^N\frac{q_i{q}_j}{\left|{\mathbf{r}}_{ij}+\mathbf{n}\right|}\mathrm{erfc}\left(\frac{\left|{\mathbf{r}}_{ij}+\mathbf{n}\right|}{\sqrt{2}\sigma }\right)\end{array}\text{\hspace{1em}(12)}$$ 单元系统在实空间的静电相互作用能为 $$E^{\mathrm{real}}=\sum _{j=1}^N{E}_j^{\mathrm{real}}=\frac{1}{4\pi {\epsilon }_0}\frac{1}{2}\sum _{\mathbf{n}}^{\ast }\sum _{i,j=1}^N\frac{q_i{q}_j}{\left|{\mathbf{r}}_{ij}+\mathbf{n}\right|}\mathrm{erfc}\left(\frac{\left|{\mathbf{r}}_{ij}+\mathbf{n}\right|}{\sqrt{2}\sigma }\right)\text{\hspace{1em}(13)}$$ 此外，还需单独计算点电荷与其对应的反 Gauss 电荷的静电相互作用能。Gauss 电荷在 ${\mathbf{r}}_i$ 处产生的静电势为 $${\varphi }_i^{\mathrm{G}}({\mathbf{r}}_i)=\frac{1}{4\pi {\epsilon }_0}\frac{q_i}{\sigma }\sqrt{\frac{2}{\pi }}\text{\hspace{1em}(14)}$$ 自相互作用能为 $$E^{\mathrm{self}}=\frac{1}{2}\sum _{i=1}^N{q}_i[-{\varphi }_i^{\mathrm{G}}({\mathbf{r}}_i)]=-\frac{1}{4\pi {\epsilon }_0}\sum _{i=1}^N\frac{q_i^2}{\sqrt{2\pi }\sigma }\text{\hspace{1em}(15)}$$

倒易空间

系综中 Gauss 电荷在空间中的密度分布为 $$\begin{array}{rl}{\rho }^{\mathrm{G}}(\mathbf{r})& =\sum _{\mathbf{n}}\sum _{i=1}^N{\rho }_i^{\mathrm{G}}(\mathbf{r}+\mathbf{n})\\ & =\sum _{\mathbf{n}}\sum _{i=1}^N{q}_i{\mathrm{G}}_{\sigma }(\mathbf{r}-{\mathbf{r}}_i+\mathbf{n})\\ & =\sum _{\mathbf{n}}\sum _{i=1}^N\frac{q_i}{(2\pi {\sigma }^2{)}^{3/2}}\exp \left(-\frac{{\left|\mathbf{r}-{\mathbf{r}}_i+\mathbf{n}\right|}^2}{2{\sigma }^2}\right)\end{array}\text{\hspace{1em}(16)}$$ 可知 Gauss 电荷在空间中周期性分布，对其进行 Fourier 变换： $$\begin{array}{rl}{\tilde{\rho }}^{\mathrm{G}}(\mathbf{k})& ={\int }_V{\rho }^{\mathrm{G}}(\mathbf{r})\exp (-2\pi i\mathbf{k}\cdot \mathbf{r})\mathrm{d}\mathbf{r}\\ & ={\int }_V\mathrm{d}\mathbf{r}\exp (-2\pi i\mathbf{k}\cdot \mathbf{r})\sum _{\mathbf{n}}\sum _{j=1}^N\frac{q_j}{(2\pi {\sigma }^2{)}^{3/2}}\exp \left(-\frac{{\left|\mathbf{r}-{\mathbf{r}}_j+\mathbf{n}\right|}^2}{2{\sigma }^2}\right)\\ & ={\int }_{{\mathbb{R}}^3}\mathrm{d}\mathbf{r}\exp (-2\pi i\mathbf{k}\cdot \mathbf{r})\sum _{j=1}^N\frac{q_j}{(2\pi {\sigma }^2{)}^{3/2}}\exp \left(-\frac{{\left|\mathbf{r}-{\mathbf{r}}_j\right|}^2}{2{\sigma }^2}\right)\\ & =\sum _{j=1}^N\frac{q_j}{(2\pi {\sigma }^2{)}^{3/2}}{\int }_{{\mathbb{R}}^3}\mathrm{d}\mathbf{r}\exp (-2\pi i\mathbf{k}\cdot \mathbf{r})\exp \left(-\frac{{\left|\mathbf{r}-{\mathbf{r}}_j\right|}^2}{2{\sigma }^2}\right)\\ & =\sum _{j=1}^N{q}_j\exp (-2\pi i\mathbf{k}\cdot {\mathbf{r}}_j)\exp (-2{\pi }^2{\sigma }^2{k}^2)\end{array}\text{\hspace{1em}(17)}$$ 其中倒易格矢量（reciprocal lattice vector）$\mathbf{k}=k_1{\mathbf{a}}_1^{\ast }+k_2{\mathbf{a}}_2^{\ast }+k_3{\mathbf{a}}_3^{\ast }$（$k_1,k_2,k_3\in \mathbb{Z}$），$k=\left|\mathbf{k}\right|$，倒易矢量（reciprocal vector）${\mathbf{a}}_{\alpha }^{\ast }$（$\alpha =1,2,3$）由 ${\mathbf{a}}_{\alpha }^{\ast }\cdot {\mathbf{a}}_{\beta }={\delta }_{\alpha \beta }$（Kronecker delta 函数）定义。 同理，Gauss 电荷产生的静电势 ${\varphi }^{\mathrm{G}}(\mathbf{r})$ 也在空间中周期性分布。对 ${\varphi }^{\mathrm{G}}(\mathbf{r})$ 做 Fourier 变换： $${\tilde{\varphi }}^{\mathrm{G}}(\mathbf{k})={\int }_V{\varphi }^{\mathrm{G}}(\mathbf{r})\exp (-2\pi i\mathbf{k}\cdot \mathbf{r})\mathrm{d}\mathbf{r}\text{\hspace{1em}(18)}$$ 两者的逆 Fourier 变换分别为 $${\rho }^{\mathrm{G}}(\mathbf{r})=\frac{1}{V}\sum _{\mathbf{k}}{\tilde{\rho }}^{\mathrm{G}}(\mathbf{k})\exp (2\pi i\mathbf{k}\cdot \mathbf{r})\text{\hspace{1em}(19)}$$ $${\varphi }^{\mathrm{G}}(\mathbf{r})=\frac{1}{V}\sum _{\mathbf{k}}{\tilde{\varphi }}^{\mathrm{G}}(\mathbf{k})\exp (2\pi i\mathbf{k}\cdot \mathbf{r})\text{\hspace{1em}(20)}$$ 将方程 $19$、$20$ 代入方程 $8$ 即 Poisson 方程得 $$\begin{array}{rl}{\nabla }^2{\varphi }^{\mathrm{G}}(\mathbf{r})& =-\frac{{\rho }^{\mathrm{G}}(\mathbf{r})}{{\epsilon }_0}\\ \frac{{\partial }^2}{\partial {\mathbf{r}}^2}\frac{1}{V}\sum _{\mathbf{k}}{\tilde{\varphi }}^{\mathrm{G}}(\mathbf{k})\exp (2\pi i\mathbf{k}\cdot \mathbf{r})& =-\frac{1}{{\epsilon }_0}\frac{1}{V}\sum _{\mathbf{k}}{\tilde{\rho }}^{\mathrm{G}}(\mathbf{k})\exp (2\pi i\mathbf{k}\cdot \mathbf{r})\\ -4{\pi }^2{k}^2{\tilde{\varphi }}^{\mathrm{G}}(\mathbf{k})& =-\frac{{\tilde{\rho }}^{\mathrm{G}}(\mathbf{k})}{{\epsilon }_0}\end{array}$$ $${\tilde{\varphi }}^{\mathrm{G}}(\mathbf{k})=\frac{{\tilde{\rho }}^{\mathrm{G}}(\mathbf{k})}{4{\pi }^2{\epsilon }_0{k}^2}\text{\hspace{1em}(21)}$$ 方程 $21$ 指出，Gauss 电荷与其静电势在倒易空间的关系非常简洁。和将方程 $21$ 代入方程 $20$，得 $$\begin{array}{rl}{\varphi }^{\mathrm{G}}(\mathbf{r})& =\frac{1}{V}\sum _{\mathbf{k}\ne 0}\frac{{\tilde{\rho }}^{\mathrm{G}}(\mathbf{k})}{4{\pi }^2{k}^2{\epsilon }_0}\exp (2\pi i\mathbf{k}\cdot \mathbf{r})\\ & =\frac{1}{4\pi {\epsilon }_0}\frac{1}{\pi V}\sum _{\mathbf{k}\ne 0}\sum _{j=1}^N\frac{q_j}{k^2}\exp (-2\pi i\mathbf{k}\cdot {\mathbf{r}}_j)\exp (2\pi i\mathbf{k}\cdot \mathbf{r})\exp (-2{\pi }^2{\sigma }^2{k}^2)\\ & =\frac{1}{4\pi {\epsilon }_0}\frac{1}{\pi V}\sum _{\mathbf{k}\ne 0}\sum _{j=1}^N\frac{q_j}{k^2}\exp (2\pi i\mathbf{k}\cdot (\mathbf{r}-{\mathbf{r}}_j))\exp (-2{\pi }^2{\sigma }^2{k}^2)\end{array}\text{\hspace{1em}(22)}$$ 注意到方程 $22$ 仅对 $\mathbf{k}\ne 0$ 成立。当单元系统为电中性即 ${\sum }_{i=1}^N{q}_i=0$ 时，$\mathbf{k}=0$ 项的贡献为零。系综中 Gauss 电荷对单元系统内点电荷的静电相互作用能的总和为 $$\begin{array}{rl}{E}^{\mathrm{recip}}& =\frac{1}{2}\sum _{i=1}^N{q}_i{\varphi }^{\mathrm{G}}(r_i)\\ & =\frac{1}{4\pi {\epsilon }_0}\frac{1}{2\pi V}\sum _{\mathbf{k}\ne 0}\sum _{i,j=1}^N\frac{q_i{q}_j}{k^2}\exp (2\pi i\mathbf{k}\cdot ({\mathbf{r}}_i-{\mathbf{r}}_j))\exp (-2{\pi }^2{\sigma }^2{k}^2)\\ & =\frac{1}{4\pi {\epsilon }_0}\frac{1}{2\pi V}\sum _{\mathbf{k}\ne 0}\frac{\exp (-2{\pi }^2{\sigma }^2{k}^2)}{k^2}\sum _{i,j=1}^N{q}_i\exp (2\pi i\mathbf{k}\cdot {\mathbf{r}}_i)q_j\exp (-2\pi i\mathbf{k}\cdot {\mathbf{r}}_j)\\ & =\frac{1}{4\pi {\epsilon }_0}\frac{1}{2\pi V}\sum _{\mathbf{k}\ne 0}\frac{\exp (-2{\pi }^2{\sigma }^2{k}^2)}{k^2}S(\mathbf{k})S(-\mathbf{k})\\ & =\frac{1}{4\pi {\epsilon }_0}\frac{1}{2\pi V}\sum _{\mathbf{k}\ne 0}\frac{\exp (-2{\pi }^2{\sigma }^2{k}^2)}{k^2}{\left|S(\mathbf{k})\right|}^2\end{array}\text{\hspace{1em}(23)}$$ 其中结构因子 $S(\mathbf{k})$ 的定义为 $$S(\mathbf{k})=\sum _{j=1}^N{q}_j\exp (2\pi i\mathbf{k}\cdot {\mathbf{r}}_j)\text{\hspace{1em}(24)}$$

本章小结

在 Ewald 求和中，系统的静电相互作用能的总和可以写为实空间、倒易空间和自相互作用的静电相互作用能三项之和： $$E^{\mathrm{elec}}=E^{\mathrm{real}}+E^{\mathrm{recip}}+E^{\mathrm{self}}\text{\hspace{1em}(25)}$$ 其中 $$E^{\mathrm{real}}=\frac{1}{4\pi {\epsilon }_0}\frac{1}{2}\sum _{\mathbf{n}}^{\ast }\sum _{i,j=1}^N\frac{q_i{q}_j}{\left|{\mathbf{r}}_{ij}+\mathbf{n}\right|}\mathrm{erfc}\left(\frac{\left|{\mathbf{r}}_{ij}+\mathbf{n}\right|}{\sqrt{2}\sigma }\right)\text{\hspace{1em}(26)}$$ $$E^{\mathrm{recip}}=\frac{1}{4\pi {\epsilon }_0}\frac{1}{2\pi V}\sum _{\mathbf{k}\ne 0}\frac{\exp (-2{\pi }^2{\sigma }^2{k}^2)}{k^2}{\left|S(\mathbf{k})\right|}^2\text{\hspace{1em}(27)}$$ $$E^{\mathrm{self}}=-\frac{1}{4\pi {\epsilon }_0}\sum _{i=1}^N\frac{q_i^2}{\sqrt{2\pi }\sigma }\text{\hspace{1em}(28)}$$ 静电相互作用能的实空间项 $E^{\mathrm{real}}$ 快速衰减，因此在实际计算中可以通过设置截断距离（cut-off distance）$r_{\mathrm{cut}-\mathrm{off}}$ 得到带电粒子的近邻列表（neighbor list），并在列表内直接求和得到。若点电荷 $q_i$ 的近邻列表为集合 ${\mathbb{L}}_i$，有 $$E^{\mathrm{real}}=\frac{1}{4\pi {\epsilon }_0}\frac{1}{2}\sum _{i=1}^N\sum _{j\in {\mathbb{L}}_i}\frac{q_i{q}_j}{\left|{\mathbf{r}}_{ij}\right|}\mathrm{erfc}\left(\frac{\left|{\mathbf{r}}_{ij}\right|}{\sqrt{2}\sigma }\right)\text{\hspace{1em}(29)}$$ Gauss 分布的标准差 $\sigma$ 的选择可以调节实空间和倒易空间所需算力的相对大小，从而影响 Ewald 求和的计算效率，但 Ewald 求和的最终结果与 $\sigma$ 的选择无关。Jackson 和 Catlow 通过预设每个实空间和倒易空间项的计算量相等，提出以下 $\sigma$ 的选取4： $$\sigma =\frac{1}{\sqrt{2\pi }}{\left(\frac{V^2}{N}\right)}^{\frac{1}{6}}\text{\hspace{1em}(30)}$$ 其中 $V$ 为单元系统的体积。此外，$\alpha \equiv 1/(2{\sigma }^2)$ 和 $\beta \equiv 1/(\sqrt{2}\sigma )$ 也是文献中常用的参数。

PME 求和

基本思想

前文提到，经典 Ewald 求和计算倒易空间部分能量的时间复杂度为 $\mathcal{O}(N^2)$，这是由离散 Fourier 变换（discrete Fourier transform，DFT）的基本性质决定的。对于较大体量的模拟系统，使用基本的离散 Fourier 变换处理的计算效率较低，如果可以将点电荷通过插值算法分布到空间中的均匀网格上，就取得了使用快速 Fourier 变换（fast Fourier transform，FFT）的条件，将时间复杂度降为 $\mathcal{O}(N{\log }_{2}N)$，从而极大地提升计算效率。PME 求和算法正是在这种思想的指导下被开发出来。

Lagrange 插值

为了近似静电结构因子 $S(\mathbf{k})$ ，需要对方程 $24$ 中出现的复指数 $\exp (2\pi i\mathbf{k}\cdot \mathbf{r})$ 进行插值。给定正整数 $K_1,K_2,K_3$ 和单元系统中的位置 $\mathbf{r}$，用 $u_1,u_2,u_3$ 表示其比例分数坐标（scaled fractional coordinates），其中 $u_{\alpha }=K_{\alpha }{\mathbf{a}}_{\alpha }^{\ast }\cdot \mathbf{r}$（$\alpha =1,2,3$）。由于周期性边界条件，可以假设 $0⩽u_{\alpha }<K_{\alpha }$。有 $$\exp \left(2\pi i\mathbf{k}\cdot \mathbf{r}\right)=\exp \left(2\pi i\frac{k_1{u}_1}{K_1}\right)\cdot \exp \left(2\pi i\frac{k_2{u}_2}{K_2}\right)\cdot \exp \left(2\pi i\frac{k_3{u}_3}{K_3}\right)\text{\hspace{1em}(31)}$$ 对于实数 $u$，设 $[u]$ 表示 $u$ 的整数部分，即 $[u]$ 是满足 $[u]⩽u⩽[u]+1$的唯一整数。可以使用线性插值（linear interpolation）对方程 $31$ 右侧的各个指数进行近似： $$\exp \left(2\pi i\frac{k_{\alpha }{u}_{\alpha }}{K_{\alpha }}\right)\approx (1-(u_{\alpha }-[u_{\alpha }]))\cdot \exp \left(2\pi i\frac{k_{\alpha }}{K_{\alpha }}[u_{\alpha }]\right)+(u_{\alpha }-[u_{\alpha }])\cdot \exp \left(2\pi i\frac{k_{\alpha }}{K_{\alpha }}([u_{\alpha }]+1)\right)\text{\hspace{1em}(32)}$$ 将方程 $32$（$\alpha =1,2,3$）分别代入方程 $31$，将方程 $31$ 右侧的乘积三线性插值（trilinear interpolation）为8项的和。用 $W_2(u)$ 表示以下线性帽函数（linear hat function）： $$W_2(u)=\{\begin{array}{ll}1-\left|u\right|& (-1⩽u⩽1)\\ 0& (\left|u\right|>1)\end{array}\text{\hspace{1em}(33)}$$ 可以将方程 $32$ 改写为 $$\exp \left(2\pi i\frac{k_{\alpha }{u}_{\alpha }}{K_{\alpha }}\right)\approx \sum _{m=-\infty }^{+\infty }{W}_2(u_{\alpha }-m)\cdot \exp \left(2\pi i\frac{k_{\alpha }}{K_{\alpha }}m\right)\text{\hspace{1em}(34)}$$ 其中 $m\in \mathbb{Z}$。注意到方程 $34$ 中的求和实际上是有限的，因为 $W_2$ 是有界支撑的（bounded support）。复指数的高阶分段 Lagrange 插值可以用类似的方式表示。考虑在点 $[u]-p+1,[u]-p+2,...,[u]+p$ 上的分段 $2p$ 阶 Lagrange 插值，定义分配函数： $$W_{2p}=\{\begin{array}{ll}\frac{{\prod }_{n=-p,n\ne m}^{p-1}(u+n-m)}{{\prod }_{n=-p,n\ne m}^{p-1}(n-m)}& (m⩽u⩽m+1,m=-p,-p+1,...,p-1)\\ 0& (\left|u\right|>p)\end{array}\text{\hspace{1em}(35)}$$ 注意到 $p=1$ 时，上式与 $W_2$ 的定义一致。利用该函数，可以写出复指数的一般分段 Lagrange 插值公式： $$\exp \left(2\pi i\frac{k_{\alpha }{u}_{\alpha }}{K_{\alpha }}\right)\approx \sum _{m=-\infty }^{+\infty }{W}_{2p}(u_{\alpha }-m)\cdot \exp \left(2\pi i\frac{k_{\alpha }}{K_{\alpha }}m\right)\text{\hspace{1em}(36)}$$ 根据 Lagrange 插值的标准误差估计，该近似中的误差有界于 $${\sigma }_{2p}=\frac{(2p)!}{(p!{)}^2}{\left(\frac{\pi \left|k_{\alpha }\right|}{2K_{\alpha }}\right)}^{2p}\text{\hspace{1em}(37)}$$

物理图像

下面讨论 Lagrange 插值的物理图像。为了简便性，考虑如图3所示的一维网格，假定网格宽度为 $1$，$\times$ 表示格点，$|$ 表示格子边界，即两个相邻格点的中点。

图3. 一维网格示意图。

假设 $P(P\in {\mathbb{Z}}_{+})$ 为 Lagrange 插值阶数。若 $P$ 为奇数，定义 $\bar{x}$ 为距离粒子最近的格点坐标；若 $P$ 为偶数，定义 $\bar{x}$ 为距离粒子最近的格点中点坐标。将粒子的电荷以分配函数 $W_i^{(P)}(x-\bar{x})$ 分配到以 $\bar{x}$ 为中心、索引分别为 $i=-\frac{P-1}{2},-\frac{P-1}{2}+1,...,\frac{P-1}{2}$ 的 $P$ 个格点上。下面给出 $P=3$ 和 $P=4$ 时的 Lagrange 插值分配函数 $W_i^{(P)}(x)$ 以供参考。

当 $P=3$ 时： $$\begin{array}{rl}& W_{-1}^{(3)}(x)=\frac{1}{8}(1-4x+4x^2)\\ & W_0^{(3)}(x)=\frac{1}{4}(3-4x^2)\\ & W_1^{(3)}(x)=\frac{1}{8}(1+4x+4x^2)\end{array}$$ 当 $P=4$ 时： $$\begin{array}{rl}& W_{-\frac{3}{2}}^{(4)}(x)=\frac{1}{48}(1-6x+12x^2-8x^3)\\ & W_{-\frac{1}{2}}^{(4)}(x)=\frac{1}{48}(23-30x-12x^2+24x^3)\\ & W_{\frac{1}{2}}^{(4)}(x)=\frac{1}{48}(23+30x-12x^2-24x^3)\\ & W_{\frac{3}{2}}^{(4)}(x)=\frac{1}{48}(1+6x+12x^2+8x^3)\end{array}$$ 分配函数容易扩展到三维空间： $$W(\mathbf{r})=W(x)W(y)W(z)\text{\hspace{1em}(38)}$$ 在周期性边界条件中，格点 $(m_1,m_2,m_3)$ 上的电荷量为 $$Q(m_1,m_2,m_3)=\sum _{i=1}^N\sum _{n_1,n_2,n_3}{q}_i\prod _{\alpha =1}^{3}W(u_{\alpha i}-m_{\alpha }-n_{\alpha }{K}_{\alpha })\text{\hspace{1em}(39)}$$

图4. 将点电荷分配到二维网格上（$P=4$）。蓝点表示初始点电荷，红点表示分配到网格上的电荷。圆点的大小指示电荷量。

本章小结

利用方程 $36$ 对静电结构因子 $S(\mathbf{k})$ 进行近似，有 $$\begin{array}{rl}S(\mathbf{k})\approx \tilde{S}(\mathbf{k})& =\sum _{i=1}^N{q}_i\sum _{m_1=-\infty }^{\infty }\sum _{m_2=-\infty }^{\infty }\sum _{m_3=-\infty }^{\infty }\prod _{\alpha =1}^{3}W(u_{\alpha i}-m_{\alpha })\exp \left(2\pi i\frac{k_{\alpha }}{K_{\alpha }}{m}_{\alpha }\right)\\ & =\sum _{m_1=0}^{K_1-1}\sum _{m_2=0}^{K_2-1}\sum _{m_3=0}^{K_3-1}Q(m_1,m_2,m_3)\exp \left(2\pi i\left(\frac{k_1{m}_1}{K_1}+\frac{k_2{m}_2}{K_2}+\frac{k_3{m}_3}{K_3}\right)\right)\\ & =F(Q)(m_1,m_2,m_3)\end{array}\text{\hspace{1em}(40)}$$ 其中 $F(Q)$ 表示数组 $Q$ 的离散 Fourier 变换。将方程 $40$ 代入方程 $27$，有 $$\begin{array}{rl}{E}^{\mathrm{recip}}\approx {\tilde{E}}^{\mathrm{recip}}& =\frac{1}{4\pi {\epsilon }_0}\frac{1}{2\pi V}\sum _{\mathbf{k}\ne 0}\frac{\exp (-2{\pi }^2{\sigma }^2{k}^2)}{k^2}F(Q)(m_1,m_2,m_3{)}^2\\ & =\frac{1}{4\pi {\epsilon }_0}\frac{1}{2\pi V}\sum _{k_1=0}^{K_1-1}\sum _{k_2=0}^{K_2-1}\sum _{k_3=0}^{K_3-1}\frac{\exp (-2{\pi }^2{\sigma }^2{k}^2)}{k^2}F(Q)(m_1,m_2,m_3{)}^2\end{array}\text{\hspace{1em}(41)}$$ 注意到对倒易格矢量 $\mathbf{k}$ 的分量 $k_1,k_2,k_3$ 的求和与对网格求和相同，这是由处理均匀取样信号的离散 Fourier 变换的性质决定的。

PME 求和算法实现

简介

• 使用C语言对 PME 求和计算倒易空间的经典相互作用能的方法进行代码实现。
• 代码调用 FFTW 库对网格上的分配电荷进行快速 Fourier 变换。
• 代码力求可读性，不具备计算效率。

示例代码

示例输出

Size: 17.150000 x 17.160000 x 56.000000 (Å)

Mesh: 42 x 42 x 140

Volume: 16.480463 (nm^3)

Atoms: 698

beta: 0.207441

Totol charge in the box: 4.164000e-06

Totol charge on the mesh: -5.851899e-03

Reciprocal energy = 15.606500 (eV)

Program ended with exit code: 0

参考文献

[1] Paul P. Ewald. "Die Berechnung optischer und elektrostatischer Gitterpotentiale." Annalen der physik 369.3 (1921): 253-287. (https://doi.org/10.1002/andp.19213690304)

[2] Tom Darden, Darrin York, and Lee Pedersen. "Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systems." The Journal of chemical physics 98.12 (1993): 10089-10092. (https://doi.org/10.1063/1.464397)

[3] Ulrich Essmann, Lalith Perera, Max L. Berkowitz, Tom Darden, Hsing Lee, and Lee G. Pedersen. "A smooth particle mesh Ewald method." The Journal of chemical physics 103.19 (1995): 8577-8593. (https://doi.org/10.1063/1.470117)

[4] R. A. Jackson, and C. R. A. Catlow. "Computer simulation studies of zeolite structure." Molecular Simulation 1.4 (1988): 207-224. (https://doi.org/10.1080/08927028808080944)