Workshop - 01 - 分子动力学模拟简介#

可以参考这篇综述：

如何定义原子间的相互作用#

在分子动力学中，每个分子都被近似为一个由经典质点或相互作用位点构成的系统，通过求解该系统的经典牛顿力学方程，我们可以得到各粒子的运动规律（蛋白质各分子的运动也符合牛顿经典力学的规律）

模型简化#

• 玻恩-奥本海默近似：原子核的质量比电子大数千倍，电子运动极快，能瞬间适应原子核位置的改变。
• 物理简化：力场把“原子核 + 电子云”打包成一个经典刚性质点，赋予其固定的质量 $m$ 和等效净电荷 $q$。
• 分子间的化学键不再是波函数的重叠，而被简化为符合经典物理规律的弹簧（简谐振子），原子间非共价碰撞则被看作带有引力和斥力的刚球。

体系的总势能 $V_{\mathrm{tot}}(\mathbf{r}_1, \dots, \mathbf{r}_N)$ 是一个输入所有原子坐标、输出一个标量能量值的函数。若体系有 $N$ 个原子，每个原子有 $(x, y, z)$ 三个自由度，整个系统就处于一个 $3N$ 维的能量超曲面（Energy Landscape） 上：

• 山峰（高势能区）：两个原子距离过近产生剧烈电子斥力，或者化学键被极度拉伸。
• 山谷（低势能区）：化学键处于平衡长度、带相反电荷的原子处于稳定吸引距离。

$$\nabla_i V = \left( \frac{\partial V}{\partial x_i}, \mathbf{\hat{x}} + \frac{\partial V}{\partial y_i} \mathbf{\hat{y}} + \frac{\partial V}{\partial z_i} \mathbf{\hat{z}} \right)$$

梯度的数学意义：函数在当前点上升最快的方向（最陡峭的上坡方向）。

物理意义与负号：自然界的所有物理系统都有自发退回能量最低点的趋势。小球在斜坡上绝不会自己往上爬，只会顺着最陡的方向滚向谷底。

负梯度 $-\nabla_i V$：代表当前位置下降最快的方向（最陡峭的下坡方向）。这个坡度越陡，原子被推向低能态的驱动力（即受力 $\mathbf{F}_i$）就越大。

轨迹的计算#

确定了力 $\mathbf{F}_i$ 后，牛顿第二定律给出加速度：

$$\mathbf{a}_i(t) = \frac{\mathbf{F}_i(t)}{m_i}$$

在宏观天文学中，二体问题可以通过微积分算出解析轨道公式；但在包含数万个原子的微观体系中，原子间受力错综复杂，不存在解析解，计算机只能依靠数值时间积分（Numerical Time Integration）：

1. 时间切片：选定一个极短的时间步长 $\Delta t$（通常为 $1 \sim 2\ \text{fs}$，即 $10^{-15}\ \text{s}$）。在这样短暂的瞬间，假定加速度几乎恒定。
2. 位置与速度推演（以基础运动学为例）：

$$\mathbf{r}_i(t + \Delta t) \approx \mathbf{r}_i(t) + \mathbf{v}_i(t)\Delta t + \frac{1}{2}\mathbf{a}_i(t)\Delta t^2$$ $$\mathbf{v}_i(t + \Delta t) \approx \mathbf{v}_i(t) + \mathbf{a}_i(t)\Delta t$$

3. 循环迭代：推演到新坐标后，原子间距离改变 $\to$ 重新计算势能与力 $\to$ 得到新一帧的加速度 $\to$ 再推演下一帧。循环数百万到数十亿次，就拼合成了分子随时间运动的连续轨迹

在这里，我们定义体系中的总势能为 $V_{tot}$，即模拟系统的总势能由键合相与非键合相组成：

$$V_{\mathrm{tot}} = V_{\mathrm{bonded}} + V_{\mathrm{non\text{-}bonded}} \tag{1}$$

其中， $V_{bonded}$ 是与化学键、键角以及由分别包含2、3和4个粒子的基团之间的二面角相关的势能之和，也就是成键势能

$$V_{\mathrm{bonded}}=V_{\mathrm{bond}}+V_{\mathrm{angle}}+V_{\mathrm{dihedral}}+V_{\mathrm{improper\ dihedral}}\tag{2}$$ $$\boxed{\text{总成键势能} = \text{键伸缩} + \text{键角弯曲} + \text{二面角扭转} + \text{非正常二面角}}$$

