查看“︁埃瓦尔德求和”︁的源代码

{{NoteTA|G1=化学|G2=物理学}}
{{向量字體常規}}
'''埃瓦尔德求和'''（{{lang-en|Ewald summation}}），是一种计算{{link-en|分子的周期性系统|Periodic systems of small molecules|周期性系统}}中长程力（如[[静电力]]）的方法，以德国物理学家[[保罗·彼得·埃瓦尔德]]命名。埃瓦尔德求和最初用于计算[[离子晶体]]的[[电势能]]，现在用于[[计算化学]]中计算长程力。埃瓦尔德求和是[[泊松求和公式]]的特殊形式，用[[倒易点阵|倒空间]]中的等效求和代替[[位置空间与动量空间|实空间]]中{{link-en|相互作用能|Interaction energy}}的总和。埃瓦尔德求和将{{link-en|相互作用势|Interatomic potential}}分为短程力和无[[奇点 (数学)|奇点]]的长程力两部分，短程力在[[位置空间与动量空间|实空间]]中计算，长程力用[[傅里叶变换]]计算。与直接求和相比，此方法的优势为能量能够快速[[收敛]]，这意味着此方法在计算长程力时具有较高的精度和合理的速度，是计算{{link-en|分子的周期性系统|Periodic systems of small molecules|周期性系统}}中长程力的标准方法。此方法需要分子系统的电中性，以准确计算总[[库仑力]]<ref>{{Cite book|url=https://books.google.com.hk/books?hl=zh-CN&lr=&id=5qTzldS9ROIC&oi=fnd&pg=PP2&dq=Understanding+molecular+simulation:+from+algorithms+to+applications&ots=nFTQUsXaXf&sig=S4_yA-K1hPor-XJf1Xjz2BpZ52I&redir_esc=y&hl=zh-CN&sourceid=cndr#v=onepage&q=Understanding%20molecular%20simulation:%20from%20algorithms%20to%20applications&f=false|last=Frenkel|first=Daan|last2=Smit|first2=Berend|date=2001-10-19|publisher=Academic Press|isbn=9780080519982|language=en|title=Understanding Molecular Simulation: From Algorithms to Applications}}</ref>。

==推导==

埃瓦尔德求和将{{link-en|相互作用势|Interatomic potential}}表示为两部分之和：

:<math>\varphi(\mathbf{r}) \ \stackrel{\mathrm{def}}{=}\  \varphi_{sr}(\mathbf{r}) + \varphi_{\ell r}(\mathbf{r})</math>,

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

长程力的[[能量]]是中央原胞的[[电荷]]与晶格所有电荷间{{link-en|相互作用能|Interaction energy}}之和，因此可以表示为[[原胞]]和[[晶格]]的[[电荷密度]]的[[双重积分]]：

:<math>
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)
</math>

其中[[原胞]]的[[电荷密度]]<math>\rho_{uc}(\mathbf{r})</math>是中央原胞中位置<math>\mathbf{r}_k</math>上的[[电量]]<math>q_k</math>之和：

:<math>
\rho_{uc}(\mathbf{r}) \ \stackrel{\mathrm{def}}{=}\  \sum_{\mathrm{charges}\ k} q_k \delta(\mathbf{r} - \mathbf{r}_k)
</math>

总[[电荷密度]]<math>\rho_\text{TOT}(\mathbf{r})</math>是[[原胞]]及其[[镜像 (几何)|镜像]]电量<math>q_{k}</math>之和：

:<math>
\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)
</math>

这里，<math>\delta(\mathbf{x})</math>表示[[狄拉克δ函数]]，<math>\mathbf{a}_1</math>、<math>\mathbf{a}_2</math>、<math>\mathbf{a}_3</math>表示晶格矢量，<math>n_1</math>、<math>n_2</math>、<math>n_3</math>的范围为所有整数。总电荷密度<math>\rho_\text{TOT}(\mathbf{r})</math>可以表示为<math>\rho_{uc}(\mathbf{r})</math>与晶格函数<math>L(\mathbf{r})</math>的[[卷积]]：

:<math>
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)
</math>

由于<math>\rho_\text{TOT}(\mathbf{r})</math>为[[卷积]]，其[[傅里叶变换]]为一个[[积]]：

:<math>
\tilde{\rho}_\text{TOT}(\mathbf{k}) = \tilde{L}(\mathbf{k}) \tilde{\rho}_{uc}(\mathbf{k})
</math>

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

:<math>
\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)
</math>

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

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

:<math>
v(\mathbf{r}) \ \stackrel{\mathrm{def}}{=}\  \int d\mathbf{r}^{\prime}\, \rho_{uc}(\mathbf{r}^\prime) \ \varphi_{\ell r}(\mathbf{r} - \mathbf{r}^\prime)
</math>

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

