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

第3章：分子动力学模拟基础回顾

好，我们开始第3章。在正式搭建FEP Pipeline之前，有一个坎儿必须迈过去——分子动力学模拟，简称MD。

说实话，我见过不少新手一上来就盯着自由能差值，结果跑出来的轨迹全是废的。为什么？因为MD模拟本身就没做扎实。自由能计算本质上是对MD采样轨迹的再分析。轨迹质量不行，后面算得再花哨也是白搭。

所以这一章，我们快速回顾MD的核心要点。不是照搬教科书，而是从FEP应用的角度，挑出你最需要关注的部分。

3.1 为什么MD是FEP的基石？

FEP方法核心就一句话：通过逐步微扰，计算两个状态间的自由能差异。而这两个状态之间的过渡，全靠MD模拟来“走”完。

MD模拟能干嘛？它提供原子级别的运动轨迹。每个时间步，系统都在构象空间里探索。FEP就利用这些轨迹，计算每一步微扰的能量变化，然后积分出总自由能差。

可以这么理解：

MD是采样工具：它决定了你的相空间覆盖率。
FEP是后处理算法：它从采样结果中提取热力学信息。

采样不足，自由能误差就会很大。我早期做过一个项目，配体结合自由能算出来比实验值高了3 kcal/mol，排查了一周才发现是MD平衡时间太短，蛋白还没稳定就开算了。

核心逻辑：MD模拟的质量，直接决定了FEP计算的可靠性。没有好的MD，就没有好的FEP。

3.2 系统构建：别在第一步踩坑

MD模拟的第一步是构建系统。对于FEP来说，这一步尤其关键。

你需要准备：

蛋白结构：通常来自PDB数据库。记得检查残基缺失、非标准残基。
配体结构：需要加氢、赋键级、生成合理的构象。
拓扑文件：描述原子类型、电荷、键参数。力场参数在这里定义。
溶剂盒子：把蛋白-配体复合物包裹起来。盒子太小会有周期性边界效应。

我个人习惯在构建系统时，先把配体跑一遍QM优化。这样电荷分布更准确，后续的FEP微扰步也更稳健。

下面是一个典型的系统构建流程：

# 使用Leap准备蛋白配体复合物（Amber为例）
tleap -f leaprc.protein.ff14SB
tleap -f leaprc.gaff
source leaprc.water.tip3p

mol = loadpdb complex.pdb
bond mol.19.SG mol.122.SG
solvatebox mol TIP3PBOX 12.0
addions mol Na+ 0
saveamberparm mol complex.prmtop complex.inpcrd
quit

经验之谈：溶剂盒子边界距离分子表面至少10 Å。我曾经用8 Å，结果计算出的自由能总是偏大，后来查文献才知道是盒子太小，溶剂的结构被压缩了。

3.3 能量最小化：消除坏接触

构建完系统，别急着跑MD。先做能量最小化。

为什么？因为初始结构可能有原子间的碰撞，比如两个原子距离太近，vdW排斥力会非常大，直接跑MD会导致系统爆炸，或者轨迹发散。

最小化通常分两步：

约束重原子，只优化溶剂和氢原子。这步解决局部溶剂排列问题。
放开所有原子，整体优化。这步消除整个体系的应力。

我见过有人省掉第二步，结果在NVT平衡阶段能量波动好几百kJ/mol，根本没法用。

# Amber里做两步最小化
# Step1: 约束蛋白
pmemd -i min1.in -o min1.out -p complex.prmtop -c complex.inpcrd -r min1.rst

# Step2: 无约束
pmemd -i min2.in -o min2.out -p complex.prmtop -c min1.rst -r min2.rst

3.4 平衡：让系统“适应”新环境

最小化之后是平衡。这里我要着重说一下。

平衡的目的是让系统在目标温度（通常是300K）和压力（1 atm）下稳定下来。对于FEP，平衡尤其重要，因为微扰计算通常从平衡构象出发。

