分子动力学模拟原理：三个核心支柱

做自由能计算，说白了就是让原子动起来。你让原子怎么动？它凭什么动？这就是分子动力学要解决的事。

我个人习惯把这部分拆成三个问题：力从哪来、时间怎么推、环境怎么控。搞明白这三块，模拟就跑得稳了。

牛顿运动方程：模拟的发动机

你想想看，每个原子受到一个力，它就会加速。力从哪来？从势能的梯度来。公式很简单：

F = -∇U(r)
a = F / m

U(r) 是势能，∇U 就是它沿着三个方向的导数。F 是力向量，m 是原子质量，a 是加速度。

嗯，这里要注意：势能函数越复杂，求梯度就越贵。GPCR 体系动不动就是几十万个原子，每一步都得算一次梯度，这就是模拟主要卡脖子的地方。

我记得有一次分析 GPCR 构象演化，发现配体进入结合口袋的瞬间，关键残基侧链扭转了约 30 度。这个细节要是没有牛顿方程驱动，根本抓不到。

所以可以这样理解：

• 势能曲面是底层地图
• 梯度给出了下山方向
• 牛顿方程决定了下山速度

积分算法：让时间离散化

理论方程是连续的，但计算机只能算离散步。你需要一个数值方案，把每个时间步的位置和速度算出来。

常用的有两个：Verlet 和 Velocity-Verlet。我推荐直接用后者，稳定且简单。

算法	特点	我在项目中遇到过的坑
Verlet	只用位置信息，速度靠差分算	能量漂移大，尤其体系剧烈运动时
Velocity-Verlet	位置、速度同时更新，对称性好	这是我现在默认的选择
Leap-Frog	速度和位置交错半步	对刚体约束处理不如 VV 方便

Velocity-Verlet 的核心操作就三步：

1. 用当前加速度，更新半步行速度
2. 用当前速度+半速，更新全步位置
3. 在新位置重新算力，再用新加速度把速度补完

写成伪代码就是这样：

v(t + 0.5*dt) = v(t) + 0.5 * a(t) * dt
r(t + dt) = r(t) + v(t + 0.5*dt) * dt
计算 a(t + dt) = -∇U(r(t + dt)) / m
v(t + dt) = v(t + 0.5*dt) + 0.5 * a(t + dt) * dt

温度控制：模拟生理环境

你在 300 K 下做实验，模拟也得体感相同。温度的本质是原子平均动能，所以控温就是调节速度。

主流方法有这几种：

• Berendsen 恒温器：把速度和目标温度做指数松弛。快，但统计系综不严格，适合平衡阶段
• Andersen 热浴：随机碰撞原子，替换速度。引入真正随机性，适合采样，但会破坏动力学连续性
• Nose-Hoover 链：引入虚拟自由度，生成严格 NVT 系综。我做药物设计时，一般用这个

为什么会这样？Berendsen 说白了是个凑合。它虽然收敛快，但在 GPCR 配体结合自由能计算这种对系综精度要求高的场景，容易低估涨落。我曾经尝过苦头：用 Berendsen 跑了 50 ns，势能均值看起来稳，但算自由能时方差偏差 20% 以上。后来换成 Nose-Hoover 链，结果才对上。

实际参数设定可以参考：

; GROMACS .mdp 中的温度控制示例
tcoupl      = Nose-Hoover
tc-grps     = Protein  Non-Protein
tau-t       = 0.5  0.5
ref-t       = 300  300

压力控制：模拟真实体系体积

生理环境不仅恒温，还得恒压（通常是 1 bar）。压力控不好，盒子体积乱变，计算的自由能就废了。

方法对比：

方法	适用场景	避坑
Berendsen 压浴	平衡阶段快速稳定	产生伪涨落，不适合成品数据采集
Parrinello-Rahman	长时间采样	对快速振荡敏感，耦合时间别小于 1 ps
Andersen 压浴	各向异性体系（如膜蛋白）	盒子形状可变，但参数调起来费劲

我曾经调一个 GPCR 嵌入 POPC 双分子层的体系，用 Parrinello-Rahman 时忘了把盒子各向异性关掉，结果跑着跑着双分子层厚度缩水了 10 Å。后来一步错，重跑 1 μs，血亏。各向异性控制是这样的：

; 正确设置膜体系控压
pcoupl      = Parrinello-Rahman
pcoupltype  = semiisotropic
tau-p       = 5.0  5.0
ref-p       = 1.0  1.0
compressibility = 4.5e-5  4.5e-5

半各向异性（semiisotropic）意味着 xy 面内各向同性，z 方向独立。膜体系必须这样设，其他地方用各向同性就行。

典型氮气-膜-蛋白质系统的工作流

把温压控制串起来，一个完整的模拟阶段长这样：

1. 能量最小化（~50000 步，陡降+共轭梯度）
2. NVT 平衡（~100 ps，保持位置约束，升温到 300 K）
3. NPT 平衡（~1 ns，去掉位置约束，稳定体系密度）
4. 成品 NPT 采样（10-100 ns，保存轨迹和能量）

我建议在第 3 步结束后看一眼盒子尺寸。如果体积抖动超过 2%，说明平衡不够，继续跑。不急，不然算出来的自由能会有系统误差。