:<math>
\tilde{V}(\mathbf{k}) \ \stackrel{\mathrm{def}}{=}\  \tilde{\rho}_{uc}(\mathbf{k}) \tilde{\Phi}(\mathbf{k})
</math>

其中定义了[[傅里叶变换]]：

:<math>
\tilde{V}(\mathbf{k}) = \int d\mathbf{r} \ v(\mathbf{r}) \ e^{-i\mathbf{k} \cdot \mathbf{r}}
</math>

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

:<math>
E_{\ell r} = \int d\mathbf{r} \ \rho_\text{TOT}(\mathbf{r}) \ v(\mathbf{r}) 
</math>

使用[[帕塞瓦尔定理]]，能量亦可于[[倒易点阵|倒空间]]中求和：

:<math>
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})
</math>

其中<math>\mathbf{k} = m_1 \mathbf{b}_1 + m_2 \mathbf{b}_2 + m_3 \mathbf{b}_3</math>是最终的和值。

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

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

在[[电子计算机|计算机]]普及前，埃瓦尔德求和是[[理论物理]]的理论。然而，自20世纪70年代以来，埃瓦尔德求和在粒子系统的[[计算机模拟]]中被广泛使用，尤其是遵守[[平方反比定律]]的粒子相互作用，如[[重力]]和[[静电力]]。最近，粒子网格埃瓦尔德方法也用于计算[[兰纳-琼斯势]]的<math>r^{-6}</math>部分，以消除{{link-en|截断|Truncation}}产生的{{link-en|伪影|Artifact (error)}}<ref>{{Cite journal|title=A Stochastic Algorithm for the Isobaric–Isothermal Ensemble with Ewald Summations for All Long Range Forces|url=http://dx.doi.org/10.1021/acs.jctc.5b00648|last=Di Pierro|first=Michele|last2=Elber|first2=Ron|date=2015-12-08|journal=Journal of Chemical Theory and Computation|issue=12|doi=10.1021/acs.jctc.5b00648|volume=11|pages=5624–5637|issn=1549-9618|pmc=4890727|pmid=26616351|last3=Leimkuhler|first3=Benedict}}</ref><ref>{{Cite journal|title=Cutoff Errors in the Ewald Summation Formulae for Point Charge Systems|url=http://dx.doi.org/10.1080/08927029208049126|last=Kolafa|first=Jiri|last2=Perram|first2=John W.|date=1992-01-01|journal=Molecular Simulation|issue=5|doi=10.1080/08927029208049126|volume=9|pages=351–368|issn=0892-7022}}</ref>。其应用包括[[等离子体]]、[[星系]]及[[分子]]的[[模拟]]<ref>{{Cite journal|title=New tricks for modelers from the crystallography toolkit: the particle mesh Ewald algorithm and its use in nucleic acid simulations|url=http://linkinghub.elsevier.com/retrieve/pii/S0969212699800331|last=Darden|first=Tom|last2=Perera|first2=Lalith|date=1999-03-15|journal=Structure|issue=3|doi=10.1016/s0969-2126(99)80033-1|volume=7|pages=R55–R60|language=en|issn=0969-2126|last3=Li|first3=Leping|last4=Pedersen|first4=Lee|access-date=2017-07-03|archive-date=2022-01-22|archive-url=https://web.archive.org/web/20220122075419/https://linkinghub.elsevier.com/retrieve/pii/S0969212699800331|dead-url=no}}</ref>。

在粒子网格埃瓦尔德方法中，和标准埃瓦尔德求和相同，{{link-en|相互作用势|Interatomic potential}}被分为两部分<math>\varphi(\mathbf{r}) \ \stackrel{\mathrm{def}}{=}\  \varphi_{sr}(\mathbf{r}) + \varphi_{\ell r}(\mathbf{r})</math>，其基本思想是用[[位置空间与动量空间|实空间]]中短程力的直接求和<math>E_{sr}</math>（粒子部分），及[[倒易点阵|倒空间]]中长程力的求和（埃瓦尔德部分），代替点粒子间相互作用的能量的直接求和：

:<math>
E_\text{TOT} = \sum_{i,j} \varphi(\mathbf{r}_{j} - \mathbf{r}_i) = E_{sr} + E_{\ell r}
</math>

:<math>
E_{sr} = \sum_{i,j} \varphi_{sr}(\mathbf{r}_j - \mathbf{r}_i)
</math>

:<math>
E_{\ell r} = \sum_{\mathbf{k}} \tilde{\Phi}_{\ell r}(\mathbf{k}) \left| \tilde{\rho}(\mathbf{k}) \right|^2
</math>

其中<math>\tilde{\Phi}_{\ell r}</math>和<math>\tilde{\rho}(\mathbf{k})</math>表示[[力]]和[[电荷密度]]的[[傅里叶变换]]。由于两个求和分别在[[位置空间与动量空间|实空间]]和[[倒易点阵|倒空间]]中迅速[[收敛]]，它们可能被精确{{link-en|截断|Truncation}}，且所需计算时间大幅减少。计算[[电荷密度]]的[[傅里叶变换]]<math>\tilde{\rho}(\mathbf{k})</math>可使用[[快速傅里叶变换]]，需在空间中的{{link-en|离散格子|Lattice (discrete subgroup)}}上（即网格部分）估计[[电荷密度]]。

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

[[电荷密度]]对网格的限制，使得粒子网格埃瓦尔德方法对[[电荷密度]]或[[向量勢|势函数]]平滑变化的系统更有效。利用{{link-en|快速多极子方法|Fast multipole method}}可以更有效地处理[[局部系统]]或[[电荷密度]]波动较大的系统。

==偶极子==
[[极性]][[晶体]]（即[[原胞]]中具有净[[偶极子]]<math>\mathbf{p}_{uc}</math>的晶体）的[[静电]]能為[[条件收敛]]，即取决于求和顺序。例如，若中央原胞的偶极与不断增加的[[立方体]]上的原胞偶极相互作用，则其能量收敛值並不會与考慮不斷增大的球面時相等。大致来说，这种条件收敛是因为在半径为<math>R</math>的壳上的偶极子数約為<math>R^{2}</math>；偶极-偶极相互作用的强度約為<math>\frac{1}{R^{3}}</math>；而兩者相乘的結果是發散的[[调和级数]]<math>\sum_{n=1}^{\infty} \frac{1}{n}</math>。

這看似令人驚訝的結果並不與現實晶體能量有限的事實相違背，因為現實晶體並非無限，具有特定邊界。具体而言，极性晶体的边界的有效[[表面电荷]]密度为<math>\sigma = \mathbf{P} \cdot \mathbf{n}</math>，其中<math>\mathbf{n}</math>为表面[[法向量]]，<math>\mathbf{P}</math>为单位体积的净[[偶极矩]]。則中央原胞之偶极子與表面电荷密度<math>\sigma</math>的相互作用能<math>U</math>可寫為<ref>{{Cite journal|title=The electrostatic surface term: (I) Periodic systems|url=http://aip.scitation.org/doi/10.1063/1.2714527|last=Herce|first=Henry David|last2=Garcia|first2=Angel Enrique|date=2007-03-26|journal=The Journal of Chemical Physics|issue=12|doi=10.1063/1.2714527|volume=126|pages=124106|issn=0021-9606|last3=Darden|first3=Thomas|access-date=2017-07-03|archive-date=2020-07-11|archive-url=https://web.archive.org/web/20200711121201/https://aip.scitation.org/doi/10.1063/1.2714527|dead-url=no}}</ref>：

:<math>
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}
</math>

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

:<math>
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}
</math>
负号来自于<math>\mathbf{r}</math>的定义，其指向电荷方向为正方向。

==历史==

埃瓦尔德求和由德国物理学家[[保罗·彼得·埃瓦尔德]]于1921年发表，用于确定[[离子晶体]]的静电能及[[马德隆常数]]<ref>{{Cite journal|title=Die Berechnung optischer und elektrostatischer Gitterpotentiale|url=http://onlinelibrary.wiley.com/doi/10.1002/andp.19213690304/abstract|last=Ewald|first=P. P.|date=1921-01-01|journal=Annalen der Physik|issue=3|doi=10.1002/andp.19213690304|volume=369|pages=253–287|language=en|issn=1521-3889|access-date=2017-07-03|archive-date=2018-01-11|archive-url=https://web.archive.org/web/20180111102639/http://onlinelibrary.wiley.com/doi/10.1002/andp.19213690304/abstract|dead-url=no}}</ref>。

==复杂度==
不同的埃瓦尔德求和具有不同的[[时间复杂度]]。直接求和的时间复杂度为<math>O(N^2)</math>，其中<math>N</math>为系统中原子数。粒子网格埃瓦尔德方法的时间复杂度为<math>O(N\,\log N)</math><ref name="Darden1993">{{Cite journal|title=Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systems|url=http://aip.scitation.org/doi/10.1063/1.464397|last=Darden|first=Tom|last2=York|first2=Darrin|date=1993-06-15|journal=The Journal of Chemical Physics|issue=12|doi=10.1063/1.464397|volume=98|pages=10089–10092|issn=0021-9606|last3=Pedersen|first3=Lee|access-date=2017-07-03|archive-date=2020-12-08|archive-url=https://web.archive.org/web/20201208150524/https://aip.scitation.org/doi/10.1063/1.464397|dead-url=no}}</ref>。

==参见==

* [[保罗·彼得·埃瓦尔德]]
* [[馬德隆常數]]
* [[泊松求和公式]]
* [[分子建模]]
* {{link-en|Wolf求和|Wolf summation}}

==参考文献==
{{Reflist}}
[[Category:电磁学]]
[[Category:位势论]]
[[Category:傅里叶分析]]
[[Category:计算化学]]
[[Category:理论化学]]