埃瓦尔德求和

埃瓦尔德求和（英语：Ewald summation），是一种计算周期性系统（英语：Periodic systems of small molecules）中长程力（如静电力）的方法，以德国物理学家保罗·彼得·埃瓦尔德命名。埃瓦尔德求和最初用于计算离子晶体的电势能，现在用于计算化学中计算长程力。埃瓦尔德求和是泊松求和公式的特殊形式，用倒空间中的等效求和代替实空间中相互作用能（英语：Interaction energy）的总和。埃瓦尔德求和将相互作用势（英语：Interatomic potential）分为短程力和无奇点的长程力两部分，短程力在实空间中计算，长程力用傅里叶变换计算。与直接求和相比，此方法的优势为能量能够快速收敛，这意味着此方法在计算长程力时具有较高的精度和合理的速度，是计算周期性系统（英语：Periodic systems of small molecules）中长程力的标准方法。此方法需要分子系统的电中性，以准确计算总库仑力^[1]。

本条目中，矢量与标量分别用粗体与斜体显示。例如，位置矢量通常用 ${\displaystyle \mathbf {r} \,\!}$ 表示；而其大小则用 ${\displaystyle r\,\!}$ 来表示。

推导

埃瓦尔德求和将相互作用势（英语：Interatomic potential）表示为两部分之和：

