GROMACS水取向、电流与空间分布实操

先下载完整输入、操作脚本与实际结果并解压。gromacs/reference和各LAMMPS输入目录中的矩阵为参考结果；请在新建工作目录复现。

1 使用现有轨迹做分析#

环境为GROMACS 2026.4 P100 thread-MPI构建，分析命令在CPU执行，不需要额外启动MD。水输入为884个SPC/E水分子的已完成轨迹；电流分析使用此前的盐水实算。每个脚本读取前20 ps，分别建立自己的工作目录。

cd gromacs
module purge
bash run-h2order.sh
bash run-current.sh
bash run-spatial.sh

这是短数据处理示例，不是收敛的取向、介电或电导率研究。

2 沿盒子方向计算水偶极取向#

gmx h2order -s md.tpr -f md.xtc -n index.ndx \
  -e 20 -d Z -sl 10 -o orientation.xvg

交互选择0，即完整水分子组。索引必须按每个水 O–H–H 的顺序组织。-sl 10沿z方向划分10层，每个水按氧位置归入层中。输出为水偶极与z方向夹角余弦的平均值，范围应在−1至1；它不同于二阶Legendre取向序参量。

本例是均匀水盒，不能把随机层间波动解释成界面取向。输出见work/h2order/orientation.xvg。

3 由含速度轨迹计算电流相关#

gmx current -s md.tpr -f md.trr -e 20 -sh 1 \
  -bfit 1 -efit 10 -bvit .1 -evit 1 \
  -o current.xvg -caf caf.xvg -dsp displacement.xvg \
  -md dipole.xvg -mj translation.xvg

选择0，即完整体系。md.trr必须包含速度，不能把仅含压缩坐标的XTC当作速度轨迹。这里电流输出单位为e·nm/ps，caf.xvg为相关函数及相关处理列。

本例只有21个速度采样帧；相关曲线尾部和拟合高度依赖时间窗口，不据此报告材料电导率。正式研究需要适合的电荷模型、足够密的速度输出、更长采样、统计误差及有限尺寸检查。

4 空间分布图先确认坐标系与网格范围#

gmx spatial -s md.tpr -f md.xtc -n index.ndx \
  -e 20 -bin .2 -nab 32 -nodiv

两次都选择0，即水氧组。-bin .2为0.2 nm网格，-nodiv输出每帧占据计数。本例不做蛋白拟合，得到的是实验室坐标系下的水氧占据分布；溶剂化壳层分析需要先围绕目标溶质做合适的居中和拟合。

在当前2026.4版本中，导出的grid.cube会省略最高端部分占据网格。本例导出区域积分为 873.651523个氧/帧，不是总数884；不要直接据此做全体系粒子数或质量守恒统计。

需要完整计数时，可从同一轨迹提取坐标，再用包内独立脚本做完整网格直方图：

cp complete_grid.py work/spatial/
cd work/spatial
printf "0\n" | gmx traj -s md.tpr -f md.xtc -n index.ndx -e 20 -ox positions.xvg
module purge
module load ase/3.29.0
python complete_grid.py

本次完整统计为 884.000000个氧/帧。脚本网格适用于当前约3 nm水盒；换体系后必须扩展边界，并重新检查是否包含全部坐标。对网格内位置的处理不能代替正确的周期边界和分子重建。

图右下为独立完整网格统计；左下的蛋白链回转半径见蛋白链与散射工具。

参考：GROMACS 2026.4命令手册、spatial输出实现。