📚 课程目录第4章模拟准备

准备开跑：能量最小化与NVT、NPT平衡的批处理

上一章咱们把体系搭建好了，拓扑文件也有了。但你别急着直接跑生产模拟——那会出大问题。我刚开始入行那会儿，有次把没平衡过的体系直接扔进生产MD，跑了三天发现蛋白都散架了。老板脸都绿了。

说白了，模拟准备的核心就三件事：能量最小化、NVT平衡、NPT平衡。这三步一步都不能省。今天咱们就聊聊怎么写批处理脚本，把这套流程自动化跑通。

核心原则：先消除原子间不合理接触，再让溶剂和温度达到平衡，最后调好压力。顺序错了，体系必爆。

为什么非要做能量最小化？

你构建体系的时候，水分子加进来了，但紧密排列的原子之间可能有严重碰撞。这些不合理接触会在分子动力学第一步就产生巨大的排斥力——模拟直接就崩了。我曾经见过一个体系，用 steep 算法跑了两千步，最大力还在 5000 kJ/mol/nm 以上，结果检查发现有个水分子嵌到了蛋白侧链里。

GROMACS 提供了两种主要的 EM 算法：steep（最速下降法）和 l-bfgs（拟牛顿法）。前者适合粗调，后者精调收敛更快。我的习惯是先跑 5000 步 steep，再切到 l-bfgs 直到最大力小于 10 kJ/mol/nm。

来看一个典型的 EM 参数设置文件：

; minim.mdp - 能量最小化
integrator               = steep      ; 最速下降
nsteps                   = 5000       ; 步数
emtol                    = 10.0       ; 力收敛阈值 (kJ/mol/nm)
emstep                   = 0.01       ; 初始步长 (nm)
nstxout                  = 100        ; 每100步输出坐标
nstlog                   = 100        ; 每100步输出日志
nstenergy                = 100        ; 每100步输出能量

小窍门：如果用 steep 收敛不理想，可以复制一份 .mdp，把 integrator 改成 l-bfgs，用前一步的 .gro 做输入继续跑。不用新建体系，省时间。

NVT 和 NPT：温度与压力的双人舞

能量最小化之后，体系只是没了碰撞，但还没处于热力学平衡态。NVT 平衡负责把体系的温度调到目标值。我一般选 300 K，用 V-rescale 耦合算法。温度耦合的时间常数设 0.1 ps 就够，太快会让体系震荡，太慢起不到调节作用。

NPT 平衡更关键。它调整的是盒子的尺寸，让压力的平均值稳定在 1 bar。这里要注意：必须用各向同性的压力耦合，除非你做的是膜蛋白或各项异性体系。Parrinello-Rahman 算法虽然收敛慢，但更准确，适合后续生产模拟。Berendsen 算法快，但算出来的压力涨落不对，我只用它做预平衡。

; nvt.mdp - NVT 平衡
integrator               = md        ; 分子动力学
nsteps                   = 50000     ; 100 ps (dt=0.002 ps)
dt                       = 0.002     ; 时间步长 (ps)
nstxout                  = 500       ; 每1 ps输出坐标
nstlog                   = 1000
nstenergy                = 1000

; 温度耦合
tcoupl                   = V-rescale
tc-grps                  = Protein  Water_and_ions
tau_t                    = 0.1       0.1
ref_t                    = 300       300

; 压力耦合（NVT阶段关掉）
pcoupl                   = no

; 约束相关
constraints              = h-bonds   ; 约束氢键
constraint-algorithm     = LINCS

; 周期性边界条件
pbc                      = xyz

; npt.mdp - NPT 平衡
; 继承 NVT 设置（省略重复部分）
; 只修改以下条目

pcoupl                   = Parrinello-Rahman
pcoupltype               = isotropic
tau_p                    = 2.0
ref_p                    = 1.0        ; 参考压力 (bar)
compressibility          = 4.5e-5    ; 水的压缩率 (1/bar)

; 建议跑更长的时间才能稳定
nsteps                   = 250000    ; 500 ps

; 生成速度
gen_vel                  = yes
gen_temp                 = 300
gen_seed                 = -1        ; 随机种子

注意！NPT 平衡的步数至少 200 ps。我曾经偷懒只跑了 50 ps，结果盒子的密度还在一直下降，等发现的时候已经浪费了两天做生产模拟。检查 log 文件中的密度和压力波动，稳定了再往下走。

批处理脚本：一次跑完三件套

手动一步步敲 gmx 命令太费时间还容易出错。我习惯写一个 bash 脚本，把三个阶段的预处理和运行串起来。下面是我的参考模板：

#!/bin/bash
# prepare_md.sh - 三阶段模拟准备批处理

# 设置环境变量
GMX=gmx_mpi
TOP=../topol_amber.top
STRUCT=../step6.0_minimization.gro

# 阶段1：能量最小化
echo "运行能量最小化..."
$GMX grompp -f minim.mdp -c $STRUCT -p $TOP -o em.tpr
$GMX mdrun -deffnm em -v

# 检查EM是否收敛
tail -5 em.log | grep "F-max"
if [ $? -ne 0 ]; then
    echo "EM 可能未收敛，建议检查日志。"
fi

# 阶段2：NVT平衡
echo "运行NVT平衡..."
$GMX grompp -f nvt.mdp -c em.gro -p $TOP -o nvt.tpr
$GMX mdrun -deffnm nvt -v

# 阶段3：NPT平衡
echo "运行NPT平衡..."
$GMX grompp -f npt.mdp -c nvt.gro -p $TOP -o npt.tpr
$GMX mdrun -deffnm npt -v

echo "模拟准备完成！可用结构：npt.gro"

脚本简单，但有个细节：每次 grompp 之后都要检查生成的 .tpr 文件。如果输入参数有问题，grompp 会报 WARNING 甚至 ERROR。我个人习惯加一行 check_tpr.sh 来验证：

# 快速检查 tpr 文件
echo "检查 npt.tpr 的内容..."
echo "5" | $GMX dump -s npt.tpr | head -30
# 或者用更现代的方法：
gmx check -f npt.tpr -n -c

表：三阶段参数对比一览

阶段	核心积分器	典型步数	压力耦合	温度耦合	检查指标
EM（steep）	steep	5000	无	无	最大力 < 10 kJ/mol/nm
EM（l-bfgs）	l-bfgs	~2000	无	无	势能平坦
NVT	md	50000 (100 ps)	无	V-rescale	温度~300 K
NPT	md	250000 (500 ps)	Parrinello-Rahman	V-rescale	密度、压力稳定

最终检查清单：

EM 日志中 F-max 小于 10，且不再下降
NVT 的温度 RMSD 不超过 10 K
NPT 的密度和压力在最后 100 ps 没有明显趋势
所有 .tpr 文件都能被 gmx dump 正常读取
盒子的形状没有扭曲（pbc 正常）

做完这三步，你的体系才算真正准备好跑生产模拟。别急，一步一个脚印，后面才能跑得稳。我每次给别人审核模拟流程，第一件事就是看他们有没有完整的三阶段日志。缺了任何一步，我都会让他们回去补上——因为那是最常见的翻车原因。