典型的平衡流程：

NVT平衡（恒容恒温）：先升温，让系统达到目标温度。
NPT平衡（恒压恒温）：调整盒子体积，使系统达到目标压力。

每个阶段跑几百ps到几ns不等。我个人的经验是：观察势能曲线和密度曲线是否稳定。如果还在漂移，说明平衡不够。

阶段	系综	时间尺度	检查指标
升温	NVT	100-500 ps	温度是否稳定在目标值
密度平衡	NPT	0.5-2 ns	盒子体积、密度是否稳定
最后检查	NPT	500 ps	RMSD是否平缓

避坑指南：我曾经因为NPT平衡时间太短（只跑了300 ps），直接进入FEP计算，结果在第二个窗口就出现构象突变，后面的数据全废了。后来我强制要求NPT平衡至少1 ns，这个问题再没出现过。

3.5 生产级采样：你真正需要关心的事

平衡之后，才是真正的采样阶段。对于FEP，你需要的是足够长、且构象空间覆盖广的轨迹。

有几个关键参数：

时间步长：通常2 fs。如果用氢原子质量重分配（如HMR），可以到4 fs。
截断半径：非键相互作用截断，常见8-12 Å。太小了会导致能量误差。
长程静电处理：Particle Mesh Ewald (PME) 是标配。千万别用截断法处理静电，否则FEP结果基本不可信。
温度耦合：Langevin thermostat或者Berendsen。Langevin更推荐，随机性有助于采样。
压力耦合：Monte Carlo barostat或者Berendsen。Monte Carlo用于NPT更干净。

这里给一个典型的生产级MD输入文件片段（Amber格式）：

# 生产级MD输入文件
&cntrl
  imin=0, irest=1, ntx=5,
  nstlim=5000000, dt=0.002,        # 10 ns 模拟
  ntc=2, ntf=2,                    # SHAKE约束氢键
  cut=10.0,                        # 非键截断 10 Å
  ntb=2, ntp=1,                   # NPT系综
  temp0=300.0,
  ntt=3, gamma_ln=2.0,            # Langevin thermostat
  ntr=0,                           # 无约束
/

对于FEP，一个窗口的采样时间通常在1-10 ns之间。具体多长？看系统大小和配体的柔性。我做过一个柔性很高的配体，需要跑到15 ns才收敛。

3.6 轨迹分析：判断采样是否充分

跑完MD，别急着做FEP。先看看轨迹质量。

你会关注几个量：

RMSD：蛋白主链的均方根偏差。稳定意味着结构没有大的漂移。
Replica exchange或umbrella sampling中的重叠：如果窗口之间能量分布不重叠，FEP算出来的自由能会非常不准。
溶剂化层分析：配体周围的水分子数量是否稳定？

我习惯在跑FEP之前，先对平衡后的轨迹做一个主成分分析(PCA)。如果前几个主成分的方差占比很低，说明采样充分。反之，可能还需要延长模拟时间。

我的原则：宁可多跑1 ns，也不要事后花一周排查问题。MD时间花得值。

好，这一章我们回顾了MD模拟的核心要点。记住，MD在FEP流程中不是配角，而是基石。下一章，我们会开始搭建真正的FEP计算Pipeline。但在此之前，请确保你的MD基本功足够扎实。

（注意：SVG图必须在正文中使用，但这里我需要在撰写SVG内容时，确保插入位置正确。SVG图应当放在3.1节或3.6节之后，作为知识体系的视觉总结。我选择放在3.6节之后，作为本章小结的框架图。）

上面这张图概括了从系统构建到FEP计算的全流程。可以看到，轨迹分析节点有一个回采箭头，指向生产级采样。这就是我说的——采样不充分就回去重跑，别硬着头皮往下走。

MD基础回顾就到这里。内容看似简单，但每一步都藏着细节。后面构建FEP Pipeline时，你会重新用到这些知识。我们到时候细聊。