3. 分子动力学模拟简介

做自由能景观图之前，你得先搞清楚分子动力学模拟是怎么回事。说白了，MD 就是用计算机把原子在时间上的运动「拍」下来。我们不是真的做实验，而是让原子遵循牛顿定律去跑轨迹。

我个人习惯把 MD 类比成拍电影。你给每个原子一个初始位置和速度，然后告诉电脑它们之间怎么相互作用——这就是力场。电脑每算一小步，就把所有原子的新位置记下来。连续拍几千几万帧，就得到一条轨迹。

3.1 MD 的核心思想

MD 的根基是牛顿第二定律：

F = m · a
a = d²r/dt²

加速度由力决定，力又由势能梯度给出：

F = -∇U(r)

所以只要知道势能函数 U(r)，就能算出每个原子受到的力，进而更新位置和速度。重复这个步骤，轨迹就出来了。道理很简单，但真正做起来全是坑。

MD 模拟的基本循环：

1. 由当前坐标计算势能 U 和力 F
2. 用 F 更新原子速度
3. 用速度更新原子坐标
4. 把新坐标写入轨迹文件
5. 回到第 1 步，直到达到设定的步数

这个循环每步只做一次。步长通常取 1~2 飞秒。你想想看，模拟 1 纳秒就要跑 50 万步。所以 MD 的计算量很大，尤其体系大了以后。

3.2 力场——模拟的「剧本」

力场定义了原子间怎么相互作用。没有好力场，MD 就是瞎跑。常见的力场有 CHARMM、AMBER、OPLS 等，各有各的适用场景。

力场通常写成这种形式：

U = U_bond + U_angle + U_dihedral + U_nonbonded

U_bond = Σ k_b (r - r₀)²
U_angle = Σ k_θ (θ - θ₀)²
U_dihedral = Σ k_φ [1 + cos(nφ - δ)]
U_nonbonded = Σ (4ε[(σ/r)¹² - (σ/r)⁶] + q_i q_j / (4πε₀ r))

前几项是键长、键角、二面角的简谐振动，最后一项是范德华（Lennard-Jones）和静电（库仑）相互作用。选力场时要注意：你的分子体系里有没有金属离子？有没有特殊共价键？我遇到过有人用通用力场跑含配位键的蛋白，结果配位键直接断了，轨迹完全不可用。

实际上，做自由能景观图最常用的是全原子力场，因为精度比较靠谱。粗粒化虽然跑得快，但丢失了原子细节，做自由能图容易失真。

3.3 积分算法——时间怎么走

MD 里最经典的积分方法是 Verlet 算法。其实原理很简单：用当前位置和之前位置，加上当前加速度，就能推出下一步位置。

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

后来人们觉得这样算速度不方便，就改进了 Velocity Verlet：

1. 计算 r(t+Δt) = r(t) + v(t) · Δt + ½ a(t) · Δt²
2. 由新坐标计算 a(t+Δt)
3. 更新 v(t+Δt) = v(t) + ½ [a(t) + a(t+Δt)] · Δt

这个算法时间可逆、能量守恒好，是目前 GROMACS、NAMD 等主流软件的标准做法。要理解自由能景观，你得知道轨迹是由这些算法一步步算出来的，不是猜的。

3.4 模拟的典型流程

一段完整的 MD 模拟不是点一下就跑。你需要按顺序做：

1. 准备结构：从 PDB 或自己建模，补氢、加溶剂、加离子中和电荷。
2. 能量最小化：消除原子间的不合理接触。我见过有人跳过这步，结果一跑就崩溃。
3. 升温：从 0 K 慢慢加热到目标温度（比如 300 K）。升太快会导致热不均匀。
4. 预平衡：固定蛋白骨架，让溶剂先在周边平衡。我习惯做两轮：先用 NVT，再用 NPT。
5. 生产模拟：放开所有约束，跑长轨迹用于分析。

每一步都有对应的控制参数。下面是一个典型的 GROMACS mdp 片段，用于生产过程：

integrator               = md        ; 使用 Velocity Verlet
dt                       = 0.002     ; 步长 2 fs
nsteps                   = 50000000  ; 总步数，对应 100 ns
nstxout-compressed       = 5000      ; 每 10 ps 输出一次坐标
nstvout-compressed       = 5000      ; 每 10 ps 输出一次速度
tcoupl                   = V-rescale ; 恒温器
ref-t                    = 300       ; 目标温度
pcoupl                   = Parrinello-Rahman ; 恒压器
ref-p                    = 1.0       ; 目标压力 1 bar

3.5 从轨迹到自由能景观

MD 模拟最终输出的是轨迹文件，里面存了每一步每个原子的坐标和速度。你可以从中提取各种反应坐标（比如两个残基之间的距离、某个二面角、RMSD），然后统计这些坐标的分布。

自由能景观图本质上就是：F(ξ) = -k_BT ln P(ξ)，其中 P(ξ) 是反应坐标 ξ 的概率密度。这个概率密度就是从轨迹上直接数出来的。

但常规 MD 很难翻越能垒，所以需要增强采样的方法（如伞形采样、元动力学）才能得到完整的自由能图。这些我会在后面的章节展开。

这里先记住一句话：MD 轨迹是自由能景观图的数据原料，没有好的 MD 采样，再漂亮的图也是假的。

3.6 一个完整的 MD 模拟流程框架

下面这张图把整个流程串了起来。你看，从初始结构开始，经过力场参数化、能量最小化、升温平衡，最后进入生产模拟。生产模拟的输出就是轨迹，而轨迹又能通过反应坐标统计得到自由能景观图。