📚 课程目录第4章分子动力学模拟基础

4. 分子动力学模拟基础：从牛顿方程到时间步长

做构象采样，咱们得先搞清楚分子动力学（MD）到底在算什么。

说白了，MD模拟就是在解一个巨大的微分方程组。每个原子的运动都服从牛顿第二定律。你用计算机一步步推演，就能看到整个分子体系随时间演化的轨迹。这条轨迹里，就藏着我们想要的构象信息。

我刚开始接触MD时，总觉得这跟解高中物理题差不多。后来发现，嗯，没那么简单。数值积分算法的选择、时间步长的设定，都会直接影响你的采样结果。今天咱们就把这些基础夯实了。

4.1 牛顿运动方程：MD模拟的物理根基

MD模拟的核心就两个方程：

Fi = mi * ai
ai = d²ri / dt²

Fi 是作用在第 i 个原子上的合力，mi 是它的质量，ai 是加速度。ri 是原子的位置向量。

你想想看，有了力，就能算加速度。有了加速度，就能更新速度和位置。就这么一直迭代下去。

问题在于，体系的势能函数特别复杂。它包含键伸缩、键角弯曲、二面角扭转，还有非键相互作用（范德华力、静电势）。所以合力 Fi 的计算，是整个模拟中最耗时的部分。

我的习惯：开始一个新项目前，我会先跑几皮秒的NVE系综（微正则系综）测试。看看总能量是否守恒。如果能量漂移超过千分之一，那大概率是力场参数或时间步长有问题。这个诊断方法我用了很多年，非常管用。

4.2 积分算法：如何在计算机里解方程

牛顿方程是连续的，但计算机只能处理离散的数据。我们需要一个“积分器”，把连续运动拆成一个个小时间片。

积分器的好坏，直接决定模拟的稳定性和效率。我个人最常用的两个算法是 Verlet 和 Leapfrog。

4.2.1 Verlet 算法

Verlet 算法的思路很直接。它用当前位置 r(t) 和上一时刻位置 r(t-Δt)，加上加速度 a(t)，来推算下一时刻位置 r(t+Δt)。

r(t+Δt) = 2r(t) - r(t-Δt) + a(t) * Δt²

这个公式里没有速度项。速度是通过相邻时刻的位置差算出来的：

v(t) = [r(t+Δt) - r(t-Δt)] / (2Δt)

Verlet 算法的优点是时间可逆，能量守恒性好。但它有个小毛病：处理速度时有点“滞后”。在需要实时调控温度或计算动能的时候，会多一步计算。

4.2.2 Leapfrog 算法

Leapfrog 是我个人更偏爱的一种算法。它的名字很形象——速度和位置像青蛙一样交替跳跃前进。

v(t+Δt/2) = v(t-Δt/2) + a(t) * Δt
r(t+Δt) = r(t) + v(t+Δt/2) * Δt

你看，速度在 t+Δt/2 时刻更新，位置在 t+Δt 时刻更新。两者错开半个时间步。

速度 Verlet 算法（最常用）

1. 用当前位置和速度推算中间位置
r(t+Δt) = r(t) + v(t) * Δt + 0.5 * a(t) * Δt²

2. 计算新位置下的力 F(t+Δt) 和加速度 a(t+Δt)

3. 更新速度
v(t+Δt) = v(t) + 0.5 * [a(t) + a(t+Δt)] * Δt

这个版本的实现非常对称，速度也同步，是做约束动力学和恒温器耦合时的好选择。我目前所有的生产级模拟都用的速度 Verlet。

4.2.3 算法对比

算法	时间可逆性	能量守恒	速度获取	适用场景
Verlet	优秀	优秀	需后处理	NVE 系综, 基础研究
Leapfrog	优秀	优秀	错半步长	NVE, NVT 系综
速度 Verlet	优秀	优秀	同步获取	通用, 建议首选

4.3 时间步长：一个看似简单却至关重要的问题

时间步长 Δt 的选择，是每个 MD 模拟者都绕不开的决策。

你想想看，时间步长太小，模拟时间不够，采样不充分。时间步长太大，体系能量会漂移，甚至直接“炸掉”。

为什么会“炸掉”？ 因为分子的最高振动频率（比如 C-H 键伸缩，周期约 10 飞秒）决定了时间步长的上限。根据采样定理，你的步长必须低于最高频率周期的一半。实践经验告诉我们，时间步长通常取最高频率周期的 1/5 到 1/10。

我遇到过最典型的事故：一个新手用 GROMACS 跑模拟，把时间步长设成 5 fs，还挂在日志文件里大摇大摆地说要跑 100 ns。

结果模拟刚开始几百步，能量直接爆表，体系原子飞到了盒子对面。嗯，这就是典型的“不稳定性爆炸”。

4.3.1 常见步长选择指南

模拟体系	建议时间步长	理由
全原子模拟 (含H原子)	1.0 - 2.0 fs	必须解析C-H键的快速振动
全原子 (应用SHAKE/RATTLE约束)	2.0 fs	约束住H原子振动，步长可翻倍
粗粒化模拟	5.0 - 20.0 fs	少了高频振动，步长可以很大
隐式溶剂模型	1.0 - 2.0 fs	溶剂的摩擦和随机力改变了稳定性

反复校验：我曾经吃过一次亏，在优化一个蛋白质配体体系时，把时间步长从 1 fs 改到 2 fs，感觉模拟跑快了一倍。但事后分析 RMSD 发现，配体的构象分布跟 1 fs 跑出来的结果有明显差异。后来我花了更大的代价重新补跑了 5 条轨迹。所以我现在的习惯是：在项目初期，至少用两个不同的时间步长跑 200 ps 的预测试验，检查关键物理量（势能、回旋半径、主链 RMSD）的一致性。

4.4 知识体系结构图

我把本章的核心逻辑总结成一张图，方便你快速理清思路：

从这张图上你可以看到，牛顿方程是源头。我们通过积分算法把它离散化。算法的“准确性”和“稳定性”直接决定了你能选多大的时间步长。而时间步长，反过来又通过采样效率影响整个模拟的可信度。

4.5 实践建议

实际做项目时，我一般按这个流程来：

确定约束策略：是否使用 SHAKE 或 LINCS 约束键长。如果用了，步长可以安全地设在 2 fs。
设置测试模拟：跑一段 10-20 ps 的 NVT 平衡。
监控能量：盯着 Potential Energy 和 Total Energy 的波动。如果出现持续上升或下降趋势，那大概率是步长太大了。
检查温度振荡：温度应该是围绕目标值上下波动的，波动幅度跟体系大小有关。如果温度异常飙升，赶紧停下来检查。

一个小技巧：我习惯在输入文件里启用 nstxout 和 nstvout 参数，每几百步输出一次坐标和速度的轨迹。这样如果模拟中途挂了，我还有前面半段的轨迹可用，不至于归零重来。

时间步长的选择不是一成不变的。我见过有人用 多重时间步长 技术，把短程相互作用和长程相互作用的更新时间分开。比如短程用 2 fs，长程用 6 fs，能进一步提速。但那是更进阶的话题了，咱们先把基础打牢。

嗯，记住一句话：稳定，永远比速度重要。一个跑得稳但慢 20% 的模拟，产出的数据比一个快但不可靠的模拟有价值得多。