技术背景
分子动力学模拟中，计算周期性边界条件的静电势常被视作计算的瓶颈之一。形式上是比较容易的，例如不考虑周期性边界条件的话，静电势能就是：
E=14πϵ0N−2∑i=0N−1∑j=i+1qiqjrij
如果考虑周期性边界条件，那么静电势能变为：
E=14πϵ0∑nN−2∑i=0N−1∑j=i+1qiqj∣∣rij+nL∣∣+14πϵ0∑|n|>0q2i|nL|,n=(n1,n2,...,nd)∈Z
且不说对无穷个盒子的叠加，就算是单个的盒子，也是O(N2)
的计算复杂度，这也是求解困难的原因。在分子动力学中，为了简化这个计算量，使用了三种思想：傅里叶变换、Ewald Summation以及Particle Mesh的方法，本文主要涉及傅里叶变换与Ewald求和计算的部分。
周期性电势
首先我们从物理概念上理解静电势能项的含义：单一的电荷qi
会在空间中形成一个电势Vi，当另一个电荷qj位于Vi对应的电场中时，就会受到静电相互作用，其能量为Eij。因为无穷远处的静电势能为0，这也就以为着，如果我们需要将qj从原始的位置推到无穷远之外，就需要Eij的能量。这个思路告诉我们，可以先计算电势Vi，再计算电势能Eij
，也就是这个东西：
Vi(r)=14πϵ0qi|r−ri|
如果考虑上周期性边界条件就是：
Vi(r)=14πϵ0∑nqi|r−ri+nL|
因为这个无穷多的求和没办法直接计算，只能明确电势具有周期性：Vi(r)=Vi(r+nL)
。
静电场泊松方程
空间中的电荷i
在r
处的电势，由泊松方程给出：
ΔVi(r)=−ρi(r)ϵ0
其中Δ
是拉普拉斯算子，ρi表示电荷密度。如果考虑一个点电荷的的情况，那么就有：ρi(r)=qiδ(r−ri)
。进而写出欧几里得空间中的泊松方程：
∇2Vi(r)=−qiδ(r−ri)ϵ0
其中∇2Vi(r)=∂2Vi(r)∂x2+∂2Vi(r)∂y2+∂2Vi(r)∂z2
。再考虑周期性边界条件，每个盒子中都有一个点电荷qi
，于是方程应该写为：
∇2Vi(r)=−∑nqiδ(r−ri+nL)ϵ0
傅里叶空间泊松方程
对上述单点形式的泊松方程的两边同时进行傅里叶变换（关于傅里叶变换的理解，可以参考前序文章1和前序文章2，有比较详细的原理介绍和相关代码实现）有：
∫(∂2Vi(r)∂x2+∂2Vi(r)∂y2+∂2Vi(r)∂z2)e−2πjkrdr=−qiϵ0∫δ(r−ri)e−2πjk⋅rdr
先使用分部积分计算左边中的一项：
∫∂2Vi(r)∂x2e−2πjk⋅rdr=∂Vi(r)∂xe−2πjk⋅r∣∣∣∞−∞+2πjkx∫∂Vi(r)∂xe−2πjk⋅rdr
需要注意的是，这里∂Vi(r)∂x=Fi(r)
的物理意义是作用力，那么我们就可以取第一类边界条件：
limx→∞Vi(x)=0
这样根据微分的定义有：
limx→∞∂Vi(x)∂x=limx→∞limϵ→0Vi(x+ϵ)−Vi(x)ϵ=0
上面这个极限代入了Vi(r)
的单点形式，其意义为，位于ri处的点电荷qi，在无穷远处生成的电场为0，对无穷远处的电荷q∞
的作用力也是0，这里不考虑周期性边界条件。则有：
∫∂2Vi(r)∂x2e−2πjk⋅rdr=2πjkx∫∂Vi(r)∂xe−2πjk⋅rdr=2πjkx[Vi(r)e−2πjk⋅r∣∣x=∞x=−∞+2πjkx∫Vi(r)e−2πjk⋅rdr]=−4π2k2xVi(k)
同理可得泊松方程左侧形式为：
∫∇2Vi(r)e−2πjk⋅rdr=−4π2k2Vi(k)
而右侧形式需要用到狄拉克函数的抽样特性：
∫∞−∞δ(t)f(t)dt=∫∞−∞δ(t)f(0)dt=f(0)∫∞−∞δ(t)dt=f(0)
即：
∫δ(r−ri)e−2πjk⋅rdr=e−2πjk⋅ri
得到傅里叶空间的单点泊松方程：
4π2k2Vi(k)=qiϵ0e−2πjk⋅ri
倒易空间
涉及到傅里叶空间，我们很自然的想到使用固体物理学的倒易空间变换，也就是把周期性盒子当作一个原胞。根据倒易空间（也叫k
空间）晶格矢（倒格矢）定义有：
k=m1k1+m2k2+m3k3
其中：
ki⋅Lj=δi,jδi,j={1,i=j0,i≠j
按照我们常用的长方体周期性边界条件：
L=(Lx,Ly,Lz)=L1+L2+L3L1=(Lx,0,0)L2=(0,Ly,0)L3=(0,0,Lz)
可以计算得：
k1=2π(L2×L3)L1⋅(L2×L3)=2πΩ(L2×L3)=πLyLzΩLxL1=2πL2xL1
其中Ω
表示周期性盒子的体积，类似的有：
k2=2π(L3×L1)L2⋅(L3×L1)=2πL2yL2
以及
k3=2π(L1×L2)L3⋅(L1×L2)=2πL2zL3
经过倒易空间变换之后，原胞体积从Ω=LxLyLz
变成：Ω∗=(2π)3LxLyLz。因为在前一步傅里叶空间的泊松方程中我们注意到k前面总是带了一个2π，这里不妨使用倒易晶格矢的定义对k
的形式做一个简化：
k=(2πLx,2πLy,2πLz)
这样一来傅里叶空间的泊松方程可以简写为：
Vi(k)=qiϵ0k2e−jk⋅ri
其中k2=|k|2=4π2m21L2x+4π2m22L2y+4π2m23L2z
，可以实现实空间到倒易空间的变换。
衰减函数构造
对于上述傅里叶变换之后的单点电势的形式，即使我们对整个k
空间进行积分，也是一个发散的结果。所以这里用到了一个非常特别的思想，由Edwald提出，把静电能量项分为远程相互作用项和短程相互作用项，分别在倒易空间和实空间收敛，这样就可以精确计算静电能。实际操作的时候有不同的推导过程，我们这里引用一种比较“数学”的推导方法(参考链接1)。
首先我们构造一个衰减函数e−k2t
，这个衰减函数有个特性是：
∫∞0e−k2tdt=(−1k2e−k2t)∣∣∣∞0=1k2
这样我们就可以用这个积分形式替换掉傅里叶-泊松方程中的1k2
项：
Vi(k)=qiϵ0e−jk⋅ri∫∞0e−k2tdt
因为这里使用的是从0到无穷大的一个积分形式，那么我们就可以实现一个截断，将其划分成两个积分的加和，假如我们在η
处做一个截断，则有：
Vi(k)=qiϵ0e−jk⋅ri(∫η0e−k2tdt+∫∞ηe−k2tdt)
这里取短程(Short Term)相互作用为：
VSi(k)=qiϵ0e−jk⋅ri∫η0e−k2tdt
以及长程(Long Term)相互作用为：
VLi(k)=qiϵ0e−jk⋅ri∫∞ηe−k2tdt=qiϵ0k2e−jk⋅rie−ηk2
短程作用项计算
按照预期，划分了短程作用项和长程作用项之后，应该可以得到一个实空间收敛的短程相互作用，我们对短程作用做一个逆傅里叶变换：
VSi(r)=1kxkykz∑kVSi(k)eik⋅r=qikxkykzϵ0∑kejk⋅(r−ri)∫η0e−k2tdt=qikxkykzϵ0∫η0∑kejk⋅(r−ri)e−k2tdt
很明显，积分内的求和项是一个指数平方函数的离散傅里叶变换。而我们可以知道，正态分布函数f(ξ)=1√2πσe−ξ22σ2
的傅里叶变换和逆傅里叶变换不改变其分布形式：
F(k)=∫f(ξ)e−jkξdξ=1√2πσ∫e−ξ2−2jkξσ22σ2dξ=1√2πσ∫e−ξ2−2jkξσ2−(jkσ2)2+(jkσ2)22σ2dξ=1√2πσe(jkσ2)22σ2∫e−(ξ+jkσ2)22σ2dξ
由于积分项只是一个实空间的积分，其本质还是一个正态分布函数的积分，我们知道其积分结果是一个常数，所以有：
F(k)=e−k2σ22
也是一个正态分布，只是其均方差从σ
变成了1σ
，也就是其积分结果为：
∫F(k)dk=√2πσ
同样的道理我们也可以计算得，正态分布函数得逆傅里叶变换结果也依然是一个正态分布函数，其均方差也是1σ
。
那么回到短程静电势，先做个变量替换t=σ22,σ≥0，则有dt=σdσ,σt=η=√2η,σt=0=0
。此时短程相互作用势为：
VSi(r)=qikxkykzϵ0∫√2η0(∑kejk⋅(r−ri)e−k2σ22)σdσ=qiϵ0∫√2η0e−(r−ri)22σ2Ncoeσdσ
这里Ncoe
是一个用于归一化正态分布的常数：
Ncoe=∫re−∣∣r−ri∣∣2σ2dr=∫∞−∞e−(z−zi)22σ2[∫∞−∞e−(y−yi)22σ2(∫∞−∞e−(x−xi)22σ2dx)dy]dz=(2πσ2)32
所以有：
VSi(r)=qiϵ0∫√2η0e−(r−ri)22σ2(2πσ2)32σdσ=qi(2π)32ϵ0∫√2η0e−(r−ri)22σ2σ2dσ
这里用一个变量替换：y=−|r−ri|√2σ
，则有：σ=−|r−ri|√2y，其微分变换形式为：dσ=|r−ri|√2y2dy，其边界为：yσ=√2η=−|r−ri|2√η,yσ=0=−∞
，代入得：
VSi(r)=qi(2π)32ϵ0∫−|r−ri|2√η−∞√2e−y2|r−ri|dy=qi2π32ϵ0|r−ri|∫−|r−ri|2√η−∞e−y2dy=qi2π32ϵ0|r−ri|∫+∞|r−ri|2√ηe−y2dy
此时使用一个误差函数Erfc(y)=2√π∫∞ye−x2dx
代入进行替换：
VSi(r)=qi4πϵ0|r−ri|Erfc(|r−ri|2√η)
因为η
是我们手动引入的一个常数参量，如果考虑η=σ22
，那么形式就变成了：
VSi(r)=qi4πϵ0|r−ri|Erfc(|r−ri|√2σ)
这个形式的短程相互作用势，表示的是单个盒子内的单个带电粒子qi
在r
处的电势，如果考虑周期性边界条件，则形式需要变为：
VSi(r)=∑nqi4πϵ0|r−ri+nL|Erfc(|r−ri+nL|√2σ)
这个形式的相互作用势相比于原始形式Vi(r)=14πϵ0∑nqi|r−ri+nL|
而言，使用了一个误差函数对实空间做了一个截断：
VSi(r)≈∑n′qi4πϵ0|r−ri+n′L|Erfc(|r−ri+n′L|√2σ),1−Erfc(|r−ri+n′L|√2σ)<ϵ
从而只需要计算有限n′
的周期性盒子即可。因为这里截断的是距离dn=|r−ri+nL|，可以用达朗贝尔判别法证明短程相互作用势在实空间的收敛性（按一般性取法先令一个δ>0）：limln→∞Erfc(ln+δ√2σ)ln(ln+δ)Erfc(ln√2σ)=limln→∞Erfc(ln+δ√2σ)Erfc(ln√2σ)=e−δ22σ2limln→∞e−lnδσ2=e−δ22σ2<1
。即：电势的短程相互作用在实空间收敛。得到短程相互作用电势的形式之后，可以进一步计算短程相互作用的电势能：
ES=∑nN−2∑i=0N−1∑j=i+1qiqj4πϵ0|rj−ri+nL|Erfc(|rj−ri+nL|√2σ)+∑|n|>0q2i4πϵ0|nL|Erfc(|nL|√2σ)
长程作用项计算
前面得到长程相互作用电势形式为：
VLi(k)=qiϵ0|k|2e−jkrie−η|k|2
同样使用短程作用项中的取值η=σ22
。做一个逆傅里叶变换变回实空间：
VLi(r)=qikxkykzϵ0∑|k|>0e−jkrie−σ2|k|22|k|2ejk⋅r
类似的可以根据达朗贝尔判别方法证明该式收敛。并且参考前面倒易空间中的k
的定义有：
ejk⋅nL=e2j|n|π=1
也就是长程相互作用项可以略去周期性盒子的求和项，因此长程相互作用电势能的形式为：
EL=12kxkykzϵ0∑|k|>0e−σ2k22k2N−1∑i=0qie−jkriN−1∑l=0qlejk⋅rl
后面这两个求和的内容形式是两个平面波函数的内积，其物理意义为把实空间的一个固定电荷按照概率幅分配到不同的倒易空间的格点上，可以定义一个倒易空间的电荷分布函数：
S(k)=N−1∑i=0qiejk⋅ri
则可以进一步简化长程相互作用势能的写法：
EL=12kxkykzϵ0∑|k|>0e−σ2k22k2|S(k)|2
需要注意的一点是，虽然这里的长程相互作用势能是收敛的，但是其中包含了点电荷i
对ri
处产生的电势能（需要去除）。而此前计算短程相互作用时可以得到的误差函数形式的长程相互作用形式：
VLi(r)=∑nqi4πϵ0|r−ri−nL|−VSi(r)=∑nqi4πϵ0|r−ri−nL|Erf(|r−ri−nL|√2σ)
虽然不收敛，但是如果在这里取一个r=ri,nL=0
，也就是前面提到的电荷自我相互作用，则长程相互作用形式变为：
VLi(ri)=qi4πϵ0√2π1σ
则可以得到在倒易空间的长程相互作用项中需要扣除的自我相互作用项为：
Eself=12N−1∑i=0qiVLi(ri)=14πϵ01√2πσN−1∑i=0q2i
总电势能
box.respectator.com）
[box.njhxpx.com）
[box.npjr.net）
[box.oilfriends.com）
[box.oydj.net）
[box.nykh6.com）
[box.qdljjx.com）
经过前面的计算，我们已经分别得到了实空间收敛的短程相互作用项、倒易空间收敛的长程相互作用项以及长程相互作用项中需要扣除的一个自我相互作用项，那么可以汇总电势能为：
E=ES+EL−Eself=∑nN−2∑i=0N−1∑j=i+1qiqj4πϵ0|rj−ri+nL|Erfc(|rj−ri+nL|√2σ)+∑|n|>0q2i4πϵ0|nL|Erfc(|nL|√2σ)+12kxkykzϵ0∑|k|>0e−σ2k22k2|S(k)|2−14πϵ01√2πσN−1∑i=0q2i
总结概要
本文介绍了Ewald求和计算方法在周期性边界条件下计算静电势能的方法。周期性的静电势函数并不是一个空间收敛的函数，通过Ewald求和可以将静电势切分为短程相互作用和长程相互作用，两项分别在实空间和倒易空间（或称傅里叶空间、k空间等）收敛。然后就可以进一步进行截断，用更少的代价获得更高精度的电势能计算结果。