${\displaystyle \varphi (\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \varphi _{sr}(\mathbf {r} )+\varphi _{\ell r}(\mathbf {r} )}$,

其中，${\displaystyle \varphi _{sr}(\mathbf {r} )}$表示实空间中和值快速收敛的短程势，${\displaystyle \varphi _{\ell r}(\mathbf {r} )}$表示倒空间中和值快速收敛的长程势。所有量（如r）的长程部分是有限的，但可能有简易的数学形式，如高斯分布。该方法假设短程势容易求和，因此需要重点考虑的是长程势。由于使用了傅里叶级数，该方法将周期性边界条件作为假设，此周期性系统的重复单元称为原胞，选择一个原胞作为中央原胞作为参考，其余单元称为镜像。

长程力的能量是中央原胞的电荷与晶格所有电荷间相互作用能（英语：Interaction energy）之和，因此可以表示为原胞和晶格的电荷密度的双重积分：

${\displaystyle E_{\ell r}=\iint d\mathbf {r} \,d\mathbf {r} ^{\prime }\,\rho _{\text{TOT}}(\mathbf {r} )\rho _{uc}(\mathbf {r} ^{\prime })\ \varphi _{\ell r}(\mathbf {r} -\mathbf {r} ^{\prime })}$

其中原胞的电荷密度${\displaystyle \rho _{uc}(\mathbf {r} )}$是中央原胞中位置${\displaystyle \mathbf {r} _{k}}$上的电量${\displaystyle q_{k}}$之和：

${\displaystyle \rho _{uc}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{\mathrm {charges} \ k}q_{k}\delta (\mathbf {r} -\mathbf {r} _{k})}$

总电荷密度${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}$是原胞及其镜像电量${\displaystyle q_{k}}$之和：

${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{n_{1},n_{2},n_{3}}\sum _{\mathrm {charges} \ k}q_{k}\delta (\mathbf {r} -\mathbf {r} _{k}-n_{1}\mathbf {a} _{1}-n_{2}\mathbf {a} _{2}-n_{3}\mathbf {a} _{3})}$

这里，${\displaystyle \delta (\mathbf {x} )}$表示狄拉克δ函数，${\displaystyle \mathbf {a} _{1}}$、${\displaystyle \mathbf {a} _{2}}$、${\displaystyle \mathbf {a} _{3}}$表示晶格矢量，${\displaystyle n_{1}}$、${\displaystyle n_{2}}$、${\displaystyle n_{3}}$的范围为所有整数。总电荷密度${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}$可以表示为${\displaystyle \rho _{uc}(\mathbf {r} )}$与晶格函数${\displaystyle L(\mathbf {r} )}$的卷积：

${\displaystyle L(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{n_{1},n_{2},n_{3}}\delta (\mathbf {r} -n_{1}\mathbf {a} _{1}-n_{2}\mathbf {a} _{2}-n_{3}\mathbf {a} _{3})}$

由于${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}$为卷积，其傅里叶变换为一个积：

${\displaystyle {\tilde {\rho }}_{\text{TOT}}(\mathbf {k} )={\tilde {L}}(\mathbf {k} ){\tilde {\rho }}_{uc}(\mathbf {k} )}$

其中晶格函数的傅里叶变换是狄拉克δ函数的另一个和：

${\displaystyle {\tilde {L}}(\mathbf {k} )={\frac {\left(2\pi \right)^{3}}{\Omega }}\sum _{m_{1},m_{2},m_{3}}\delta (\mathbf {k} -m_{1}\mathbf {b} _{1}-m_{2}\mathbf {b} _{2}-m_{3}\mathbf {b} _{3})}$

其中定义倒空间矢量为${\displaystyle \mathbf {b} _{1}\ {\stackrel {\mathrm {def} }{=}}\ {\frac {\mathbf {a} _{2}\times \mathbf {a} _{3}}{\Omega }}}$（周期性排列），其中${\displaystyle \Omega \ {\stackrel {\mathrm {def} }{=}}\ \mathbf {a} _{1}\cdot \left(\mathbf {a} _{2}\times \mathbf {a} _{3}\right)}$为中心原胞的体积（几何形状通常为平行六面体），${\displaystyle L(\mathbf {r} )}$和${\displaystyle {\tilde {L}}(\mathbf {k} )}$为实函数和偶函数。

为了简洁起见，定义有效单粒子势能：

${\displaystyle v(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \int d\mathbf {r} ^{\prime }\,\rho _{uc}(\mathbf {r} ^{\prime })\ \varphi _{\ell r}(\mathbf {r} -\mathbf {r} ^{\prime })}$

因为其亦为卷积，其傅里叶变换是一个积：

${\displaystyle {\tilde {V}}(\mathbf {k} )\ {\stackrel {\mathrm {def} }{=}}\ {\tilde {\rho }}_{uc}(\mathbf {k} ){\tilde {\Phi }}(\mathbf {k} )}$

其中定义了傅里叶变换：

${\displaystyle {\tilde {V}}(\mathbf {k} )=\int d\mathbf {r} \ v(\mathbf {r} )\ e^{-i\mathbf {k} \cdot \mathbf {r} }}$

现在，长程力的能量可以表示为单个电荷密度的积分：

${\displaystyle E_{\ell r}=\int d\mathbf {r} \ \rho _{\text{TOT}}(\mathbf {r} )\ v(\mathbf {r} )}$

使用帕塞瓦尔定理，能量亦可于倒空间中求和：

${\displaystyle E_{\ell r}=\int {\frac {d\mathbf {k} }{\left(2\pi \right)^{3}}}\ {\tilde {\rho }}_{\text{TOT}}^{*}(\mathbf {k} ){\tilde {V}}(\mathbf {k} )=\int {\frac {d\mathbf {k} }{\left(2\pi \right)^{3}}}{\tilde {L}}^{*}(\mathbf {k} )\left|{\tilde {\rho }}_{uc}(\mathbf {k} )\right|^{2}{\tilde {\Phi }}(\mathbf {k} )={\frac {1}{\Omega }}\sum _{m_{1},m_{2},m_{3}}\left|{\tilde {\rho }}_{uc}(\mathbf {k} )\right|^{2}{\tilde {\Phi }}(\mathbf {k} )}$

其中${\displaystyle \mathbf {k} =m_{1}\mathbf {b} _{1}+m_{2}\mathbf {b} _{2}+m_{3}\mathbf {b} _{3}}$是最终的和值。

计算出${\displaystyle {\tilde {\rho }}_{uc}(\mathbf {k} )}$后，${\displaystyle \mathbf {k} }$的和值或积分是显然的，可以很快地收敛。不能收敛的最常见原因是原胞不太明确，其必须为电中性，以避免无穷大的和。

粒子网格埃瓦尔德（PME）方法

在计算机普及前，埃瓦尔德求和是理论物理的理论。然而，自20世纪70年代以来，埃瓦尔德求和在粒子系统的计算机模拟中被广泛使用，尤其是遵守平方反比定律的粒子相互作用，如重力和静电力。最近，粒子网格埃瓦尔德方法也用于计算兰纳-琼斯势的${\displaystyle r^{-6}}$部分，以消除截断产生的伪影（英语：Artifact (error)）^[2][3]。其应用包括等离子体、星系及分子的模拟^[4]。

在粒子网格埃瓦尔德方法中，和标准埃瓦尔德求和相同，相互作用势（英语：Interatomic potential）被分为两部分${\displaystyle \varphi (\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \varphi _{sr}(\mathbf {r} )+\varphi _{\ell r}(\mathbf {r} )}$，其基本思想是用实空间中短程力的直接求和${\displaystyle E_{sr}}$（粒子部分），及倒空间中长程力的求和（埃瓦尔德部分），代替点粒子间相互作用的能量的直接求和：

${\displaystyle E_{\text{TOT}}=\sum _{i,j}\varphi (\mathbf {r} _{j}-\mathbf {r} _{i})=E_{sr}+E_{\ell r}}$

${\displaystyle E_{sr}=\sum _{i,j}\varphi _{sr}(\mathbf {r} _{j}-\mathbf {r} _{i})}$

${\displaystyle E_{\ell r}=\sum _{\mathbf {k} }{\tilde {\Phi }}_{\ell r}(\mathbf {k} )\left|{\tilde {\rho }}(\mathbf {k} )\right|^{2}}$

其中${\displaystyle {\tilde {\Phi }}_{\ell r}}$和${\displaystyle {\tilde {\rho }}(\mathbf {k} )}$表示力和电荷密度的傅里叶变换。由于两个求和分别在实空间和倒空间中迅速收敛，它们可能被精确截断，且所需计算时间大幅减少。计算电荷密度的傅里叶变换${\displaystyle {\tilde {\rho }}(\mathbf {k} )}$可使用快速傅里叶变换，需在空间中的离散格子（英语：Lattice (discrete subgroup)）上（即网格部分）估计电荷密度。

由于埃瓦尔德方法隐含的周期性假设，粒子网格埃瓦尔德方法于物理系统中的应用需施加周期性。因此，该方法最适合用于空间范围内可以模拟为无限的系统。在分子动力学模拟中，常构造可以无限平铺形成镜像的电中性原胞；然而，为了正确解释这种近似效应，这些镜像被重新并入原始模拟原胞中，这种整体效应被称为周期性边界条件。 想象一个单位立方体，上表面与下表面有效接触，右侧面与左侧面有效接触，前表面与后表面有效接触。因此，原胞的尺寸必须足够大，以避免两个接触面间不正确的运动相关性，但仍需足够小以便计算。短程力与长程力间截断的定义也可以引入伪影（英语：Artifact (error)）。

电荷密度对网格的限制，使得粒子网格埃瓦尔德方法对电荷密度或势函数平滑变化的系统更有效。利用快速多极子方法（英语：Fast multipole method）可以更有效地处理局部系统或电荷密度波动较大的系统。

偶极子

极性晶体（即原胞中具有净偶极子${\displaystyle \mathbf {p} _{uc}}$的晶体）的静电能为条件收敛，即取决于求和顺序。例如，若中央原胞的偶极与不断增加的立方体上的原胞偶极相互作用，则其能量收敛值并不会与考虑不断增大的球面时相等。大致来说，这种条件收敛是因为在半径为${\displaystyle R}$的壳上的偶极子数约为${\displaystyle R^{2}}$；偶极-偶极相互作用的强度约为${\displaystyle {\frac {1}{R^{3}}}}$；而两者相乘的结果是发散的调和级数${\displaystyle \sum _{n=1}^{\infty }{\frac {1}{n}}}$。

这看似令人惊讶的结果并不与现实晶体能量有限的事实相违背，因为现实晶体并非无限，具有特定边界。具体而言，极性晶体的边界的有效表面电荷密度为${\displaystyle \sigma =\mathbf {P} \cdot \mathbf {n} }$，其中${\displaystyle \mathbf {n} }$为表面法矢量，${\displaystyle \mathbf {P} }$为单位体积的净偶极矩。则中央原胞之偶极子与表面电荷密度${\displaystyle \sigma }$的相互作用能${\displaystyle U}$可写为^[5]：

${\displaystyle U={\frac {1}{2V_{uc}}}\int {\frac {\left(\mathbf {p} _{uc}\cdot \mathbf {r} \right)\left(\mathbf {p} _{uc}\cdot \mathbf {n} \right)dS}{r^{3}}}}$

其中，${\displaystyle \mathbf {p} _{uc}}$和${\displaystyle V_{uc}}$分别为原胞的净偶极矩和体积，${\displaystyle dS}$为晶面上的无穷小区域，${\displaystyle \mathbf {r} }$为中央原胞到无穷小区域的矢量。此公式来自于对能量${\displaystyle dU=-\mathbf {p} _{uc}\cdot \mathbf {dE} }$积分，其中${\displaystyle d\mathbf {E} }$表示无穷小电场，由无穷小的表面电荷${\displaystyle dq{\stackrel {\mathrm {def} }{=}}\sigma dS}$产生（库仑定律）：

${\displaystyle d\mathbf {E} \ {\stackrel {\mathrm {def} }{=}}\ \left({\frac {-1}{4\pi \epsilon }}\right){\frac {dq\ \mathbf {r} }{r^{3}}}=\left({\frac {-1}{4\pi \epsilon }}\right){\frac {\sigma \,dS\ \mathbf {r} }{r^{3}}}}$

负号来自于${\displaystyle \mathbf {r} }$的定义，其指向电荷方向为正方向。

历史

埃瓦尔德求和由德国物理学家保罗·彼得·埃瓦尔德于1921年发表，用于确定离子晶体的静电能及马德隆常数^[6]。

复杂度

不同的埃瓦尔德求和具有不同的时间复杂度。直接求和的时间复杂度为${\displaystyle O(N^{2})}$，其中${\displaystyle N}$为系统中原子数。粒子网格埃瓦尔德方法的时间复杂度为${\displaystyle O(N\,\log N)}$^[7]。

参见

• 保罗·彼得·埃瓦尔德
• 马德隆常数
• 泊松求和公式
• 分子建模
• Wolf求和（英语：Wolf summation）