Typesetting math: 100%

Ewald求和在分子静电势能计算中的应用

技术背景

分子动力学模拟中,计算周期性边界条件的静电势常被视作计算的瓶颈之一。形式上是比较容易的,例如不考虑周期性边界条件的话,静电势能就是:

E=14πϵ0i=0N2j=i+1N1qiqjrij

如果考虑周期性边界条件,那么静电势能变为:

E=14πϵ0ni=0N2j=i+1N1qiqj|rij+nL|+14πϵ0|n|>0qi2|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|rri|

如果考虑上周期性边界条件就是:

Vi(r)=14πϵ0nqi|rri+nL|

因为这个无穷多的求和没办法直接计算,只能明确电势具有周期性:Vi(r)=Vi(r+nL)

静电场泊松方程

空间中的电荷ir处的电势,由泊松方程给出:

ΔVi(r)=ρi(r)ϵ0

其中Δ是拉普拉斯算子,ρi表示电荷密度。如果考虑一个点电荷的的情况,那么就有:ρi(r)=qiδ(rri)。进而写出欧几里得空间中的泊松方程:

2Vi(r)=qiδ(rri)ϵ0

其中2Vi(r)=2Vi(r)x2+2Vi(r)y2+2Vi(r)z2。再考虑周期性边界条件,每个盒子中都有一个点电荷qi,于是方程应该写为:

2Vi(r)=nqiδ(rri+nL)ϵ0

傅里叶空间泊松方程

对上述单点形式的泊松方程的两边同时进行傅里叶变换(关于傅里叶变换的理解,可以参考前序文章1前序文章2,有比较详细的原理介绍和相关代码实现)有:

(2Vi(r)x2+2Vi(r)y2+2Vi(r)z2)e2πjkrdr=qiϵ0δ(rri)e2πjkrdr

先使用分部积分计算左边中的一项:

2Vi(r)x2e2πjkrdr=Vi(r)xe2πjkr|+2πjkxVi(r)xe2πjkrdr

需要注意的是,这里Vi(r)x=Fi(r)的物理意义是作用力,那么我们就可以取第一类边界条件:

limxVi(x)=0

这样根据微分的定义有:

limxVi(x)x=limxlimϵ0Vi(x+ϵ)Vi(x)ϵ=0

上面这个极限代入了Vi(r)的单点形式,其意义为,位于ri处的点电荷qi,在无穷远处生成的电场为0,对无穷远处的电荷q的作用力也是0,这里不考虑周期性边界条件。则有:

2Vi(r)x2e2πjkrdr=2πjkxVi(r)xe2πjkrdr=2πjkx[Vi(r)e2πjkr|x=x=+2πjkxVi(r)e2πjkrdr]=4π2kx2Vi(k)

同理可得泊松方程左侧形式为:

2Vi(r)e2πjkrdr=4π2k2Vi(k)

而右侧形式需要用到狄拉克函数的抽样特性:

δ(t)f(t)dt=δ(t)f(0)dt=f(0)δ(t)dt=f(0)

即:

δ(rri)e2πjkrdr=e2πjkri

得到傅里叶空间的单点泊松方程:

4π2k2Vi(k)=qiϵ0e2πjkri

倒易空间

涉及到傅里叶空间,我们很自然的想到使用固体物理学的倒易空间变换,也就是把周期性盒子当作一个原胞。根据倒易空间(也叫k空间)晶格矢(倒格矢)定义有:

k=m1k1+m2k2+m3k3

其中:

