朗之万动力学

在物理学中，朗之万动力学（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