其中， $V_{\mathrm{non\text{-}bonded}}$ 反映了任意一对粒子 i,j 之间的相互作用，非键相互作用能是范德华相互作用与库仑相互作用之和，其中范德华相互作用通常采用 Lennard-Jones (L-J) 势来建模（ $V_{L-J}$）：

$$V_{\text{non-bonded}} = V_{\text{van der Waals}} + V_{\text{Coulomb}}$$ $$V_{\text{L--J}}(r) = 4\varepsilon_{ij} \left[ \left(\frac{\sigma_{ij}}{r}\right)^{12} - \left(\frac{\sigma_{ij}}{r}\right)^{6} \right]$$ $$V_{\text{coulomb}}(r) = \frac{q_i q_j}{4\pi\varepsilon_0 \varepsilon_r r}$$

力场的区分#

• 力场（FF） = 函数形式（Functional forms） + 参数集（Parameters）。它包含了计算体系成键和非键势能所需的全部信息。
• 根据描述分子时的空间分辨率（Spatial resolution），力场被划分为全原子力场（All-atom）、联合原子力场（United-atom）与粗粒化力场（Coarse-grained）。

空间分辨率划分#

全原子力场（All-atom, AA）#

• 做法：体系中的每一个原子（包括所有的轻氢原子）都被独立看作一个力学质点求解
• 适用与局限：精度最高，能捕捉微观细节，但计算极其昂贵，通常只能模拟几万到几十万个原子、微秒（$\mu\text{s}$）时间尺度的中小型体系

联合原子力场（United-atom, UA）#

• 做法：保留所有非氢“重原子”，但将非极性基团上的氢原子（如 $-\mathrm{CH}_3$、$-\mathrm{CH}_2-$）隐式合并到与其相连的碳原子中，把整个基团视为单个复合质点（如 GROMOS 力场）。
• 为什么这样做？

1. 显著减少系统粒子总数；
2. 消除分子中频率最高的 C-H 键伸缩振动，从而允许使用更大的积分步长，大幅降低采样计算成本。

粗粒化力场（Coarse-grained, CG）#

• 做法：为了在更长时间尺度上模拟更大规模的生物复合物（如光合系统 PSII 超大复合体、大面积生物膜），进一步做“镜头拉远（Zoom out）”：将 3～6 个非氢重原子直接打包融合成一个“粗粒化珠子”（CG bead）。
• 代表力场（Martini FF）：当前应用最广泛的粗粒化力场，大约将 4 个重原子归为一个物理化学珠子（例如把一个酯基整体当作一个交互位点），参数主要根据宏观实验的油水分配系数（Partitioning data）进行标定。

动力学模拟的工作流#

构建模拟盒子 (Building the Simulation Box)#

1. 结构清洗与修补
2. 质子化状态判定与加氢（生成拓扑）：蛋白质在生理 pH（7.4）下，酸碱氨基酸带电状态不同。最典型的是组氨酸（His），需要判定它是 $\delta$ 位加氢（HID）、$\varepsilon$ 位加氢（HIE），还是两处都加氢带正电（HIP），以及二硫键，根据力场补全所有氢原子

能量最小化 (Energy Minimization)#

• 在初始搭建拼装时，原子间难免存在微小的几何重叠（clash）。
• 若直接进行动力学积分，短程排斥力（L-J 势的 $r^{-12}$ 项）会产生极其巨大的排斥力，导致系统速度发散爆炸（blow-up）。
• 通常采用最速下降法（Steepest Descent）或共轭梯度法优化原子空间坐标，让体系滑落到局部的低能量平缓状态。

分步弛豫（Stepwise Relaxation）#

• 加水算法的妥协：算法会往体系里面加水，为了防止水原子与蛋白质重叠引发极大的范德华排斥力，只要是水分子和蛋白原子重叠，就会被删掉
• 空腔与表面断裂：这会导致蛋白质深层的活性口袋、离子通道的孔道内部、以及蛋白质凹凸不平的表面缝隙中，留下许多原本在生理状态下应该充满水、但此时却是绝对真空的空洞
• 无保护下的灾难：
  • 塌陷变形，蛋白质原有结构发生破坏
  • 剪切冲击引发发散，加速度非常巨大，体系被破坏

