GROMACS 压力相关与黏度实操

下载完整水盒输入、真实轨迹、压力数据与绘图脚本。解压后进入资料包目录，先查看已有结果：

bash analyze.sh

会在 reference/figures 生成4张高清图，在 reference/analysis 生成汇总CSV。脚本调用集群既有Python3.11后处理环境，包含MDAnalysis2.9.0、NumPy1.26.4和SciPy1.13.1，无需自己找绘图依赖。

重新计算整组，在新解压目录执行：

sbatch submit.slurm
squeue -u "$USER"

数组模板为16个小任务，每个占1张P100和2个CPU线程。核时用户使用模板中的P100核时分区；套餐用户按自己的账户规则修改分区。等全部任务结束后：

bash analyze.sh work

已有 work 的目录不会被覆盖。reference 是已有参考，work 是你新算的结果；两者不要混在一起。

环境为GROMACS2026.4 P100 thread-MPI构建，GNU13.3、CUDA12.8、FFTW3.3.10；动力学的短程非键和PME使用GPU，陡降优化的PME使用CPU。本例不涉及PLUMED。水采用SPC/E固定电荷模型，拓扑内部的SETTLE保持刚性；constraints = none 不会关闭拓扑里定义的SETTLE。

1 从建模到生产，先做哪些事#

本例顺序是：建立纯水盒 → 能量最小化 → 20 ps NPT → 20 ps NVT → 固定盒子的60 ps NVE → 提取每步压力 → 积分与画图。4种初始盒长2.2、2.6、3.0、3.4 nm，每种4个速度种子，分别得到348、578、884、1284个水分子。16项是采样对照，不是16种不同性质。

下面手动操作要在已分配的GPU资源中执行。先从附件复制一个空白输入目录：

cp -r inputs/water-L2p2-s290112 water
cd water
module purge
module load gromacs/2026.4-p100-threadmpi
gmx solvate -cs spc216.gro -box 2.2 2.2 2.2 -o solvated.gro -p topol.top
cat topol.top

topol.top 引用 oplsaa.ff/forcefield.itp 和 oplsaa.ff/spce.itp；本例是纯水，没有额外有机分子的OPLS参数。gmx solvate 会在分子表中写入实际水数量。

2 完整输入怎么写#

能量最小化 em.mdp#

integrator = steep
nsteps = 3000
emtol = 1000
emstep = 0.01
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 0.9
rvdw = 0.9
DispCorr = EnerPres
fourierspacing = 0.12
pme-order = 4
pbc = xyz
constraints = none
nstlist = 20
verlet-buffer-tolerance = 0.005

NPT预平衡 npt.mdp#

integrator = md
dt = 0.002
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 0.9
rvdw = 0.9
DispCorr = EnerPres
fourierspacing = 0.12
pme-order = 4
pbc = xyz
constraints = none
nstlist = 20
verlet-buffer-tolerance = 0.005
nstenergy = 100
nstcalcenergy = 10
nstlog = 1000
nstxout = 0
nstvout = 0
nstfout = 0
nstxout-compressed = 0
comm-mode = Linear
nstcomm = 100
nsteps = 10000
continuation = no
tcoupl = v-rescale
tc-grps = System
tau-t = 0.2
ref-t = 300
pcoupl = C-rescale
pcoupltype = isotropic
tau-p = 2.0
ref-p = 1.0
compressibility = 4.5e-5
gen-vel = yes
gen-temp = 300
gen-seed = 290112
ld-seed = 290212

固定盒子热化 nvt.mdp#

integrator = md
dt = 0.002
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 0.9
rvdw = 0.9
DispCorr = EnerPres
fourierspacing = 0.12
pme-order = 4
pbc = xyz
constraints = none
nstlist = 20
verlet-buffer-tolerance = 0.005
nstenergy = 100
nstcalcenergy = 10
nstlog = 1000
nstxout = 0
nstvout = 0
nstfout = 0
nstxout-compressed = 0
comm-mode = Linear
nstcomm = 100
nsteps = 10000
continuation = yes
tcoupl = v-rescale
tc-grps = System
tau-t = 0.2
ref-t = 300
pcoupl = no
gen-vel = no
ld-seed = 290312

压力采样 md.mdp#

