分析梳理--分子动力学模拟的常规步骤五（Gromacs）

作者，Evil Genius

人生有很多挫折的事，每个人都一样，大家能做的就是面对，想要杜绝是不可能的。

今天我们继续分子动力学：NVT和NPT预平衡。

也是分两步走：第一步NVT预平衡。

能量最小化计算可保证我们的初始结构在几何构型和溶剂分子取向等方面都趋于合理化。为了开始真正的动力学模拟，我们必须对蛋白质周围的溶剂和离子进行平衡。如果我们直接进行非限制的动力学模拟，体系会崩溃从而得不到合理的结果。原因在于我们能量最小化时基本上只是优化了溶剂分子自身，而没有考虑溶质。需要将体系置于设定的模拟温度下，以确定溶质(蛋白质)的合理取向。达到正确的温度。模拟时我们使用的温度为300K，它接近于大多数实验条件的温。有些人会在310K下进行模拟，因为这个温度更接近于体温或生理温度。

为什么需要NVT预平衡？

初始构型可能“太冷”或“太热”：无论是从晶体结构、PDB文件构建，还是通过packmol填充得到的初始构型，其原子速度通常为0（如果未指定）或随机分配。如果不加控制，模拟初期会出现剧烈震荡甚至崩溃。

消除不合理的内应力：初始模型可能存在原子间距离过近、键角扭曲、二面角不合理等情况，导致巨大的排斥力。NVT下，通过控温器将动能逐步引入，体系可以自然地“松弛”掉这些局部应力。

逐步升温：有些体系（如蛋白质）不能直接从0 K跳到300 K，否则会因热冲击而变性。NVT预平衡可以通过逐步升温（退火）或直接使用Berendsen/Bussi控温器温和地达到目标温度。

NVT预平衡的核心控制参数

在NVT模拟中，主要控制温度恒定。常用的控温方法有：

实践建议：

先用 Berendsen 快速消除温度偏差（例如100 ps），然后用 Velocity-rescaling 或 Nosé-Hoover 进行后续NVT平衡。

NVT预平衡的典型时长

小分子溶液体系：100 ps - 500 ps 通常足够。

生物大分子（蛋白质、核酸）：500 ps - 2 ns 常见。有时需要更长的平衡（如5 ns）让侧链和溶剂充分弛豫。

周期性固体/界面体系：50 ps - 200 ps 通常足够，因为结构初始就较为合理。

一个典型的NVT预平衡的mdp文件

; NVT equilibration
title       = NVT equilibration 
define      = -DPOSRES        ; 可选项：对重原子（如蛋白质骨架）施加位置限制，防止剧烈形变
integrator  = md              ; 分子动力学步进
dt          = 0.002           ; 2 fs 时间步长
nsteps      = 50000           ; 总步数 = 时间 / dt，例如 100 ps = 50000步（若dt=0.002）
nstxout     = 5000            ; 坐标输出频率（低输出以节省空间）
nstvout     = 5000            ; 速度输出频率
nstenergy   = 1000            ; 能量输出频率

; Temperature coupling
tcoupl      = V-rescale       ; 或 berendsen（初期），后切换为 v-rescale
tc-grps     = System          ; 整个体系一起控温；也可分蛋白、溶剂组
tau_t       = 0.1             ; 时间常数，典型值0.1-0.5 ps
ref_t       = 300             ; 目标温度 300 K

; Pressure coupling - 无，NVT下不开启
pcoupl      = no

; 其他：PME长程静电，Verlet截断等标准设置

位置限制：如果体系包含柔性大分子（如蛋白质），NVT初期往往建议对重原子施加较弱的位置限制（force constant ~1000 kJ/mol/nm²），让溶剂先弛豫，然后逐渐减弱或移除。

时间步长：使用2 fs（含氢键约束时）。如果体系初始能量极高，可先进行能量最小化，再开始NVT。

