埃瓦爾德求和

埃瓦爾德求和（英語：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）

參考文獻

1. Frenkel, Daan; Smit, Berend. Understanding Molecular Simulation: From Algorithms to Applications. Academic Press. 2001-10-19. ISBN 9780080519982 （英語）.
2. Di Pierro, Michele; Elber, Ron; Leimkuhler, Benedict. A Stochastic Algorithm for the Isobaric–Isothermal Ensemble with Ewald Summations for All Long Range Forces. Journal of Chemical Theory and Computation. 2015-12-08, 11 (12): 5624–5637. ISSN 1549-9618. PMC 4890727 . PMID 26616351. doi:10.1021/acs.jctc.5b00648.
3. Kolafa, Jiri; Perram, John W. Cutoff Errors in the Ewald Summation Formulae for Point Charge Systems. Molecular Simulation. 1992-01-01, 9 (5): 351–368. ISSN 0892-7022. doi:10.1080/08927029208049126.
4. Darden, Tom; Perera, Lalith; Li, Leping; Pedersen, Lee. New tricks for modelers from the crystallography toolkit: the particle mesh Ewald algorithm and its use in nucleic acid simulations. Structure. 1999-03-15, 7 (3): R55–R60 [2017-07-03]. ISSN 0969-2126. doi:10.1016/s0969-2126(99)80033-1. （原始內容存檔於2022-01-22） （英語）.
5. Herce, Henry David; Garcia, Angel Enrique; Darden, Thomas. The electrostatic surface term: (I) Periodic systems. The Journal of Chemical Physics. 2007-03-26, 126 (12): 124106 [2017-07-03]. ISSN 0021-9606. doi:10.1063/1.2714527. （原始內容存檔於2020-07-11）.
6. Ewald, P. P. Die Berechnung optischer und elektrostatischer Gitterpotentiale. Annalen der Physik. 1921-01-01, 369 (3): 253–287 [2017-07-03]. ISSN 1521-3889. doi:10.1002/andp.19213690304. （原始內容存檔於2018-01-11） （英語）.
7. Darden, Tom; York, Darrin; Pedersen, Lee. Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systems. The Journal of Chemical Physics. 1993-06-15, 98 (12): 10089–10092 [2017-07-03]. ISSN 0021-9606. doi:10.1063/1.464397. （原始內容存檔於2020-12-08）.