从刘维尔方程到Velocity-Verlet算法

技术背景

我们说分子动力学模拟是一个牛顿力学的过程,在使用量子化学的手段或者深度学习的方法或者传统的力场方法,去得到某个时刻某个位置的受力之后,就可以获取下一步的整个系统的状态信息。这个演化的过程所使用的算法,也在历史上演化了多次,从1967年的Verlet算法,到后来的Leap-Frog算法,再到Velocity-Verlet算法。我们可以先看一看这三种算法的形式,再从刘维尔方程出发,看看Velocity-Verlet算法的由来。

Verlet算法

我们把r(t+δt)r(tδt)看做是两个函数,分别对r(t+δt)r(tδt)在时刻t处进行二阶泰勒展开有:

r(t+δt)=r(t)+v(t)δt+F(t)2mδt2r(tδt)=r(t)v(t)δt+F(t)2mδt2

将上面第二个公式中的v(t)δt项移到左侧,把r(tδt)移到右侧,再代入到第一个公式中,就可以得到下一步的坐标:

r(t+δt)=2r(t)r(tδt)+F(t)mδt2+O(δt4)

然后再把1、2两个公式相减,就可以得到当前时刻的速度:

v(t)=r(t+δt)r(tδt)2δt+O(δt2)

到这里就得到了Verlet算法的更新步骤,过程也非常的简单。但是有个比较致命的问题是,Verlet算法的更新中,不仅仅需要到上一步的坐标位置,还需要用到上上一步的坐标位置,这就有可能导致两个问题:

  1. 第一步的更新,没有上上一步的信息;
  2. 在算法执行的过程中,每次迭代不仅仅要存储上一步的坐标位置,还需要额外存储上上一步的位置,更新较为麻烦,也会占据额外的空间。

目前这种传统的Verlet算法应用已经较少,主要还是使用接下来要讲到的Leap-Frog算法和Velocity-Verlet算法。

Leap-Frog算法

在蛙跳法中,引入了另外一个概念:用两点之间的中间时刻的速度近似为两个点之间的平均速度,这样就可以得到一个坐标更新公式:

r(t+δt)=r(t)+v(t+δt2)δt

其中半步的速度是基于上一个半步的速度来更新的:

v(t+δt2)=v(tδt2)+F(t)δtm=v(tδt2)Vr(t)δtm

在上面的方程中已经用势能对坐标的偏导来替代力的计算,这也跟哈密顿力学中只有势能项显含坐标有关。虽然到这里我们已经完成了坐标和速度的更新,但是速度和坐标之间是不同步的,为此我们还需要用两个半步速度取平均来计算一个时间同步的速度:

v(t)=v(t+δt2)v(tδt2)2

由于这里只涉及到前半步的速度,而不涉及到前一步的坐标,因此Leap-Frog算法在实际应用场景中有着较为广泛的使用。

刘维尔方程与Velocity Verlet

首先我们看一下刘维尔方程的形式:

dρ(p,q,t)dt=ρt+L^ρ=0

其中刘维尔算子L^的形式为:

L^=L1^+L2^=i=13N(HpiqiHqipi)

其中写成L1^L2^的形式也是为了方便后面做Trotter分解:

L1^=i=13NHpiqiL2^=i=13NHqipi

写完刘维尔方程的表述之后,我们再回过头看看刘维尔方程的物理含义,这里的密度函数ρ(p,q,t)是指在相空间中粒子的分布密度,对于整体的积分有:

N=ρ(p,q,t)dqdp

这里的N所表示的就是整个系统的总粒子数。因此,实际上刘维尔方程所表述的内容就是:分布函数沿着相空间的任何轨迹是常数。

Trotter-Suzuki分解

我们首先需要回顾一个知识点,虽然对于两个常数a,b来说,其加和的指数可以等于指数的乘积,即ea+b=eaeb,但如果是对于两个矩阵A,B的话,类似的等式往往是不成立的。而Trotter-Suzuki公式,将其表示为一个显式的误差:

ej=1mHjt=j=1meHjt+O(m2t2)

