GROMACS轨迹筛选与水环境分析

下载小体系输入、真实结果与绘图数据，在自己的目录解压：

unzip small-examples-*.zip
cd small-examples-20261027

1 准备轨迹与索引#

使用GROMACS 2026.4分析程序，输入为884个SPC/E水分子的100 ps真实轨迹；多数示例取前20 ps。CPU即可做这些分析。索引必须每行适当换行，不能把数千个编号拼成一行。

module load gromacs/2026.4-p100-threadmpi
mkdir -p work/water-tools
cp gromacs/shared/md.* work/water-tools/
cp gromacs/sorient/index.ndx work/water-tools/
cd work/water-tools

索引中0为第一个水分子的氧，1为其他883个完整水分子，2为全部水。这里研究指定水分子周围的环境，不将它称为蛋白溶剂化层。

2 轨迹频率筛选#

gmx filter -f md.xtc -s md.tpr -n index.ndx -nf 5 -ol lowpass.xtc -oh highpass.xtc -e 20

交互选择依次输入：2。

低通/高通轨迹适合观察不同时间尺度的运动。筛选改变坐标统计，不能将处理后的轨迹直接替代原轨迹计算热力学量。

3 旧版氢键分析#

gmx hbond-legacy -f md.xtc -s md.tpr -num hbonds.xvg -dist distance.xvg -ang angle.xvg -nthreads 1 -e 20 -nomerge

交互选择依次输入：0 → 0。

默认距离阈值0.35 nm、角度阈值30°；输出角度分布时显式使用-nomerge。新研究优先看新版hbond；旧版适合复现历史处理流程。

4 向盒子插入分子#

先复制下载包 gromacs/insert-molecules/water.gro 到当前工作目录：

cp ../../gromacs/insert-molecules/water.gro .

插入后为90原子，即30个水分子；还需更新拓扑中的分子数并做能量最小化，不能直接作为平衡水。

gmx insert-molecules -ci water.gro -box 2 2 2 -nmol 30 -seed 27171 -o inserted.gro

-seed固定随机插入种子，-box单位nm。这里是建模示例，不是液体密度计算。

5 水分子取向分布#

gmx sorient -f md.xtc -s md.tpr -n index.ndx -rmin 0.2 -rmax 0.5 -e 20 -o orientation.xvg -no normal.xvg -ro radial.xvg -co cumulative.xvg -rc counts.xvg

交互选择依次输入：0 → 1。

图中为0.2–0.5 nm壳层内水分子的取向分布，均匀取向的相对分布基线为1。

6 偶极壳层统计#

gmx spol -f md.xtc -s md.tpr -n index.ndx -rmin 0.2 -rmax 0.8 -dip 2.35 -e 20 -o polarization.xvg

交互选择依次输入：0 → 1。

SPC/E是固定电荷模型，-dip 2.35给定参考偶极矩（D）；不能把这个处理当成极化力场的诱导偶极计算。

7 近邻水排序#

gmx trjorder -f md.xtc -s md.tpr -n index.ndx -na 3 -da 1 -r 0.35 -e 20 -o ordered.xtc -nshell shell.xvg

交互选择依次输入：0 → 1。

输出按离参考氧的距离重排分子。帧与帧之间的原子编号不再对应同一个水分子，不能用重排后轨迹计算单分子扩散。壳层计数已与独立select查询逐帧对照。

8 介电谱处理流程#

printf '0\n' | gmx dipoles -f md.xtc -s md.tpr -e 100 -o dipole.xvg -eps epsilon.xvg -a average.xvg -d distribution.xvg -corr total -c dipole-acf.xvg
gmx dielectric -f dipole-acf.xvg -ffn none -eint 10 -tail 10 -bfit 1 -efit 10 -eps0 70 -epsRF 0 -o dielectric.xvg -d derivative.xvg -c cole.xvg

先从轨迹得到总偶极自相关。这里-eps0 70是明确给定的归一化参考值，不是这条短轨迹计算出的收敛静态介电常数；-eint 10截取短相关窗口，仅演示频域处理，不能据此报告材料介电谱。

grep -v '^[#@]' hbonds.xvg | head
grep -v '^[#@]' shell.xvg | head

下载包中各工具的 reference-*.xvg 为实际输出，便于逐列比较。增大体系前先检查选择组、PBC、单位与采样间隔；统计性质还需更长采样和独立重复。

重新绘制下载包中的数据图#

回到解压包根目录执行，图片输出到新的 replot 目录：

module load ase/3.29.0
python plot.py

脚本使用NumPy和Matplotlib，参考原始数值在 results.json，轨迹和工具输出在相应软件目录。