GROMACS 正常模、Hessian 与谐振量子修正

下载本组输入、真实参考结果和绘图脚本，解压后进入资料包目录。inputs 是可以修改的输入，reference 是鸢算科技集群已经计算得到的参考结果。每次运行在 work 下建立新目录，路径写入 LAST-RUN.txt。计算请在已经分配的资源中执行；整组可用 sbatch run.sbatch 提交。

先跑通，再修改模型#

本例是三个原子的经典柔性 TIP3P 类型教学水分子。它适合看懂正常模工作流，频率不等同于 DFT 或实验水分子光谱。使用以下入口可从最小化一直运行到正常模及投影监测：

bash run.sh modes
cat LAST-RUN.txt

先阅读输入，再调整。inputs/modes/start.gro 给出初始坐标，topol.top 定义键、角和原子质量。这里没有 SETTLE，也不约束 O—H 键，否则无法得到这三个完整内部振动模。

1. 选择双精度环境#

module purge
module load gromacs/2026.4-cpu-double
gmx_d --version

这是固定 v2026.4 源码构建的 CPU 双精度入口，程序显示版本 2026.4-dev；GNU 13.3、thread-MPI/OpenMP、内部 FFTPACK/BLAS/LAPACK，未连接 PLUMED。普通 MD 与 GPU 运行仍按环境教程选择相应模块。下面的手动操作可在独立目录完成：

mkdir modes-manual
cp inputs/modes/* modes-manual/
cd modes-manual

2. 查看完整拓扑和最小化输入#

topol.top 中的长度单位为 nm，能量为 kJ/mol，质量为 u：

; Explicit flexible TIP3P-type teaching monomer; nm, kJ/mol, e units.
[ defaults ]
1 2 yes 0.5 0.8333333333
[ atomtypes ]
OW 8 15.9994 0 A 0.315061 0.636386
HW 1 1.008 0 A 0 0
[ moleculetype ]
SOL 2
[ atoms ]
1 OW 1 SOL OW 1 -0.834 15.9994
2 HW 1 SOL HW1 1 0.417 1.008
3 HW 1 SOL HW2 1 0.417 1.008
[ bonds ]
1 2 1 0.09572 502416
1 3 1 0.09572 502416
[ angles ]
2 1 3 1 104.52 628.02
[ system ]
Flexible water monomer
[ molecules ]
SOL 1

本例盒边长 3 nm，非键截断 1 nm；SOL 2 排除分子内相邻关系对应的非键相互作用。两个键与一个角产生三个内部振动自由度。

最小化输入 em.mdp：

cutoff-scheme = Verlet
nstlist = 1
rlist = 1.0
coulombtype = Cut-off
rcoulomb = 1.0
vdwtype = Cut-off
rvdw = 1.0
pbc = xyz
constraints = none
verlet-buffer-tolerance = -1
nstenergy = 1
nstcalcenergy = 1
integrator = steep
nsteps = 5000
emtol = 0.000001
emstep = 0.001

正常模对驻点和坐标精度敏感，这里使用比普通 MD 最小化更紧的力阈值。对于大体系，应先充分最小化，再研究阈值与低频模的稳定性。

gmx_d grompp -f em.mdp -c start.gro -p topol.top -o em.tpr
gmx_d mdrun -deffnm em -ntmpi 1 -ntomp 2 -nb cpu -bonded cpu -update cpu
printf '0\n' | gmx_d trjconv -f em.trr -s em.tpr -o minimum.g96 -dump 0

最后一步取出 TRR 保存的高精度最终坐标。-dump 0 选择离指定时刻最近的一帧，并输出带原子名称的单帧 G96；本例 EM TRR 只有最终帧，显示的时间标签来自最小化步数，不是 MD 演化时间。不要把最终坐标先写成默认三位小数 GRO，再用于严苛的 Hessian 检查。

3. 计算 Hessian#

完整 nm.mdp：

cutoff-scheme = Verlet
nstlist = 1
rlist = 1.0
coulombtype = Cut-off
rcoulomb = 1.0
vdwtype = Cut-off
rvdw = 1.0
pbc = xyz
constraints = none
verlet-buffer-tolerance = -1
nstenergy = 1
nstcalcenergy = 1
integrator = nm
nsteps = 0

gmx_d grompp -f nm.mdp -c minimum.g96 -p topol.top -o nm.tpr
gmx_d mdrun -deffnm nm -mtx hessian.mtx -ntmpi 1 -ntomp 2 -nb cpu -bonded cpu -update cpu
gmx_d nmeig -f hessian.mtx -s nm.tpr -first 7 -last 9 -of frequency.xvg -ol eigenvalue.xvg -v eigenvector.trr -T 300 -qc quantum-correction.xvg
gmx_d dump -mtx hessian.mtx > hessian.txt

Hessian 是 9×9 矩阵。有限、非线性的三原子分子有三个平移和三个转动近零模，其余三个为内部振动；所以本例取第 7–9 模。这个编号不能照搬到周期晶体、线性分子或带约束的大分子。

4. 看实际结果#

最小化最大力约 4.03e-7 kJ mol⁻¹ nm⁻¹。独立读取 Hessian、按质量加权并对角化得到相同频率，最大差别约 0.0035 cm⁻¹，主要来自 XVG 的输出位数。

cat frequency.xvg
cat quantum-correction.xvg

300 K 时，三个内部谐振子的焓量子修正为 50.8939 kJ/mol，热容修正为 −24.8965 J/(mol·K)。负值是相对于经典内部谐振热容的修正，不代表体系热容为负；这些数也不是完整分子的标准生成焓或完整自由能。

图中三条流程同时涉及最小化、坐标保存和数值精度，不能将它们视为严格只改变一个变量的实验。完整第一组近零模的最大绝对频率约 0.249 cm⁻¹。

检查虚频时还要看原始特征值。 当前 nmeig 的频率输出会将负特征值对应频率写为零。请补充输出全部特征值，不能只看频率文件没有负号就断言稳定：

gmx_d nmeig -f hessian.mtx -s nm.tpr -first 1 -last 9 -of frequency-all.xvg -ol eigenvalue-all.xvg -v eigenvector-all.trr
cat eigenvalue-all.xvg

接下来可按结构集合与本征动力学监测生成动画和投影。

参数释义可查 GROMACS nmeig 官方说明。