GROMACS 横向电流相关与黏度工具

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

为什么普通 XTC 不够#

tcaf 同时需要每帧坐标与速度。XTC 没有速度；只有速度的 TRR 也不能完成该分析。本例使用 216 个 SPC/E 水分子、648 原子，10 ps NVE，每 1 fs 保存坐标和速度，共 10001 帧。

bash run.sh md

需要单独运行时，在新目录复制 inputs/velocity，加载普通 GROMACS 2026.4 模块。

1. 检查完整 MD 输入#

integrator = md
dt = 0.001
nsteps = 10000
continuation = yes
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 0.8
rvdw = 0.8
pbc = xyz
tcoupl = no
tc-grps = System
tau-t = 0.1
ref-t = 300
pcoupl = no
gen-vel = no
constraints = none
nstenergy = 10
nstxout-compressed = 100
fourierspacing = 0.12

nstxout = 1
nstvout = 1

ld-seed = -67143713

nstxout=1 和 nstvout=1 保证 TRR 每步同时含坐标与速度；高频输出会明显增加文件量。本例完整 TRR 约 150 MiB，研究大体系前先估算存储。

gmx grompp -f md.mdp -c start.gro -p topol.top -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 2 -nb cpu -bonded cpu -update cpu
gmx check -f md.trr

初始结构与速度来自同模型此前的短水盒计算，没有在 NVE 段重新随机生成速度。真实研究还应先核对温度、密度和能量漂移。

2. 运行横向流相关分析#

printf '0\n' | gmx tcaf -s md.tpr -f md.trr -ot transverse-current.xvg -oa all-current.xvg -o average-current.xvg -of fits.xvg -oc cubic-average.xvg -ov viscosity.xvg -acflen 500 -wt 0.02

本例选全部原子组 0。-acflen 指相关长度，-wt 是拟合权重衰减时间，单位 ps；这些都应结合体系时间尺度重新选择。-oc 适用于本例的立方盒，按相同波矢长度平均。

第一幅图展示质量加权横向流；工具旧的坐标标签可能仅写 nm/ps，但此原子算法实际累加了质量。独立用 TRR 坐标和速度计算首个波矢的正弦/余弦分量，与原生文本前 21 帧最大差约 4.8e-4，对应文本输出精度。

3. 怎样读黏度文件#

head -25 average-current.xvg
cat viscosity.xvg

viscosity.xvg 包含多个数据段，& 分隔个别波矢与按壳层平均的结果。横坐标是 k，单位 nm⁻¹；纵坐标的 10⁻³ kg/(m·s) 等于 mPa·s。本例壳层表观值约 0.21–0.23 mPa·s，只是这段小盒短轨迹与所选拟合窗口的输出。

它不是已经外推的宏观水黏度。 真正的 k→0 极限需要波矢/盒长、相关窗口、拟合和长采样检查。不要将几个有限 k 点的简单平均冒充宏观结果，也不要把异常的长相关时间拟合当作自动可信的统计误差。

方法与参数见 GROMACS tcaf。