朗之萬動力學

在物理學中，朗之萬動力學（）是一種利用朗之萬方程對分子系統動力學進行數學建模的方法。該方法最初由法國物理學家保羅·朗之萬提出。其特點在於採用簡化模型，同時透過隨機微分方程來考量被忽略的自由度。朗之萬動力學模擬屬於蒙地卡羅模擬的一種。

概述
現實世界中的分子系統通常存在於空氣或溶劑中，而非在真空環境下孤立存在。溶劑或空氣分子的碰撞會產生摩擦，而偶發的高速碰撞則會擾動系統。朗之萬動力學試圖將分子動力學延伸，以納入這些效應。此外，朗之萬動力學允許像恆溫器那樣控制溫度，從而近似正則系綜。

朗之萬動力學模擬了溶劑的黏性特性。它並未完整建模；具體而言，該模型既未考慮靜電屏蔽效應，也未考慮疏水效應。對於密度較高的溶劑，朗之萬動力學無法捕捉其流體動力學相互作用。

對於由 N 個質量為 M 的粒子所組成的系統，其座標 X=X(t) 構成一個隨時間變化的隨機變量，由此得到的朗之萬方程為
:M\,\ddot{\mathbf{X}} = - \mathbf{\nabla} U(\mathbf{X}) - \gamma\,M\,\dot{\mathbf{X}} + \sqrt{2\,M\,\gamma\,k_{\rm B} T}\,\mathbf{R}(t)\,,
其中 U(\mathbf{X}) 為粒子相互作用勢；\nabla 為梯度算子，故 -\mathbf{\nabla} U(\mathbf{X}) 為根據粒子相互作用勢能計算出的力；點運算子為時間導數，故 \dot{\mathbf{X}} 為速度，而 \ddot{\mathbf{X}} 為加速度；\gamma 為阻尼常數（單位為時間導數），亦稱為；T 為溫度，k_{\rm B} 為波茲曼常數；而 \mathbf{R}(t) 為均值為零的δ相關平穩高斯過程，稱為高斯白雜訊，滿足
:\left\langle \mathbf{R}(t) \right\rangle = 0
:\left\langle \mathbf{R}(t)\cdot\mathbf{R}(t') \right\rangle = \delta(t - t')

此處，\delta 為狄拉克δ函數。

隨機微分方程的建模
考慮標準布朗運動或維納過程 W_t 的共變異數，我們可以得到

: \mathbb{E}(W_tW_\tau) = \min(t,\tau)

定義導數的共變異數矩陣為
:\mathbb{E}(\dot{W_t}\dot{W_\tau}) = \frac{\partial}{\partial t}\frac{\partial}{\partial \tau}\mathbb{E}(W_tW_\tau) = \frac{\partial}{\partial t}\frac{\partial}{\partial \tau}\min(t,\tau) =\delta(t-\tau)
因此，在共變異數的意義下，我們可以說

:{\rm d}W_t = \mathbf{R}(t){\rm d}t
在不失一般性的情況下，設質量 M = 1，\sigma = \sqrt{M\gamma k_{\rm B}T}，則原始的隨機微分方程將變為
: {\rm d}\dot{\mathbf{X}} = -\nabla U(\mathbf{X}){\rm d}t-\gamma {\rm d}{\mathbf{X}} +\sqrt{2}\sigma{\rm d} \mathbf{W}(t)

過阻尼的朗之萬動力學
若主要目標是控制溫度，應謹慎選用較小的阻尼常數 \gamma。隨著 \gamma 增大，系統將從慣性態延伸至擴散態。非慣性系統的朗之萬動力學極限通常被描述為。布朗動力學可視為過阻尼的朗之萬動力學，即不發生平均加速度的朗之萬動力學。在此極限下，我們有 {\rm d}\dot{X}=0，原始的隨機微分方程隨之變為

: {\rm d}{\mathbf{X}} =-\frac{1}{\gamma}\nabla U(\mathbf{X}){\rm d}t+\frac{\sqrt{2}\sigma}{\gamma}{\rm d} \mathbf{W}(t)

平移朗之萬方程可透過多種數值方法求解，這些方法在解析解的精密程度、允許的時間步長、時間可逆性（）以及零摩擦極限等方面存在差異。

朗之萬方程可推廣至分子、布朗粒子等的旋轉動力學。一種標準的（根據NIST）做法是利用基於四元數的描述來處理隨機旋轉運動。

