GROMACS键角、取向与空间分析实操

先下载本批输入、参考数据与绘图脚本，解压进入 full-function-wave1。参考数据在 cases、tools，自己的新计算在 work，避免覆盖或向旧输出追加数据。

unzip full-function-wave1-*.zip
cd full-function-wave1
python start.py

1 用已有轨迹准备索引#

本篇使用GROMACS 2026.4分析程序和此前SPC/E水的100 ps真实轨迹。模块、编译器及GPU分支见环境指南；以下分析命令在CPU上执行，不要求分配GPU。先在自己的工作目录准备输入：

mkdir analysis-work
cd analysis-work
cp ../shared/gromacs/md.xtc .
cp ../shared/gromacs/md.tpr .
cp ../tools/gromacs-angle/angles.ndx .
cp ../tools/gromacs-rotacf/pairs.ndx .
cp ../tools/gromacs-densmap/index.ndx .
module load gromacs/2026.4-p100-threadmpi

angles.ndx 每三个编号是一组H–O–H，中心原子必须是O；pairs.ndx 每两个编号是一根O–H方向矢量。index.ndx 的第0组为884个水氧，第1组为前100个水氧。换体系时重新生成编号，不能照抄本例的原子顺序。

2 键角分布和取向相关#

printf '0
' | gmx angle -f md.xtc -n angles.ndx \
  -type angle -od angle-distribution.xvg -ov angle-time.xvg
printf '0
' | gmx rotacf -f md.xtc -s md.tpr \
  -n pairs.ndx -d -P 2 -o rotacf.xvg

本例平均H–O–H角为 104.518°，与SPC/E的刚性角104.52°相符。压缩轨迹坐标精度会使测得的角度出现小幅分布，不能把它解释成模型中真实的键角振动。rotacf 的 -d 使用成对原子，-P 2 求二阶Legendre相关；它描述取向记忆，与平移MSD不同。

3 二维数密度、距离偏离与回转半径#

printf '0
' | gmx densmap -f md.xtc -s md.tpr \
  -n index.ndx -bin 0.1 -od density.dat -o density.xpm
printf '1
' | gmx rmsdist -f md.xtc -s md.tpr \
  -n index.ndx -dt 1 -o rmsdist.xvg
printf '1
' | gmx gyrate-legacy -f md.xtc -s md.tpr \
  -n index.ndx -o gyrate.xvg

密度图平均掉z方向，横纵坐标为nm，默认归一化是数密度。这里选氧原子，每个水分子贡献一个中心；不要把它当作包含所有原子的质量密度。rmsdist 比较原子对距离的变化，不做普通RMSD的刚体拟合。前100个水氧构成的是演示用子集，回转半径不代表水分子的尺寸；分析蛋白或聚合物时，应换成完整的对应分子组。

4 距离连通簇与势函数检查#

printf '0
' | gmx clustsize -f md.xtc -s md.tpr \
  -n index.ndx -cut 0.28 -dt 5 \
  -nc count.xvg -mc largest.xvg -ac average.xvg -o clusters.xpm
gmx sigeps -sig 0.316557 -eps 0.650194 -o lj.xvg

前者按0.28 nm氧—氧距离建立连通关系，输出簇数量和最大簇；它不同于按角度与距离共同定义的氢键网络。后者使用SPC/E氧的LJ参数绘制两体势，σ单位nm，ε单位kJ/mol；不包含整套水模型的静电贡献。

5 单独估计PME误差#

gmx pme_error 要求用于估计的输入满足特定邻居表条件。附件保留了可复现设置；生成单独的估计输入，不修改此前已经完成的模拟文件。

cp ../shared/gromacs/md.mdp pme-check.mdp
cat >> pme-check.mdp <<'EOF'
verlet-buffer-tolerance = -1
rlist = 1.0
EOF
gmx grompp -f pme-check.mdp -c ../shared/gromacs/md.gro \
  -t ../shared/gromacs/md.cpt -p ../shared/gromacs/topol.top -o pme-check.tpr
gmx pme_error -s pme-check.tpr -o pme-error.out -self 0.1 -seed 27131

本例实空间估计为0.05533 kJ/(mol·nm)，倒空间估计为1.854 kJ/(mol·nm)。这是特定网格、插值阶数及抽样条件下的工具估计，不是材料的误差条；正式计算应结合网格、截断和实际力的对照检查。

可查看GROMACS命令手册。附件包含上述全部数值文件，用 python plot.py 重画。