分子动力学模拟基础：力场、积分算法、周期性边界条件、温度/压力耦合

要做 GROMACS 输出能量特征提取，得先搞明白 MD 模拟在干什么。说白了，整个模拟就是算原子间的力 → 积分运动方程 → 更新坐标 → 再算力，如此循环。这里面有四个核心概念：力场、积分算法、周期性边界条件、温度/压力耦合。我当年刚接触时，觉得这些太“理论”，结果第一个体系能量直接发散……嗯，从那以后我再也不跳过了。

力场：分子世界的“物理规则”

力场就是一组数学函数和参数，描述原子核之间的相互作用。GROMACS 里最常用的有 OPLS-AA、CHARMM36、AMBER99SB-ILDN 等，每种适合不同的体系——比如蛋白质模拟我偏向用 AMBER，膜体系用 CHARMM 稍多。力场能量分成键合项和非键合项：

E_total = E_bond + E_angle + E_dihedral + E_LJ + E_coulomb

• 键伸缩 / 角弯曲：谐波近似，弹簧系数决定刚度
• 二面角：用傅里叶级数描述旋转势垒
• 范德华：Lennard-Jones 6-12 势，描述色散和排斥
• 静电：库仑相互作用，通常用 PME 处理长程项

积分算法：谁让原子动起来？

力算出来了，怎么更新位置和速度？GROMACS 默认用 蛙跳 (leap-frog) 算法，速度和位置相差半个时间步长：

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

你想想看，这样只需要存储一套速度中间值，计算效率高。另一种常用的是 Velocity-Verlet，GROMACS 也支持，但蛙跳用得更多。

周期性边界条件：有限盒子模拟无限系统

我们不可能模拟几亿个原子，所以取一小块区域（通常是个立方体或菱形十二面体），然后四周复制镜像。这就是 周期性边界条件 (PBC)。原子跑出盒子一头，就从另一头回来。同时，计算非键相互作用时要用 最小镜像约定——只和最近的那个镜像算，避免重复计算。

GROMACS 里常用截断半径（cutoff）来处理短程相互作用，一般设 1.0~1.2 nm。长程静电用 PME 网格。我个人的习惯：

• 盒子边长至少 > 2 × cutoff，这样原子不会跟自己镜像直接作用（虽然 GROMACS 会自动排除 1-2/1-3 原子对）
• 水盒子推荐用 editconf -d 1.0 加至少 1.0 nm 的溶剂层

温度/压力耦合：控制模拟环境

真实的实验有恒温恒压环境。MD 模拟通过 恒温器 和 恒压器 来模拟。GROMACS 提供了多种选择：

我的经验：做能量提取时，采样阶段务必用能产生正确系综的耦合器，否则提取的“热力学量”可能毫无意义。我一般这样搭配：

; 平衡阶段（快速稳定）
tcoupl  = V-rescale
pcoupl  = Berendsen

; 采样阶段（严格系综）
tcoupl  = Nose-Hoover
pcoupl  = Parrinello-Rahman

知识结构总览

下面这张图把本章提到的四个基础概念串起来，顺便标出它们对 能量特征提取 的影响路径，方便你后续看输出文件时知道哪些是“关键参数”。