简谐位置限制势#

在原子的初始实验坐标 $\mathbf{r}_i^0$ 处，虚拟地挂上一根刚性极强的弹簧：

$$V_{\mathrm{posres}}(\mathbf{r}_i) = \frac{1}{2} k_{\mathrm{pr}} \Vert{}\mathbf{r}_i - \mathbf{r}_i^0\Vert{}^2$$

对坐标求导后，原子受到的额外恢复力为：

$$\mathbf{F}_{\mathrm{posres}} = -k_{\mathrm{pr}} (\mathbf{r}_i - \mathbf{r}_i^0)$$

• 效果：水分子、离子和脂质分子是完全自由的（$k_{\mathrm{pr}} = 0$）；而蛋白质的骨架原子甚至侧链重原子则无法发生比较大的位移
• 原子可以在实验坐标周围发生极微小的弹性晃动，但绝对无法发生大尺度的位移或构象翻转

设定边界条件与物理系综 (Application of Boundary Conditions)#

周期性边界条件#

用于模拟无限宏观连续体系（如溶液中的蛋白质折叠、水溶液化学反应、金属晶体内部），目的是彻底消除器壁效应和表面张力

类似于吃豆人，从左边的盒子出去，然后从右边的盒子出来

速度与动量完全连续：

$$\mathbf{r}_{\text{new}} = \mathbf{r}_{\text{old}} - \mathbf{L}, \quad \mathbf{v}_{\text{new}} = \mathbf{v}_{\text{old}}$$

反射边界条件（Reflective Boundary Conditions）#

• 模拟盒子四周被赋予了不可逾越的物理外壁，空间在边界处被硬性截断
• 动量发生突变：当原子碰撞垂直于 $x$ 轴的壁面时，法向动量瞬间反向：

$$v_x \to -v_x, \quad v_y \to v_y, \quad v_z \to v_z$$

或者通过边界排斥势（如 Lennard-Jones 9-3 壁面势）施加一个向内的推力将其推回

体系中压强的计算方式：维里定理（Virial Theorem）#

$$P = \underbrace{\frac{N k_B T}{V}}_{\text{动能项（由温度决定）}} + \underbrace{\frac{1}{3V}\sum_{i<j} \mathbf{r}_{ij} \cdot \mathbf{F}_{ij}}_{\text{维里项（由原子间相互作用力决定）}}$$

• 动能项：原子跑得越快（温度 $T$ 越高），撞击产生的冲量压强就越大
• 维里项：原子之间是互相吸引还是互相排斥，产生的内应力有多大

$NVT$预平衡（恒温恒容 / 升温相）#

• 控制的物理量：粒子数 $N$（固定）、体积 $V$（固定）、温度 $T$（受控）
• 核心任务：给冰冷的原子加热，建立目标温度下的速度分布
  • 能量最小化是在 $0\ \text{K}$ 静态下完成的，所有原子的初速度全部为 0
  • 在 $NVT$阶段，开启控温器（Thermostat），计算机根据麦克斯韦-玻尔兹曼分布随机赋予每个原子一个初始热运动速度，将系统平缓加热至目标实验温度（例如 $300\ \text{K}$ 或人体体温 $310\ \text{K}$）
• 此时的状态：模拟盒子的尺寸是死死锁定的（体积 $V$ 恒定），同时对蛋白质施加位置限制（Position Restraints），仅让水分子在特定空间内受热自由运动并建立水化层

$NPT$预平衡（恒温恒压 / 密度平衡相）#

• 控制的物理量：粒子数 $N$（固定）、压强 $P$（受控）、温度 $T$（受控）
• 核心任务：让盒子自由形变，调整体系至真实的物理密度与生理压强
  • 在构建水盒子时，人工填充的水密度和真实生理环境往往存在细微偏差（可能偏稀或偏密）
  • 在这一步，保持控温器的同时，进一步开启控压器（Barostat），将压强设定在标准大气压（$1\ \text{bar}$）
  • 控压器会根据体系内部的微观维里应力，自动微调模拟盒子的大小与体积 $V$，像活塞一样缓慢压缩或微扩体系，直到溶液密度精准稳定在生理密度（如常温纯水约为 $1.0\ \text{g/cm}^3$）

