分子模擬/分子動力學

氬的分子動力學模擬

分子動力學模擬是一種模擬多體系統原子運動的工具。分子動力學用於計算經典系統的平衡和傳輸性質（粘度、熱導率、擴散、反應速率、蛋白質摺疊時間、結構和表面塗層）。要模擬分子的動力學，可以使用經典的牛頓力學。在實踐中，分子動力學模擬是透過以時間的小增量移動來執行的。它是對原子和分子的經典運動方程的求解，以獲得系統的演化。由於分子系統通常由大量粒子組成，因此不可能透過解析方法確定這些複雜系統的性質，因此 MD 被應用於這些系統，透過使用數值方法來解決這個問題。

Verlet 演算法

Verlet 演算法是一種流行的演算法，可以預測原子在未來某個時間點的位置。^[1] 該演算法使用泰勒級數逼近來近似時間 $t+\Delta t$ （從時間 t 開始的未來一個時間步長）的粒子位置。

泰勒級數基於函式在另一個點 a 的值和導數來近似點 x 處的函式。

$f(x)=f(a)+{\frac {{\textrm {d}}f(a)}{{\textrm {d}}t}}\left(x-a\right)+{\frac {1}{2!}}{\frac {{\textrm {d}}^{2}f(a)}{{\textrm {d}}t^{2}}}\left(x-a\right)^{2}+{\frac {1}{3!}}{\frac {{\textrm {d}}^{3}f(a)}{{\textrm {d}}t^{3}}}\left(x-a\right)^{3}+\ldots$

Verlet 演算法的推導

泰勒級數逼近了分子在 $t+\Delta t$ 的位置方程。

$r_{i}(t+\Delta t)\approx r_{i}(t)+{\frac {{\textrm {d}}r_{i}(t)}{dt}}(t+\Delta t-t)+{\frac {1}{2}}{\frac {{\textrm {d}}^{2}r_{i}(t)}{{\textrm {d}}t^{2}}}(t+\Delta t-t)^{2}$

$r_{i}(t+\Delta t)\approx r_{i}(t)+{\frac {{\textrm {d}}r_{i}(t)}{{\textrm {d}}t}}\Delta t+{\frac {1}{2}}{\frac {{\textrm {d}}^{2}r_{i}(t)}{{\textrm {d}}t^{2}}}\Delta t^{2}$

位置的一階和二階導數分別對應速度和加速度，加速度透過牛頓第二定律與力相關。

${\frac {{\textrm {d}}r_{i}(t)}{{\textrm {d}}t}}=v_{i}(t)$

${\frac {{\textrm {d}}^{2}r_{i}(t)}{{\textrm {d}}t^{2}}}=a_{i}(t)={\frac {F_{i}(t)}{m_{i}}}$