输出频率：预平衡阶段不需要高频输出轨迹，只需监控温度、势能和总能量是否稳定。

如何判断NVT预平衡是否完成？

检查输出文件（如GROMACS的.log或energy.xvg）：

温度：围绕目标温度（如300 K）平稳波动，无漂移趋势。

势能：稳定在某个平均值附近，不再持续下降（意味着应力已释放）。

总能量：在NVT系综下，总能量（动能+势能）应基本守恒（小范围涨落），不应有系统性上升或下降。若使用Berendsen控温，总能量可能略有漂移，这是正常的，但不应剧烈变化。

RMSD（对大分子）：相对于初始结构的RMSD快速上升后趋于平台，表明结构已弛豫。

我们来运行一下，首先产生nvt预平衡的tpr文件

gmx grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr -r em.gro

运行NVT预平衡

gmx mdrun -deffnm nvt

我们来看一下生成的文件，基于温度进行分子的重新排布

绘图看看能量变化，满足需求。

接下来第二步：NPT预平衡

核心目的是在保持温度恒定（NVT）的基础上，引入正确的压力（通常是1 bar），让体系的密度调整到合理的值，为最终的生产模拟（NPT或NVT）准备好一个压力稳定、密度正确的平衡体系。

为什么NPT预平衡如此重要？

调整密度至实验值：NVT预平衡结束时，体系的盒子体积和密度是初始设定的（比如来自溶剂填充）。这个密度很可能不是该力场模型下、目标温度压力对应的正确密度。NPT模拟允许盒子体积变化，从而找到正确的平衡密度。

释放残余应力：如果初始溶剂填充或离子放置不够完美，体系中可能存在局部压力不均。NPT模拟通过耦合恒压器，能让整个体系均匀地弛豫。

准备生产环境：大多数生物分子模拟（如蛋白质-配体结合）在恒温恒压（NPT）系综下进行，以模拟真实的生理环境。NPT预平衡确保了体系在进入生产模拟前，压力、密度和体积已经稳定。

核心控制参数（与NVT的主要区别）

NPT模拟比NVT多了一个压力耦合（恒压器），这是关键。

首先我们准备npt.mdp文件（下面是示例）

; NPT equilibration
title       = NPT equilibration
define      = -DPOSRES          ; 可选：继续对蛋白质骨架施加较弱的位置限制，让溶剂和体积先弛豫
integrator  = md
dt          = 0.002             ; 2 fs 时间步长
nsteps      = 100000            ; 200 ps (0.002 * 100000)
nstxout     = 5000              ; 坐标输出频率（每10 ps）
nstvout     = 5000
nstenergy   = 1000              ; 能量输出频率（每2 ps）
nstlog      = 1000
nstxtcout   = 5000              ; 压缩轨迹输出频率（推荐）
xtc-precision = 1000

; 温度耦合 (同NVT)
tcoupl      = V-rescale
tc-grps     = System
tau_t       = 0.1
ref_t       = 300

; 压力耦合 (新增关键部分)
pcoupl      = C-rescale         ; 预平衡首选。之后生产模拟可改用 Parrinello-Rahman
pcoupltype  = isotropic         ; 各向同性
tau_p       = 2.0               ; 压力时间常数，典型值2-5 ps
ref_p       = 1.0               ; 目标压力 1 bar
compressibility = 4.5e-5        ; 水的等温压缩系数，4.5e-5 bar^-1

; 键约束 (同NVT)
constraints    = h-bonds
constraint-algorithm = lincs

我们运行一下，首先生成npt.tpr文件

gmx grompp -f npt.mdp -c nvt.gro -p topol.top -o npt.tpr -r nvt.gro

运行NPT预平衡

gmx mdrun -deffnm npt

如何判断NPT预平衡是否完成？

你需要分析输出的 npt.edr 能量文件，主要看以下几个量是否达到稳定平台：

看一下生成的结果

生活很好，有你更好