此时我们再回顾刘维尔算子的分解:L^=L1^+L2^,再进一步将其分解为:L^=L2^2+L1^+L2^2,至于为什么用这个形式来分解,我也不懂,也许是尝试出来的。基于这个形式的分解,我们将其代入到刘维尔算子的演化中。定义一个广义参量x(t)={p(t),q(t)},则刘维尔算子对该参量的演化为:

(1)eL^tx(0)=e(L2^2+L1^+L2^2)tx(0)(2)eL2^t/2eL1^teL2^t/2x(0)(3)eL2^t/2eL1^t(1+L2^δt2)x(0)(4)=eL2^t/2eL1^t(x(0)δt2i=13NHqix(0)pi)(5)=eL2^t/2eL1^tx(1)

观察上述推导过程的倒数第二步,因为x(t)={p(t),q(t)},并且在相空间中所有的pi是正交的关系,因此x(0)pi得到的结果全为1。又因为在哈密顿力学中有Hq=dpdt,Hp=dqdt。因此,假定x(0)={pi(t0),qi(t0)},i=1,2,...,3N,则x(1)={pi(t0)+dp(t0)dtδt2,qi(t0)},i=1,2,...,3N。使用类似的方法,我们继续往下推导:

(6)eL^tx(0)eL2^t/2eL1^tx(1)(7)=eL2^t/2(x(1)+δti=13NHpix(1)qi)(8)=eL2^t/2x(2)

其中x(2)={pi(t0)+dp(t0)dtδt2,qi(t0)+dq(t0)dtδt},i=1,2,...,3N,同样的方法,再完成最后一步的分解:

(9)eL^tx(0)eL2^t/2x(2)(10)=x(2)δt2i=13NHqix(2)pi(11)=x(3)

需要注意的是,虽然在前面x(0)x(1)的演化中共轭动量项在L2^的作用下发生了变化,但是显含的动量项保持不变,因此这里的偏导项得到的结果依然是1,那么就有x(3)={pi(t0)+dp(t0)dtδt,qi(t0)+dq(t0)dtδt},i=1,2,...,3N。到这一步,问题逐渐露出端倪,我们发现在刘维尔算子的作用下,经过Trotter-Suzuki分解和Taylor展开的操作,正则坐标q和共轭动量p已经完成了一个时间单位δt的演化,正对应到分子动力学模拟中的单步演化。

Velocity Verlet算法

参考上一个章节中刘维尔算子的演化过程x(0)x(1),我们可以先将速度进行半步演化:

v(t+δt2)=v(t)+F(t)mδt2+O(δt3)

参考x(1)x(2)过程,我们可以将坐标进行单步演化:

r(t+δt)=r(t)+v(t+δt2)δt+O(δt4)

最后参考x(2)x(3)的过程,将半步演化的速度再同步到单步演化:

v(t+δt)=v(t+δt2)+F(t+δt)mδt2+O(δt2)

这个过程最漂亮的地方在于,演化的参数不依赖于上一步或者上半步的任何参数,只需要具备了v(t),r(t)即可演化得到v(t+δt),r(r+δt),当然,这里面需要用量子化学或者深度学习或者是力场参数的形式,去分别计算得tt+δt时刻的作用力。

总结概要

本文延续历史上分子动力学模拟演化算法的发展顺序,分别讲述了Verlet、LeapFrog和Velocity-Verlet三个算法的形式,并且结合刘维尔方程,推导了Velocity-Verlet算法中的三个演化步骤的内在含义。三种不同的演化算法,都有不同的初始依赖,而对于计算过程而言,越多的初始依赖,就会涉及到越多的参数存储问题。一个好的演化算法,只需要依赖于少量的参数,而具备较高的精度。

版权声明

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

作者ID:DechinPhy

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

打赏专用链接:https://www.cnblogs.com/dechinphy/gallery/image/379634.html

腾讯云专栏同步:https://cloud.tencent.com/developer/column/91958

本文作者:Dechin的博客

本文链接:https://www.cnblogs.com/dechinphy/p/liouville.html

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

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

喜欢请打赏

扫描二维码打赏

了解更多