kiLj=δi,jδi,j={1,i=j0,ij

按照我们常用的长方体周期性边界条件:

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πLx2L1

其中Ω表示周期性盒子的体积,类似的有:

k2=2π(L3×L1)L2(L3×L1)=2πLy2L2

以及

k3=2π(L1×L2)L3(L1×L2)=2πLz2L3

经过倒易空间变换之后,原胞体积从Ω=LxLyLz变成:Ω=(2π)3LxLyLz。因为在前一步傅里叶空间的泊松方程中我们注意到k前面总是带了一个2π,这里不妨使用倒易晶格矢的定义对k的形式做一个简化:

k=(2πLx,2πLy,2πLz)

这样一来傅里叶空间的泊松方程可以简写为:

Vi(k)=qiϵ0k2ejkri

其中k2=|k|2=4π2m12Lx2+4π2m22Ly2+4π2m32Lz2,可以实现实空间到倒易空间的变换。

衰减函数构造

对于上述傅里叶变换之后的单点电势的形式,即使我们对整个k空间进行积分,也是一个发散的结果。所以这里用到了一个非常特别的思想,由Edwald提出,把静电能量项分为远程相互作用项和短程相互作用项,分别在倒易空间和实空间收敛,这样就可以精确计算静电能。实际操作的时候有不同的推导过程,我们这里引用一种比较“数学”的推导方法(参考链接1)。

首先我们构造一个衰减函数ek2t,这个衰减函数有个特性是:

0ek2tdt=(1k2ek2t)|0=1k2

这样我们就可以用这个积分形式替换掉傅里叶-泊松方程中的1k2项:

Vi(k)=qiϵ0ejkri0ek2tdt

因为这里使用的是从0到无穷大的一个积分形式,那么我们就可以实现一个截断,将其划分成两个积分的加和,假如我们在η处做一个截断,则有:

Vi(k)=qiϵ0ejkri(0ηek2tdt+ηek2tdt)

这里取短程(Short Term)相互作用为:

ViS(k)=qiϵ0ejkri0ηek2tdt

以及长程(Long Term)相互作用为:

ViL(k)=qiϵ0ejkriηek2tdt=qiϵ0k2ejkrieηk2

短程作用项计算

按照预期,划分了短程作用项和长程作用项之后,应该可以得到一个实空间收敛的短程相互作用,我们对短程作用做一个逆傅里叶变换:

ViS(r)=1kxkykzkViS(k)eikr=qikxkykzϵ0kejk(rri)0ηek2tdt=qikxkykzϵ00ηkejk(rri)ek2tdt

很明显,积分内的求和项是一个指数平方函数的离散傅里叶变换。而我们可以知道,正态分布函数f(ξ)=12πσeξ22σ2的傅里叶变换和逆傅里叶变换不改变其分布形式:

F(k)=f(ξ)ejkξdξ=12πσeξ22jkξσ22σ2dξ=12πσeξ22jkξσ2(jkσ2)2+(jkσ2)22σ2dξ=12πσe(jkσ2)22σ2e(ξ+jkσ2)22σ2dξ

由于积分项只是一个实空间的积分,其本质还是一个正态分布函数的积分,我们知道其积分结果是一个常数,所以有:

F(k)=ek2σ22

也是一个正态分布,只是其均方差从σ变成了1σ,也就是其积分结果为:

F(k)dk=2πσ

同样的道理我们也可以计算得,正态分布函数得逆傅里叶变换结果也依然是一个正态分布函数,其均方差也是1σ
那么回到短程静电势,先做个变量替换t=σ22,σ0,则有dt=σdσ,σt=η=2η,σt=0=0。此时短程相互作用势为:

ViS(r)=qikxkykzϵ002η(kejk(rri)ek2σ22)σdσ=qiϵ002ηe(rri)22σ2Ncoeσdσ

这里Ncoe是一个用于归一化正态分布的常数:

Ncoe=re|rri|2σ2dr=e(zzi)22σ2[e(yyi)22σ2(e(xxi)22σ2dx)dy]dz=(2πσ2)32

所以有:

ViS(r)=qiϵ002ηe(rri)22σ2(2πσ2)32σdσ=qi(2π)32ϵ002ηe(rri)22σ2σ2dσ

这里用一个变量替换:y=|rri|2σ,则有:σ=|rri|2y,其微分变换形式为:dσ=|rri|2y2dy,其边界为:yσ=2η=|rri|2η,yσ=0=,代入得:

ViS(r)=qi(2π)32ϵ0|rri|2η2ey2|rri|dy=qi2π32ϵ0|rri||rri|2ηey2dy=qi2π32ϵ0|rri||rri|2η+ey2dy

此时使用一个误差函数Erfc(y)=2πyex2dx代入进行替换:

ViS(r)=qi4πϵ0|rri|Erfc(|rri|2η)

因为η是我们手动引入的一个常数参量,如果考虑η=σ22,那么形式就变成了:

ViS(r)=qi4πϵ0|rri|Erfc(|rri|2σ)

这个形式的短程相互作用势,表示的是单个盒子内的单个带电粒子qir处的电势,如果考虑周期性边界条件,则形式需要变为:

ViS(r)=nqi4πϵ0|rri+nL|Erfc(|rri+nL|2σ)

这个形式的相互作用势相比于原始形式Vi(r)=14πϵ0nqi|rri+nL|而言,使用了一个误差函数对实空间做了一个截断:

ViS(r)nqi4πϵ0|rri+nL|Erfc(|rri+nL|2σ),1Erfc(|rri+nL|2σ)<ϵ

从而只需要计算有限n的周期性盒子即可。因为这里截断的是距离dn=|rri+nL|,可以用达朗贝尔判别法证明短程相互作用势在实空间的收敛性(按一般性取法先令一个δ>0):limlnErfc(ln+δ2σ)ln(ln+δ)Erfc(ln2σ)=limlnErfc(ln+δ2σ)Erfc(ln2σ)=eδ22σ2limlnelnδσ2=eδ22σ2<1。即:电势的短程相互作用在实空间收敛。得到短程相互作用电势的形式之后,可以进一步计算短程相互作用的电势能:

ES=ni=0N2j=i+1N1qiqj4πϵ0|rjri+nL|Erfc(|rjri+nL|2σ)+|n|>0qi24πϵ0|nL|Erfc(|nL|2σ)

长程作用项计算

前面得到长程相互作用电势形式为:

ViL(k)=qiϵ0|k|2ejkrieη|k|2

同样使用短程作用项中的取值η=σ22。做一个逆傅里叶变换变回实空间:

ViL(r)=qikxkykzϵ0|k|>0ejkrieσ2|k|22|k|2ejkr

类似的可以根据达朗贝尔判别方法证明该式收敛。并且参考前面倒易空间中的k的定义有:

ejknL=e2j|n|π=1

也就是长程相互作用项可以略去周期性盒子的求和项,因此长程相互作用电势能的形式为:

EL=12kxkykzϵ0|k|>0eσ2k22k2i=0N1qiejkril=0N1qlejkrl

后面这两个求和的内容形式是两个平面波函数的内积,其物理意义为把实空间的一个固定电荷按照概率幅分配到不同的倒易空间的格点上,可以定义一个倒易空间的电荷分布函数:

S(k)=i=0N1qiejkri

则可以进一步简化长程相互作用势能的写法:

EL=12kxkykzϵ0|k|>0eσ2k22k2|S(k)|2

需要注意的一点是,虽然这里的长程相互作用势能是收敛的,但是其中包含了点电荷iri处产生的电势能(需要去除)。而此前计算短程相互作用时可以得到的误差函数形式的长程相互作用形式:

ViL(r)=nqi4πϵ0|rrinL|ViS(r)=nqi4πϵ0|rrinL|Erf(|rrinL|2σ)

虽然不收敛,但是如果在这里取一个r=ri,nL=0,也就是前面提到的电荷自我相互作用,则长程相互作用形式变为:

ViL(ri)=qi4πϵ02π1σ

则可以得到在倒易空间的长程相互作用项中需要扣除的自我相互作用项为:

Eself=12i=0N1qiViL(ri)=14πϵ012πσi=0N1qi2

总电势能

经过前面的计算,我们已经分别得到了实空间收敛的短程相互作用项、倒易空间收敛的长程相互作用项以及长程相互作用项中需要扣除的一个自我相互作用项,那么可以汇总电势能为:

E=ES+ELEself=ni=0N2j=i+1N1qiqj4πϵ0|rjri+nL|Erfc(|rjri+nL|2σ)+|n|>0qi24πϵ0|nL|Erfc(|nL|2σ)+12kxkykzϵ0|k|>0eσ2k22k2|S(k)|214πϵ012πσi=0N1qi2

总结概要

本文介绍了Ewald求和计算方法在周期性边界条件下计算静电势能的方法。周期性的静电势函数并不是一个空间收敛的函数,通过Ewald求和可以将静电势切分为短程相互作用和长程相互作用,两项分别在实空间和倒易空间(或称傅里叶空间、k空间等)收敛。然后就可以进一步进行截断,用更少的代价获得更高精度的电势能计算结果。

版权声明

本文首发链接为:https://www.cnblogs.com/dechinphy/p/ewald.html

作者ID:DechinPhy

更多原著文章:https://www.cnblogs.com/dechinphy/

请博主喝咖啡:https://www.cnblogs.com/dechinphy/gallery/image/379634.html

参考链接

  1. http://micro.stanford.edu/mediawiki/images/4/46/Ewald_notes.pdf

本文作者:DECHIN

本文链接:https://www.cnblogs.com/dechinphy/p/18434945/ewald

版权声明:本作品采用CC BY-NC-SA 4.0许可协议进行许可。

posted @   DECHIN  阅读(430)  评论(0编辑  收藏  举报
点击右上角即可分享
微信分享提示
评论
收藏
关注
推荐
深色
回顶
收起

喜欢请打赏

扫描二维码打赏

了解更多