應用
朗之萬恆溫器
朗之萬恆溫器是分子動力學中的一種恆溫器演算法，用於在指定溫度下正則系綜（NVT）。它整合了以下朗之萬運動方程：

M \ddot{\mathbf{X}} = -\nabla U(\mathbf{X}) - \gamma \dot{\mathbf{X}} + \sqrt{2\gamma k_BT} \textbf{R}(t)

-\nabla U(\mathbf{X}) 為確定性力項；\gamma 為摩擦係數，而 \gamma \dot{X} 為摩擦或阻尼項；最後一項為隨機力項（k_B：玻茲曼常數，T：溫度）。此方程式使系統能與假想的「熱浴」耦合：系統的動能透過摩擦／阻尼項消散，並從隨機力／波動中獲得；耦合強度由 \gamma 控制。此方程式可透過隨機微分方程求解器進行模擬，例如歐拉-丸山法，其中在每個積分步驟中，隨機力項被高斯隨機數取代（變異數 \sigma^2 = 2\gamma k_BT/ \Delta t，\Delta t：時間步長），或朗之萬蛙跳積分法等。此方法亦稱為朗之萬積分器。

朗之萬蒙地卡羅方法
可表示為

{\rm d} \mathbf{x}_t = - \frac{D}{k_BT} \nabla_\mathbf{x} U(\mathbf{x}_t) {\rm d}t + \sqrt{2D} {\rm d} W_t

此處，D = k_B T / \gamma 為愛因斯坦關係中的擴散係數。根據福克-普朗克方程的證明，在適當條件下，\mathbf x_t 的穩態分佈為波茲曼分布 p(\mathbf{x}) \propto e^{-U(\mathbf{x})/k_BT}。

由於 \nabla \log p(\mathbf{x})=-\nabla U(\mathbf{x})/k_BT，此方程式等同於以下形式：

{\rm d} \mathbf{x}_t = \epsilon \nabla_\mathbf{x}\log p(\mathbf x_t) {\rm d}t + \sqrt{2\epsilon} {\rm d} W_t

且 \mathbf x_t (t\to \infty) 的分布遵循 p(\mathbf{x})。換言之，由於 \nabla \log p(\mathbf{x}) 項的作用，朗之萬動力學會驅使粒子沿著梯度流朝向穩態分布 p(\mathbf{x}) 移動，同時仍允許存在某些隨機波動。這提供了一種馬可夫鏈蒙地卡羅方法，可用於從目標分佈 p(\mathbf{x}) 採樣資料 \mathbf x，此方法稱為朗之萬蒙地卡羅方法。

在許多應用中，我們會遇到一個期望的分布 p(\mathbf{x})，並希望從中採樣 \mathbf x，但直接採樣可能相當困難或效率不高。朗之萬蒙地卡羅方法提供了一種替代方案，透過採樣符合朗之萬動力學的馬可夫鏈來取得 \mathbf x \sim p(\mathbf x)，其定態即為 p(\mathbf{x})。梅特羅波利斯修正朗之萬演算法（MALA）即為一例：給定當前狀態 \mathbf x_t，MALA會利用上述朗之萬動力學提出一個新狀態 \tilde{x}_{t+1}。隨後根據梅特羅波利斯-黑斯廷斯算法決定是否接受該提案。在選擇 \tilde{x}_{t+1} 時納入朗之萬動力學，能提升運算效率，因為該動力學會驅使粒子進入 p(\mathbf{x}) 機率較高的區域，因而更可能被接受。

基於分數的生成模型
朗之萬動力學是基於分數的生成模型的基礎之一：

\mathbf{x}_{i+1} \gets \mathbf{x}_i + \epsilon \nabla_\mathbf{x} \log p(\mathbf{x}_i) + \sqrt{2\epsilon} \mathbf{z}_i,
\quad i=0,1,\cdots, K

其中 \mathbf{z}_i \sim N(0,1)。當 \epsilon \to 0 且 K \to \infty 時，生成的 \mathbf{x}_K 將收斂至目標分佈 p(\mathbf x)。基於分數的模型使用 \mathbf{s}_\theta(\mathbf{x}) \approx \nabla_\mathbf{x} \log p(\mathbf{x}) 作為近似。

與其他理論的關聯
克萊因-克拉默斯方程
作為一項隨機微分方程，朗之萬動力學方程擁有其對應的偏微分方程——，這是一項特殊的福克-普朗克方程，用以描述相空間中粒子的機率分布。原始的朗之萬動力學方程可重寫為以下一階隨機微分方程：
: {\rm d}\mathbf{X} = \mathbf{P}{\rm d} t
: {\rm d} \mathbf{P} = -\gamma\mathbf{P}{\rm d}t-\nabla U(\mathbf{X}){\rm d} t+\sqrt{2}\sigma {\rm d}\mathbf{W}(t)
現在考慮以下情況及其 (\mathbf{X},\mathbf{P}) 法則：