动力学模拟的预处理部分总结#

标准分子动力学采样存在的局限性与非标准动力学模拟采样方法#

标准分子动力学（Standard MD）存在的核心问题#

• 时间尺度瓶颈与各态历经性（Ergodicity）缺失
  • 全原子标准 MD 需显式考虑轻原子（如氢原子）的高频振动，数值积分步长 $\Delta t$ 必须限制在飞秒量级（约几飞秒），使得模拟长时间轨迹（大于微秒）的计算成本极高
  • 大多数标准 MD 模拟的时长远不足以从平衡分布中充分采样，导致各态历经假说失效（即模拟时间内的物理量平均值不等于真实的系综平均值）
• 高能垒陷阱（Energy Barrier Trapping）
  • 标准 MD 容易长时间困在势能面的局部极小值中，很难自发跨越显著的高自由能垒
  • 许多关键的生物物理过程（如蛋白质构象转换、配体结合与解离、色素猝灭态转变）涉及高能中间态或过渡路径，在标准 MD 的可及时间内极少被访问到
• 固定质子化状态无法反映真实 pH 响应
  • 标准 MD 在模拟开始时必须为每个可电离残基手动指定固定的质子化状态，并在整个轨迹中保持不变
  • 这种处理将采样严格限制在单一质子化状态对应的构象子空间内，切断了质子化状态与空间构象之间的动态耦合，无法模拟生理环境下的连续 pH 响应（例如类囊体腔内酸化诱导的构象开关）
• 经典力场缺失电子自由度
  • 标准 MD 基于玻恩-奥本海默近似，将原子核作为经典质点在电子基态势能面上推演，完全忽略了电子自由度
  • 因此，它无法直接处理光合作用中的光激发、激发能量传递、电荷分离以及化学键断裂/形成（如水裂解 Kok 循环）等量子力学效应
• 空间与时间尺度的双重鸿沟
  • 面对类囊体膜上数百纳米尺度的超分子复合体组装、蛋白质侧向重排（如状态转换 State transitions）及脂质-蛋白长程扩散等宏观慢过程，标准全原子 MD 的算力完全无法支撑

增强采样（Enhanced Sampling）#

自由能景观#

热力学状态（Thermodynamic State）#

在常温溶液中，蛋白质等生物大分子并不是僵硬不动的晶体，而是在不断振动和摇摆。一个“状态”（比如“结合态”与“未结合态”，或“活性构象”与“非活性构象”）并不是单一的具体坐标，而是一群相似微观构象构成的集合（系综）。

自由能景观（Free Energy Landscape / Surface）#

这是一张多维的“地势图”。图上的“山谷”（能量极小值区域）对应着分子最容易停留的稳定或亚稳状态；图上的“山脊/垭口”（高能垒区域）对应着构象转变的过渡态。

热力学状态差异（$\Delta G$）：#

即两个山谷底部之间的相对海拔差（自由能差 $\Delta G$）

这个海拔差决定了两个状态在平衡时的数量分布比例（平衡常数 $K$）：

$$\Delta G = -k_B T \ln K$$

• 山谷与山谷之间的“山峰高度”（活化能垒 $\Delta G^\ddagger$）则决定了构象翻转或反应进行的快慢速率

为什么常规采样算不准#

常规模拟的统计缺陷主要在于高能垒处的“采样黑洞”：

指数级稀少的样本数：#

假设两个构象状态之间隔着一道 $15\ k_BT$（约 $9\text{ kcal/mol}$）的常见能垒。根据玻尔兹曼分布，体系处于山顶过渡态的概率只有底部的：

$$e^{-15} \approx 3 \times 10^{-7}$$

在常规纳秒或微秒级 MD 模拟中，哪怕走了数千万步，分子出现在山顶附近的次数很可能也是 0 次。

$\ln(0)$ 导致的统计崩溃：#

在计算公式 $F(x) \propto -k_B T \ln P(x)$ 中，如果山峰处的采样计数为 0（即 $x=0, lnP(0)$），对数没有定义；即便偶尔跑过去 1 次（即 $x=1, lnP(1)=0$），也会伴随极其巨大的统计误差。其结果是：地势图的山谷部分测得较准，但山谷与山谷之间的山峰和整个路径全是断裂的空白，导致你根本无法准确得出两座山谷之间的真实相对高度差（$\Delta G$）