integrator = md
dt = 0.002
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 0.9
rvdw = 0.9
DispCorr = EnerPres
fourierspacing = 0.12
pme-order = 4
pbc = xyz
constraints = none
nstlist = 20
verlet-buffer-tolerance = 0.005
nstenergy = 1
nstcalcenergy = 1
nstlog = 1000
nstxout = 0
nstvout = 0
nstfout = 0
nstxout-compressed = 100
comm-mode = Linear
nstcomm = 100
nsteps = 30000
continuation = yes
tcoupl = no
pcoupl = no
gen-vel = no

nstcalcenergy=1 每2 fs计算能量/压力，nstenergy=1 每步写入EDR。压力相关法需要这样的密集采样；把普通稀疏能量文件交给黏度工具，不能补回已经丢掉的高频信息。XTC每0.2 ps写坐标，供扩散分析，不包含速度。

3 逐段运行#

gmx grompp -f em.mdp -c solvated.gro -p topol.top -o em.tpr
gmx mdrun -deffnm em -ntmpi 1 -ntomp 2 -nb gpu -pme cpu -bonded cpu -update cpu -pin off

gmx grompp -f npt.mdp -c em.gro -p topol.top -o npt.tpr
gmx mdrun -deffnm npt -ntmpi 1 -ntomp 2 -nb gpu -pme gpu -bonded cpu -update cpu -pin off

gmx grompp -f nvt.mdp -c npt.gro -t npt.cpt -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt -ntmpi 1 -ntomp 2 -nb gpu -pme gpu -bonded cpu -update cpu -pin off

gmx grompp -f md.mdp -c nvt.gro -t nvt.cpt -p topol.top -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 2 -nb gpu -pme gpu -bonded cpu -update cpu -pin off

-t把上一段的速度等状态传给下一段，NVE段不重新生成速度。仅仅跑完20 ps NPT不能证明体系已经在1 bar充分平衡；更换体系后应延长平衡并查看温度、密度、压力及能量漂移。

4 提取压力张量，再计算原生黏度#

printf 'Pres-XX\nPres-XY\nPres-XZ\nPres-YX\nPres-YY\nPres-YZ\nPres-ZX\nPres-ZY\nPres-ZZ\nTemperature\nPressure\n0\n' |
  gmx energy -s md.tpr -f md.edr -o pressure.xvg -dp

gmx energy -s md.tpr -f md.edr -o selected.xvg \
  -vis viscosity.xvg -evisco einstein.xvg -eviscoi einstein-integral.xvg \
  -acflen 3001 -einstein_blocks 4 -einstein_restarts 20

固定体积EDR可能没有Volume项，工具此时会询问盒体积，单位nm³。参考算例的实际立方边长2.1799 nm，输入体积 10.358806345399 nm³。自己重算后盒长可能不同，要用自己的最终盒子；自动运行脚本已从对应 md.gro 读取，不需要手工修改。原生 selected.xvg 在缺Volume项时可能出现重复的LJ列，不要把该列当作体积。

• viscosity.xvg：第一列ps，第二列剪切积分，第三列体积压力涨落的积分；原生cP等于mPa·s。
• einstein-integral.xvg：压力时间积分的均方增长，单位Pa·s·ps。
• einstein.xvg：上述量的有限差分，单位Pa·s；换为mPa·s乘1000。

Einstein法的能量积分可以利用EDR保存的步间和，但仍要足够频繁计算压力；不能据此认为稀疏计算压力也没有影响。

5 为什么曲线没有一个自动可信的最终值#

本例第一条348水轨迹的剪切积分，积分到1、3、6 ps分别约 0.741、1.125、1.277 mPa·s，还没有可靠平台。其短轨迹剪切压力均值约−55.0、70.4、−15.4 bar，也不是严格零。原生工具的剪切相关保留这些均值；图中虚线另做“减去本段样本均值”的有限时长对照，6 ps约0.857 mPa·s。虚线不是偷偷替换原生结果，也不因更平就自动更准确。体积压力相关则按原生方法减去平均压力。

原生Einstein的稀疏差分时间标签并非实际差分区间中点；附件保留原始文本，图上按实际区间中点画相同的导数值。不要拿这种标签偏移解释流体的时间尺度。

这里60 ps NVE只是普通用户的操作参考，尚不足以报告收敛黏度。更长采样、多个种子、合理相关窗口、体系尺寸、电荷处理和误差分析都要继续比较。大幅变化、负的后段积分或未出现平台时，应保留数据并延长统计，不要挑一个好看的数值当材料结论。

本例方法依据GROMACS黏度说明及energy参数。扩散和盒长比较继续看水盒扩散与有限尺寸。