GROMACS 水层密度、界面与四面体结构

下载输入、真实结果与绘图资料包，解压进入资料包目录。inputs 可修改，reference 是集群已计算的参考结果；脚本每次在 work 下建立独立目录。计算在分配好的资源内执行，可用 sbatch run.sbatch 跑整组，结果路径写入 LAST-RUN.txt。

本例模型#

从 216 个 SPC/E 水分子的小盒出发，延长 Z 方向，在两端加入远距 12–6 壁，采用 pbc=xy 和 3dc 长程静电修正，运行 20 ps。盒约 1.863×1.863×7.450 nm³。它是小面积、短采样的水层操作例，包含壁边界，不是已收敛的自由液面或表面张力研究。

bash run.sh md

1. 查看边界和 MD 参数#

完整 inputs/slab/slab.mdp：

integrator = md
dt = 0.001
nsteps = 20000
continuation = yes
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 0.8
rvdw = 0.8
pbc = xy
nwall = 2
wall-type = 12-6
wall-atomtype = opls_116 opls_116
ewald-geometry = 3dc
tcoupl = v-rescale
tc-grps = System
tau-t = 0.1
ref-t = 300
pcoupl = no
gen-vel = no
constraints = none
nstenergy = 50
nstxout-compressed = 100
fourierspacing = 0.12

ld-seed = 1820309415

nwall=2 与 pbc=xy 配套；壁位于 z=0 和 z=Lz。3dc 是 3D reciprocal sum 加 Z 方向修正，不是完全不依赖真空厚度的精确二维静电解。壁距、真空层和面内尺寸仍需检查，不能假定壁贡献严格为零。

在独立目录复制 inputs/slab 后运行：

gmx grompp -f slab.mdp -c start.gro -p topol.top -o slab.tpr
gmx mdrun -deffnm slab -ntmpi 1 -ntomp 2 -nb cpu -bonded cpu -update cpu
printf '0\n' | gmx density -s slab.tpr -f slab.xtc -d Z -sl 100 -o mass-density.xvg

2. 按密度提取界面#

gmx select -s slab.tpr -select all -on system.ndx
printf '0\n' | gmx densorder -s slab.tpr -f slab.xtc -n system.ndx -dt 1 -bw 0.93129 -bwn 0.1 -d1 0 -d2 1000 -tblock 5 -og density-lower.xpm density-upper.xpm -or density-lower.out density-upper.out

选取全部原子、使用质量密度。-d1/-d2 给两相参考密度，单位 kg/m³；这里 0 与 1000 是演示阈值，应根据自己的实际平台密度调整。-bw 控制切向网格，-bwn 控制法向网格。小盒上的粗网格不能提供精细界面涨落谱。

本例每 1 ps 取一帧，5 帧平均一块，产生四个完整时间块，分别对应 0–4、5–9、10–14、15–19 ps。原生旧格式 raw 文件的时间字段不要直接作为这些块的物理时间，优先按实际抽帧与块编号映射；20 ps 的最后单帧不组成新的完整块。

3. 用 hydorder 看四面体结构阈值#

gmx select -s slab.tpr -select 'name OW' -on oxygen.ndx
printf '0\n' | gmx hydorder -s slab.tpr -f slab.xtc -n oxygen.ndx -dt 1 -bw 0.35 -sgang1 0 -sgang2 0.09146862 -tblock 5 -o sg-fine-lower.xpm sg-fine-upper.xpm -or sg-fine-lower.out sg-fine-upper.out

这里先独立读实际水层坐标，在中心液层测得 Chau 角序参量约 0.09147，再用它设定对比阈值。定义为 Sg=3/32×Σ六对(cosθ+1/3)²，理想四面体为零；不要与常见“理想四面体 q=1”的另一种定义混淆。

真空没有四面体结构。 工具将空格点记零只是一种数据约定，不能说真空比水更接近完美四面体。先用密度识别有物质区域，再解释其中的 Sg；hydorder 的阈值界面与密度界面也不必重合。

密度积分与 216 个水分子的总质量一致；图中 Sg 曲线只统计有原子的区域。侧面小、时间短、阈值和网格都有限，结果用于理解操作，不直接报告收敛界面厚度、毛细波或表面张力。

参考 MDP 边界与壁参数、densorder、hydorder。