为什么要做增强采样#

1. 突破时间尺度瓶颈与各态历经性限制（Ergodic Limitation） 常规全原子 MD 的时间步长通常在飞秒（fs）量级，模拟时长通常受限于纳秒（ns）到微秒（µs）级别。在这个时间尺度内，模拟体系通常无法在平衡分布中充分采样，导致不满足各态历经假说（即体系在有限时间内的长时间平均无法等同于真实的系综平均）
2. 跨越能垒，避免深陷局部能量极小值（Energy Trapping） 实际的生物分子过程（例如光合复合体中的构象转换、色素-蛋白相互作用重组、辅因子/配体的结合与解离）往往涉及由高能垒分隔的多个不同自由能极小点。在常规 MD 中，受限于室温下的热运动能量（$k_BT$），体系通常会被困在初始结构附近的局部能量凹陷中，很难自发跨越显著的能量屏障。
3. 捕获稀有事件（Rare Events）与过渡路径 许多对生物物理和生化功能至关重要的构象转变属于“稀有事件”。增强采样可以特异性地提高系统在过渡路径和高能量能面区域的采样概率，从而直接观察到分子解离、转运通道开启或构象异构化等过程。
4. 准确计算自由能景观与热力学状态差异 热力学量（如结合自由能、构象转变自由能差）依赖于足够充分的相空间采样。常规 MD 采样不足会导致统计误差极大甚至无法收敛，而增强采样通过在反应路径上获取充分的概率分布，可用于精确重构自由能面（Free Energy Surface, FES）并计算状态间的自由能差。

常见的一些增强采样方法#

导向分子动力学（Steered MD, SMD）#

• 物理机制：
  • 沿预设的反应坐标（Reaction Coordinate）向目标分子施加一个外加机械力（如恒力拉伸或移动的虚拟弹簧谐振力）。
  • 这种外力直接迫使体系克服能垒发生解离或移动。
• 特点与应用：
  • 属于非平衡态不可逆过程，能够直观揭示配体进出通道的动力学路径，常用于模拟分子解离、水分子孔道通透，并与原子力显微镜（AFM）的单分子力拉伸实验直接对接。

伞形采样（Umbrella Sampling）#

假设你想测量从山谷 A（状态 A）翻过一座险峻的山峰（高能垒过渡态）到达山谷 B（状态 B）的完整地形高度（自由能变化 PMF）。

• 常规分子动力学的困境（自由漫步的徒步者）：

在常规模拟中，分子就像一个疲惫的徒步者，绝大部分时间都在舒服的山谷 A 或山谷 B 底部打转，几乎不可能自发爬上陡峭险峻的山顶。由于山顶上几乎没有脚印（采样计数为 0），你根本无法知道山峰究竟有多高、过渡路径到底长什么样。

• 伞形采样的解决思路（栓着弹簧绳的分段科考队）：

与其放任分子自由漫步，不如人为介入：

1. 沿登山小路从山谷到山顶一路钉下 20 根木桩（划分成 20 个采样窗口）。
2. 派 20 位科考队员，每人腰上系一根弹性非常好的弹簧绳，分别拴在对应的木桩上。
3. 即使是拴在山顶上的那位队员，想往下溜也会被弹簧绳硬拽回来。这样一来，每一位队员都会在各自木桩周围被约束着走动，从而在悬崖峭壁和山峰各处都留下密密麻麻的采样脚印。
4. 因为我们事先精确知道每根弹簧的弹性系数和拉力大小，模拟结束后，只需在数学上把弹簧的拉力全部扣除，就能完美复原出没有弹簧时山脉的真实海拔起伏。

伞形采样的四个具体执行步骤#

1） 选定路线：定义反应坐标（Reaction Coordinate, $\xi$）#

伞形采样需要先验知道构象变化的大致方向。这个坐标 $\xi$ 可以是：

• 两个结构域之间的距离（例如通道的开合度）；
• 配体中心与受体结合口袋的距离（例如配体的解离脱出）；
• 某个关键残基的二面角。

2） 沿途设桩：划分采样窗口（Windows）#

沿着反应坐标 $\xi$ 的全长（比如从 $0.5\text{ nm}$ 到 $2.0\text{ nm}$），等间距划分出多个窗口（例如每隔 $0.05\text{ nm}$ 设一个参考点 $\xi_0$）