$r_{i}(t+\Delta t)\approx r_{i}(t)+v_{i}(t)\Delta t+{\frac {1}{2}}{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

考慮時間上向後推移一步

$r_{i}(t-\Delta t)\approx r_{i}(t)+{\frac {dr_{i}(t)}{dt}}(t-\Delta t-t)+{\frac {1}{2}}{\frac {d^{2}r_{i}(t)}{dt^{2}}}(t-\Delta t-t)^{2}$

$r_{i}(t-\Delta t)\approx r_{i}(t)-v_{i}(t)\Delta t+{\frac {1}{2}}{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

將前一步和下一步的位置相加

$r_{i}(t+\Delta t)+r_{i}(t-\Delta t)=r_{i}(t)+v_{i}(t)\Delta t+{\frac {1}{2}}{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}+r_{i}(t)-v_{i}(t)\Delta t+{\frac {1}{2}}{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

$r_{i}(t+\Delta t)+r_{i}(t-\Delta t)=2r_{i}(t)+{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

重新排列

Verlet 演算法

$r_{i}(t+\Delta t)=2r_{i}(t)-r_{i}(t-\Delta t)+{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

速度 Verlet 演算法推導

Verlet 方程的另一種形式，它涉及速度。透過使用當前時間步長的下一個時間步長的泰勒級數

$r_{i}(t+\Delta t)\approx r_{i}(t)+v_{i}(t)\Delta t+{\frac {1}{2}}{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

使用下一個時間步長的當前時間步長

$r_{i}(t)\approx r_{i}(t+\Delta t)-v_{i}(t+\Delta t)\Delta t+{\frac {1}{2}}{\frac {F_{i}(t+\Delta t)}{m_{i}}}\Delta t^{2}$

透過替換

速度 Verlet 演算法

$v_{i}(t+\Delta t)=v_{i}(t)+{\frac {\Delta t}{2m_{i}}}\left[F_{i}(t)+F_{i}(t+\Delta t)\right]$

另請參閱

維基百科：Verlet 積分 Chemwiki: Lennard-Jones

時間步長

大多數分子動力學模擬的時間步長在飛秒量級，這是化學鍵振動的量級。分子動力學模擬受最高頻率振動的限制，通常時間步長應低於最高頻率的十分之一。通常，有機分子分子動力學模擬使用 1-2 fs ( ${10^{-15}}$ s) 的時間步長，因為含有氫的鍵的振動速度很快。

${\text{Number of time steps}}={\frac {\text{length of simulation}}{\text{time step length}}}$

第一個問題是，一些過程，如蛋白質摺疊，需要相對較長的模擬時間才能完整地觀察到。為了在 MD 模擬中觀察這些過程，我們需要至少微秒級的模擬。

${\text{Number of time steps}}={\frac {10^{-6}}{10^{-15}}}=1{\text{ billion time steps}}$

即使對於計算機來說，這對於時間步長來說也是一個非常大的數字。為了解決這個問題，應用了一些約束條件。許多分子動力學模擬將所有含氫鍵約束為平衡值。因為氫是一種非常輕的原子，並且具有非常快的鍵振動。透過將這些含氫鍵約束為固定在平衡值，快速振動將被忽略。

另一個問題是共價鍵振動的振動能級在量子能級之間有很大的空間，因此使用經典力學對其描述很差，最好假設鍵是剛性的。

時間可逆性

一個透過牛頓運動定律傳播的孤立系統應該是時間可逆的。也就是說，如果你。在許多情況下，分子動力學演算法是時間可逆的。這意味著透過向前執行軌跡，然後向後執行相同數量的步長，可以按相同的步長序列返回起點。因此，兩個具有相同初始條件和時間步長的 MD 模擬將遵循完全相同的軌跡。

第二個結果是該系統被稱為**確定性**。如果你知道粒子在給定時間點的位移和速度，並且可以計算作用在它們身上的任何位置的力，原則上，你可以積分運動定律以確定系統在任何未來時間點的狀態。透過這種方式，我們可以根據我們現在所處的狀態來*確定*系統的未來狀態（即位置和速度）。

使用分子動力學取樣熱力學系綜

微正則系綜 (NVE)

使用近似牛頓動力學的積分器進行的分子動力學模擬（例如 Verlet 演算法）將守恆系統的總能量（即動能和勢能之和在每個時間步長保持恆定）。因此，這些模擬將對微正則系綜 (NVE) 進行取樣，因為系統的總能量 (E) 是守恆的。

Verlet 方程以及所有其他積分器都與泰勒級數展開中高階項的截斷相關聯的一些誤差。當使用更大的時間步長時，此誤差會更顯著。此外，由於在有限精度計算機（例如 32 或 64 位暫存器）上執行這些計算而導致的舍入誤差會產生數值誤差。因此，分子動力學模擬的總能量會在模擬過程中偏離其初始值。

取樣不同的系綜

要取樣不同的系綜，必須修改運動方程。這通常被描述為使用恆溫器或恆壓器進行分子動力學模擬。

正則系綜 (NVT)

在正則系綜中，物質的量 (N)、體積 (V) 和溫度 (T) 是恆定的。在這個系綜中，吸熱和放熱過程的能量與恆溫器交換。已經開發了幾種方法來修改分子動力學模擬的運動方程，以便它們對正則系綜進行取樣，例如 wikipedia::Andersen 恆溫器和 wikipedia::Nosé–Hoover 恆溫器。