1.\mathbf{{\rm d}{X}} = \mathbf{P}{\rm d} t
, \mathbf{{\rm d}{P}} = -\gamma\mathbf{P}{\rm d} t-\nabla U(\mathbf{X}){\rm d} t+\sqrt{2}\sigma {\rm d}\mathbf{W}(t) 其中 (\mathbf{X}_0,\mathbf{P}_0)\sim\rho_0

1. \frac{\partial\rho}{\partial t} = -\mathbf{P}\nabla_{\mathbf{X}}\rho+\nabla_{\mathbf{P}}(\gamma\mathbf{P}\rho+\nabla_{\mathbf{X}}U(\mathbf{X})\rho)+\nabla_{\mathbf{P}}^2(\sigma_{T}^2\rho) 其中 \rho(t=0,\mathbf{X},\mathbf{P})=\rho_0

考慮一個動量與位置的一般函數

: \Psi_t = \Psi(\mathbf{X},\mathbf{P})

該函數的期望值為

: \mathbb{E}[\Psi_t] = \int \rho(t,\mathbf{X},\mathbf{P})\Psi(\mathbf{X},\mathbf{P}){\rm d} \mathbf{P}{\rm d} \mathbf{X}

對時間 t 求導，並應用伊藤公式，可得
: \mathbb{E}[\frac{\rm d}{{\rm d}t}\Psi(\mathbf{X},\mathbf{P})]
=\mathbb{E}[\nabla_{\mathbf{X}} \Psi\frac{{\rm d}\mathbf{X}}{{\rm d}t} + \nabla_{\mathbf{P}}\Psi\frac{{\rm d}\mathbf{P}}{{\rm d}t} + \sigma_T^2\nabla_{\mathbf{P}}^2\Psi\frac{1}{{\rm d} t}({\rm d} \mathbf{W}(t))^2]
此式可簡化為
: \int (\frac{\partial}{\partial t}\rho)\Psi(\mathbf{X},\mathbf{P}){\rm d} \mathbf{X}{\rm d}\mathbf{P} =\mathbb{E}[(\nabla_{\mathbf{X}} \Psi)\mathbf{P} + \nabla_{\mathbf{P}}\Psi(-\gamma\mathbf{P}-\nabla_{\mathbf{X}}U(\mathbf{X})) + \sigma_T^2\nabla_{\mathbf{P}}^2\Psi]
對右邊進行分部積分，由於無限動量或速度下密度為零，故有
: (\frac{\partial}{\partial t}\rho)\Psi(\mathbf{X},\mathbf{P}){\rm d} \mathbf{X}{\rm d}\mathbf{P} = \int(-\mathbf{P}\nabla_{\mathbf{X}}\rho+\nabla_{\mathbf{P}}(\gamma\mathbf{P}\rho+\nabla_{\mathbf{X}}U(\mathbf{X})\rho)+\nabla_{\mathbf{P}}^2(\sigma_{T}^2\rho) )\Psi(\mathbf{X},\mathbf{P}){\rm d}\mathbf{X}{\rm d}\mathbf{P}
此方程對任意 \Psi 皆成立，因此我們要求密度滿足
: \frac{\partial \rho}{\partial t} = -\mathbf{P}\nabla_{\mathbf{X}}\rho+\nabla_{\mathbf{P}}(\gamma\mathbf{P}\rho+\nabla_{\mathbf{X}}U(\mathbf{X})\rho)+\nabla_{\mathbf{P}}^2(\sigma_{T}^2\rho)
此方程稱為，是福克-普朗克方程的一種特例。它是一個偏微分方程，描述系統在相空間中機率密度的演化。

福克-普朗克方程
對於過阻尼極限，我們有 {\rm d}\mathbf{P} = 0，因此系統的演化可簡化為位置子空間。根據類似的推理，我們可以證明位置的隨機微分方程
:{\rm d} \mathbf{X} = -\frac{1}{\gamma}\nabla U(\mathbf{X}){\rm d} t +\sqrt{2}\frac{\sigma}{\gamma}\mathbf{R}(t){\rm d} t
對應於機率密度
: \frac{\partial \rho(t,\mathbf{X})}{\partial t} = \nabla_{\mathbf{X}}(\frac{1}{\gamma}\nabla_{\mathbf{X}}U(\mathbf{X})\rho(t,\mathbf{X}))+\Delta_\mathbf{X}(\frac{\sigma^2}{\gamma^2}\rho(t,\mathbf{X}))
的福克-普朗克方程。

