從劉維爾方程到Velocity-Verlet算法

技術背景

我們說分子動力學模擬是一個牛頓力學的過程，在使用量子化學的手段或者深度學習的方法或者傳統的力場方法，去得到某個時刻某個位置的受力之后，就可以獲取下一步的整個系統的狀態信息。這個演化的過程所使用的算法，也在歷史上演化了多次，從1967年的Verlet算法，到后來的Leap-Frog算法，再到Velocity-Verlet算法。我們可以先看一看這三種算法的形式，再從劉維爾方程出發，看看Velocity-Verlet算法的由來。

Verlet算法

我們把\(r(t+\delta t)\)和\(r(t-\delta t)\)看做是兩個函數，分別對\(r(t+\delta t)\)和\(r(t-\delta t)\)在時刻\(t\)處進行二階泰勒展開有：

\[r(t+\delta t)=r(t)+v(t)\delta t+\frac{F(t)}{2m}\delta t^2\\ r(t-\delta t)=r(t)-v(t)\delta t+\frac{F(t)}{2m}\delta t^2 \]

將上面第二個公式中的\(v(t)\delta t\)項移到左側，把\(r(t-\delta t)\)移到右側，再代入到第一個公式中，就可以得到下一步的坐標：

\[r(t+\delta t)=2r(t)-r(t-\delta t)+\frac{F(t)}{m}\delta t^2+O(\delta t^4) \]

然后再把1、2兩個公式相減，就可以得到當前時刻的速度：

\[v(t)=\frac{r(t+\delta t)-r(t-\delta t)}{2\delta t}+O(\delta t^2) \]

到這里就得到了Verlet算法的更新步驟，過程也非常的簡單。但是有個比較致命的問題是，Verlet算法的更新中，不僅僅需要到上一步的坐標位置，還需要用到上上一步的坐標位置，這就有可能導致兩個問題：

1. 第一步的更新，沒有上上一步的信息；
2. 在算法執行的過程中，每次迭代不僅僅要存儲上一步的坐標位置，還需要額外存儲上上一步的位置，更新較為麻煩，也會占據額外的空間。

目前這種傳統的Verlet算法應用已經較少，主要還是使用接下來要講到的Leap-Frog算法和Velocity-Verlet算法。

Leap-Frog算法

在蛙跳法中，引入了另外一個概念：用兩點之間的中間時刻的速度近似為兩個點之間的平均速度，這樣就可以得到一個坐標更新公式：

\[r(t+\delta t)=r(t)+v(t+\frac{\delta t}{2})\delta t \]

其中半步的速度是基於上一個半步的速度來更新的：

\[v(t+\frac{\delta t}{2})=v(t-\frac{\delta t}{2})+\frac{F(t)\delta t}{m}=v(t-\frac{\delta t}{2})-\frac{\partial V}{\partial r(t)}\cdot\frac{\delta t}{m} \]

在上面的方程中已經用勢能對坐標的偏導來替代力的計算，這也跟哈密頓力學中只有勢能項顯含坐標有關。雖然到這里我們已經完成了坐標和速度的更新，但是速度和坐標之間是不同步的，為此我們還需要用兩個半步速度取平均來計算一個時間同步的速度：

\[v(t)=\frac{v(t+\frac{\delta t}{2})-v(t-\frac{\delta t}{2})}{2} \]

由於這里只涉及到前半步的速度，而不涉及到前一步的坐標，因此Leap-Frog算法在實際應用場景中有着較為廣泛的使用。

劉維爾方程與Velocity Verlet

首先我們看一下劉維爾方程的形式：

\[\frac{d\rho(p,q,t)}{dt}=\frac{\partial\rho}{\partial t}+\hat{L}\rho=0 \]

其中劉維爾算子\(\hat{L}\)的形式為：

\[\hat{L}=\hat{L_1}+\hat{L_2}=\sum_{i=1}^{3N}\left(\frac{\partial H}{\partial p_i}\frac{\partial}{\partial q_i}-\frac{\partial H}{\partial q_i}\frac{\partial}{\partial p_i}\right) \]

其中寫成\(\hat{L_1}\)和\(\hat{L_2}\)的形式也是為了方便后面做Trotter分解：

\[\hat{L_1}=\sum_{i=1}^{3N}\frac{\partial H}{\partial p_i}\frac{\partial}{\partial q_i}\\ \hat{L_2}=-\sum_{i=1}^{3N}\frac{\partial H}{\partial q_i}\frac{\partial}{\partial p_i} \]

寫完劉維爾方程的表述之后，我們再回過頭看看劉維爾方程的物理含義，這里的密度函數\(\rho(p,q,t)\)是指在相空間中粒子的分布密度，對於整體的積分有：

\[N=\int\rho(p,q,t)dqdp \]

這里的\(N\)所表示的就是整個系統的總粒子數。因此，實際上劉維爾方程所表述的內容就是：分布函數沿着相空間的任何軌跡是常數。

Trotter-Suzuki分解

我們首先需要回顧一個知識點，雖然對於兩個常數\(a,b\)來說，其加和的指數可以等於指數的乘積，即\(e^{a+b}=e^{a}e^{b}\)，但如果是對於兩個矩陣\(A,B\)的話，類似的等式往往是不成立的。而Trotter-Suzuki公式，將其表示為一個顯式的誤差：

