朗之萬動力學

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

## 概述

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

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

對於由 $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})$ 機率較高的區域，因而更可能被接受。

### 基於分數的生成模型

朗之萬動力學是基於分數的生成模型的基礎之一。根據（過阻尼）朗之萬動力學。

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

生成模型旨在產生遵循（未知資料分佈）$p(\mathbf{x})$ 的樣本。為達成此目標，基於分數的生成模型會學習一個近似分數函數 $\mathbf{s}_\theta(\mathbf{x}) \approx \nabla_\mathbf{x} \log 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 <pre> , \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$

- $\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})] <pre> =\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}$$

## 外部連結

- Langevin Dynamics (LD) Simulation