3） 拴上弹簧：施加偏置谐振势（Biasing Potential）#

在每个窗口中分别独立运行一段 MD 模拟，模拟时，给体系额外加上一个抛物线形的人工弹簧势

$$V_{\text{bias}}(\xi) = \frac{1}{2} k (\xi - \xi_0)^2$$

其中：$k$ 是弹簧劲度系数，$\xi_0$ 是当前窗口的参考中心点

当分子试图逃离这个窗口时，势能迅速升高，迫使体系只能在 $\xi_0$ 附近震荡采样， 这保证了哪怕是能量极高的过渡态（山顶），也能被强制收集到成千上万个构象样本

4）数学去偏与拼接：重构真实自由能剖面（WHAM 算法）#

每个窗口采样完毕后，都会得到一个局部的构象分布柱状图。但这些数据都被弹簧“人为扭曲”了。

• 通过加权直方图分析法（WHAM）或 MBAR 算法，利用已知的弹簧势能函数，在数学上把外加的弹簧势精确扣除（去偏，Unbiasing）
• 将各个窗口去偏后的真实概率分布无缝拼合在一起，最终换算出整条反应路径上的平均力势能曲线（PMF / 自由能剖面）

元动力学（Metadynamics）#

假设你是一个被困在深山谷（低能量稳态）里的徒步者，四周是高耸的悬崖（高能垒）。常温下的热运动能量极小，你无论怎么跳跃都翻不过去

元动力学的策略是：

1. 记忆足迹： 你不用费力去硬跳。只要你在山谷底部某个位置停留，脚下就会自动落下一小堆沙子（添加一个微小的高斯偏置能量包）
2. 填平洼地： 随着你在山谷里走动，你踩过的地方沙子越积越高，这个坑被一点一点填平
3. 自发越垒： 当沙子把整座山谷彻底填满到与山顶山脊持平的那一刻，这个原本的“深坑”对你而言已经变成了平地。在微弱的热运动下，你自然而然地漫步越过了山脊，滑落到了隔壁未知的第二个山谷中。
4. 地势复原： 当你把沿途所有遇到的山谷全部用沙子填成了一马平川的大平原时，你一共倾倒了多少沙子、沙堆各处有多厚，取一个负号，就正好严丝合缝地复原出原本山谷有多深、山峰有多高

元动力学的实现路径#

1） 选择集体变量（CVs）#

计算机无法直接对体系成千上万个原子的 $3N$ 维坐标倒沙子，必须将运动降维到 1～3 个关键的集体变量 $s$（比如某个二面角、两个结构域之间的质心距离、或者结合口袋与配体的配位数

2）沉积历史相关的偏置势（History-Dependent Bias Potential）#

在模拟运行过程中，每隔固定的时间步长 $\tau$（例如每隔 1～2 ps），算法会检测当前体系所处的 CV 坐标 $s(t)$，并在此处投放一个局部的高斯形状势能包：

$$V_{\text{bias}}(s, t) = \sum_{t' = \tau, 2\tau, \dots < t} W \exp\left( -\frac{(s - s(t'))^2}{2\sigma^2} \right)$$

• $W$：高斯势能包的高度（Hill Height），即“每次倒多少沙子”
• $\sigma$：高斯势能包的宽度（Hill Width），即“沙堆铺开的范围”
• $s(t')$：历史时刻 $t'$ 时体系所处的坐标

这些外加的高斯势能会像地毯一样累加在分子本身的势能曲面上。外加势能产生的排斥力会将分子不断推离已经访问过的区域，迫使体系去探索那些未曾涉足的高能和未知构象空间

3） 自由能的严格重构（Free Energy Surface, FES）#

经过足够长的时间，当所有的能阱都被填平后，体系在整个集体变量空间内感受到的有效合外力几乎为零，分子开始在整个相空间内进行自由扩散。

此时，外加偏置势的累积总和与体系本征的自由能 $F(s)$ 刚好大小相等、方向相反： $F(s) = -\lim_{t \to \infty} V_{\text{bias}}(s, t) + C$

不需要像伞形采样那样做复杂的加权直方图反卷积（WHAM），把累积的外加高斯势直接倒过来画，就是真实的自由能剖面（FES）