\[e^{\sum_{j=1}^{m}H_jt}=\prod_{j=1}^{m}e^{H_jt}+O(m^2t^2) \]

此時我們再回顧劉維爾算子的分解：\(\hat{L}=\hat{L_1}+\hat{L_2}\)，再進一步將其分解為：\(\hat{L}=\frac{\hat{L_2}}{2}+\hat{L_1}+\frac{\hat{L_2}}{2}\)，至於為什么用這個形式來分解，我也不懂，也許是嘗試出來的。基於這個形式的分解，我們將其代入到劉維爾算子的演化中。定義一個廣義參量\(x(t)=\{p(t),q(t)\}\)，則劉維爾算子對該參量的演化為：

\[\begin{align} e^{\hat{L}t}x(0)&=e^{(\frac{\hat{L_2}}{2}+\hat{L_1}+\frac{\hat{L_2}}{2})t}x(0)\\ &\approx e^{\hat{L_2}t/2}e^{\hat{L_1}t}e^{\hat{L_2}t/2}x(0)\\ &\approx e^{\hat{L_2}t/2}e^{\hat{L_1}t}\left(1+\hat{L_2}\frac{\delta t}{2}\right)x(0)\\ &=e^{\hat{L_2}t/2}e^{\hat{L_1}t}\left(x(0)-\frac{\delta t}{2}\sum_{i=1}^{3N}\frac{\partial H}{\partial q_i}\frac{\partial x(0)}{\partial p_i}\right)\\ &=e^{\hat{L_2}t/2}e^{\hat{L_1}t}x(1) \end{align} \]

觀察上述推導過程的倒數第二步，因為\(x(t)=\{p(t),q(t)\}\)，並且在相空間中所有的\(p_i\)是正交的關系，因此\(\frac{\partial x(0)}{\partial p_i}\)得到的結果全為1。又因為在哈密頓力學中有\(-\frac{\partial H}{\partial q}=\frac{dp}{dt},\frac{\partial H}{\partial p}=\frac{dq}{dt}\)。因此，假定\(x(0)=\{p_i(t_0),q_i(t_0)\},i=1,2,...,3N\)，則\(x(1)=\{p_i(t_0)+\frac{dp(t_0)}{dt}\frac{\delta t}{2},q_i(t_0)\},i=1,2,...,3N\)。使用類似的方法，我們繼續往下推導：

\[\begin{align} e^{\hat{L}t}x(0)&\approx e^{\hat{L_2}t/2}e^{\hat{L_1}t}x(1)\\ &=e^{\hat{L_2}t/2}\left(x(1)+\delta t\sum_{i=1}^{3N}\frac{\partial H}{\partial p_i}\frac{\partial x(1)}{\partial q_i}\right)\\ &=e^{\hat{L_2}t/2}x(2) \end{align} \]

其中\(x(2)=\{p_i(t_0)+\frac{dp(t_0)}{dt}\frac{\delta t}{2},q_i(t_0)+\frac{dq(t_0)}{dt}\delta t\},i=1,2,...,3N\)，同樣的方法，再完成最后一步的分解：

\[\begin{align} e^{\hat{L}t}x(0)&\approx e^{\hat{L_2}t/2}x(2)\\ &=x(2)-\frac{\delta t}{2}\sum_{i=1}^{3N}\frac{\partial H}{\partial q_i}\frac{\partial x(2)}{\partial p_i}\\ &=x(3) \end{align} \]

需要注意的是，雖然在前面\(x(0)\rightarrow x(1)\)的演化中共軛動量項在\(\hat{L_2}\)的作用下發生了變化，但是顯含的動量項保持不變，因此這里的偏導項得到的結果依然是1，那么就有\(x(3)=\{p_i(t_0)+\frac{dp(t_0)}{dt}\delta t,q_i(t_0)+\frac{dq(t_0)}{dt}\delta t\},i=1,2,...,3N\)。到這一步，問題逐漸露出端倪，我們發現在劉維爾算子的作用下，經過Trotter-Suzuki分解和Taylor展開的操作，正則坐標\(q\)和共軛動量\(p\)已經完成了一個時間單位\(\delta t\)的演化，正對應到分子動力學模擬中的單步演化。

Velocity Verlet算法

參考上一個章節中劉維爾算子的演化過程\(x(0)\rightarrow x(1)\)，我們可以先將速度進行半步演化：

\[v(t+\frac{\delta t}{2})=v(t)+\frac{F(t)}{m}\frac{\delta t}{2}+O(\delta t^3) \]

參考\(x(1)\rightarrow x(2)\)過程，我們可以將坐標進行單步演化：

\[r(t+\delta t)=r(t)+v(t+\frac{\delta t}{2})\delta t+O(\delta t^4) \]

最后參考\(x(2)\rightarrow x(3)\)的過程，將半步演化的速度再同步到單步演化：

\[v(t+\delta t)=v(t+\frac{\delta t}{2})+\frac{F(t+\delta t)}{m}\frac{\delta t}{2}+O(\delta t^2) \]

這個過程最漂亮的地方在於，演化的參數不依賴於上一步或者上半步的任何參數，只需要具備了\(v(t),r(t)\)即可演化得到\(v(t+\delta t),r(r+\delta t)\)，當然，這里面需要用量子化學或者深度學習或者是力場參數的形式，去分別計算得\(t\)和\(t+\delta 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