波動-耗散定理
考慮自由粒子的朗之萬動力學（即 U(\mathbf{X})=0 ），此時動量方程將變為
: {\rm d}\mathbf{P} = -\frac{1}{\gamma} \mathbf{P}{\rm d}t +\frac{\sqrt{2}\sigma}{\gamma} {\rm d}\mathbf{W}_t
此隨機微分方程的解析解為
:\mathbf{P} = \mathbf{P}_0e^{- t/\gamma}+\frac{\sqrt{2}\sigma}{\gamma}\int_0^t {\rm e}^{-(t-t')/\gamma}{\rm d}\mathbf{W}_t'
因此，動量二階矩的平均值將變為（此處我們應用）
:\mathbb{E}(\mathbf{P}^2) = \mathbf{P}^2_0{\rm e}^{-2 t/\gamma}+ \frac{\sigma^2}{\gamma}(1-{\rm e}^{-2t/\gamma})\overset{t\to\infty}{\to}\frac{\sigma^2}{\gamma}
也就是說，當時間趨近正無窮大時的極限行為顯示，該系統的動量波動與其能量耗散（摩擦項參數 \gamma ）有關。將此結果與將粒子動能的平均值與溫度相關聯的均分定理結合
: \langle v^2\rangle = k_{\rm B}T
我們便能在如朗之萬恆溫器等應用中，確定變異數 \sigma 的數值。
:\sigma^2/\gamma = k_B T \to \sigma = \sqrt{k_{\rm B} T\gamma}
這與假設 M=1 的原始定義是一致的。

路徑積分
路徑積分表述源自量子力學。但對於朗格文隨機微分方程，我們也可以推導出相應的路徑積分。考慮以下過阻尼朗格文方程，其中為保持一般性，我們取 \gamma = \sigma = 1 ，
: {\rm d}{X} = -\nabla U({X}){\rm d}t +\sqrt{2}{\rm d}W_t
進行離散化並定義 t_n = n\Delta t ，我們得到
: {X}_{n+1} -{X}_{n} + \nabla U({X})\Delta t= \sqrt{2}(W_{t_n}-W_{t_{n-1}})\sim \mathcal{N}(0,2\sqrt{\Delta t})
因此，傳播機率將為
: P({X}_{n+1}|{X}_{n}) = \int {\rm d}\xi \frac{1}{2\sqrt{\pi\Delta t}}{\rm e}^{-\frac{\xi^2}{4\Delta t}} \delta({X}_{n+1} -{X}_{n} + \nabla U({X})\Delta t-\xi)
對狄拉克δ函數進行傅立葉變換，可得
: P = \int \frac{{\rm d}k}{2\pi}{\rm e}^{{\rm i}k({X}_{n+1}-{X}_n + \nabla U({X})\Delta t)}\int {\rm d}\xi \frac{1}{2\sqrt{\pi\Delta t}}{\rm e}^{-\frac{\xi^2}{4\Delta t}}{\rm e}^{-{\rm i}k\xi}
第二項為高斯積分，其結果為
: P = \int \frac{{\rm d}k}{2\pi}{\rm e}^{{\rm i}k({X}_{n+1}-{X}_n + \nabla U({X})\Delta t)}{\rm e}^{-k^2\Delta t}
現在考慮從初始值 X_0 到終值 X_n 的機率。
: P(\mathbf{X}_n|\mathbf{X}_0) = \int \frac{1}{2\pi}\prod_i^{N-1} {\rm d}k_i {\rm e}^{({\rm i}k_i(\dot{X} + \nabla U(X))-k_i^2)\Delta t}
取 \Delta t\to 0 的極限，可得
: P(\mathbf{X}_n|\mathbf{X}_0) = \int \mathcal{D}[k] {\rm e}^{\int_0^{t_n}({\rm i}k(\dot{X} + \nabla U(X))-k^2){\rm d} t}

參見

• 哈密顿力学
• 统计力学
• 隨機微分方程
• 朗之万方程
• 朗之萬蒙地卡羅方法

*

參考資料
外部連結

• [https://web.archive.org/web/20010619210927/http://cmm.cit.nih.gov/intro_simulation/node24.html Langevin Dynamics (LD) Simulation]