GROMACS 非平衡余弦加速

GROMACS可以沿x方向施加随z位置变化的余弦加速度，生成周期性的非平衡驱动，用于黏度相关方法。本篇先完整跑通输入、提交和原生数据提取，再说明为什么一段短输出不能直接作为黏度。

下载完整输入、真实轨迹与绘图资料包，约147 MiB。资料包含32个短算例及6次同构型电场力对照。解压进入资料包目录，先查看已有结果：

bash view.sh

四张轨迹数据图位于 analysis/figures，电场力图位于 field-analysis/figures；每个条件的CSV和数值汇总位于 analysis/条件名。脚本调用集群已有MDAnalysis2.9.0、NumPy、Matplotlib环境，画已有数据时不提交计算。

本例使用348个SPC/E水分子，共1044个原子，初始坐标来自已有水盒的NVT平衡末帧。新计算重新生成速度，300 K、固定盒子NVT、步长1 fs、长度20 ps，每0.1 ps保存一帧，共201帧完整位置/速度TRR。两种速度种子13579、24680分别在P100和V100 16G计算。不同平台重复相同种子，不当作四个独立初始种子。

GROMACS2026.4 thread-MPI，GNU13.3、CUDA12.8、FFTW3.3.10。短程非键和PME使用GPU，键相互作用及坐标更新使用CPU；Colvars集体变量的计算不因此变成GPU内核。

1 准备与提交#

python3 prepare.py
cd cases/p100-cos-plus-s13579
cat md.mdp
sbatch job.slurm
squeue -u "$USER"

负向加速改进 p100-cos-minus-s13579；第二种子与V100目录也已准备。流程：同一水盒坐标 → 生成指定种子速度 → 施加位置相关加速度 → 20 ps NVT → 提取速度投影和原生逆黏度项 → 检查时间变化与重复结果。

模板每项申请1张对应GPU和1个CPU核心，使用核时分区；套餐用户按账户规则修改分区。不要在登录节点直接运行完整MD。已有 cases 或结果目录时，重新解压到新目录准备，避免覆盖。

任务实际执行的主命令如下，供已分配资源中的手动操作参考：

gmx grompp -f md.mdp -c initial.gro -p topol.top -o md.tpr   -po md-processed.mdp -pp md-processed.top
gmx mdrun -deffnm md -ntmpi 1 -ntomp 1 -nb gpu -pme gpu   -bonded cpu -update cpu

拓扑完整分子表如下；原子电荷、刚性水几何由所引用的SPC/E模型给出：

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/spce.itp"

[ system ]
SPC/E liquid transport tutorial in water

[ molecules ]
SOL               348

2 完整MDP和加速度的含义#

integrator = md
dt = .001
nsteps = 20000
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
rcoulomb = .9
rvdw = .9
DispCorr = EnerPres
fourierspacing = .12
pme-order = 4
pbc = xyz
constraints = none
nstenergy = 100
nstcalcenergy = 100
nstlog = 1000
nstxout = 100
nstvout = 100
nstfout = 0
nstxout-compressed = 0
tcoupl = v-rescale
tc-grps = System
tau-t = .5
ref-t = 300
pcoupl = no
comm-mode = Linear
nstcomm = 100
continuation = no
gen-vel = yes
gen-temp = 300
gen-seed = 13579
ld-seed = 13680
cos-acceleration = .02

关键项是：

a_x(z) = A * cos(2*pi*z/Lz)

cos-acceleration = .02 的A为0.02 nm/ps²，加速度沿x，沿盒子z方向变化。不是给每个原子同样大小的恒定x加速度。改成 -.02 会反转驱动相位。

恒温器、流动速度与温度的定义应遵循相应非平衡方法，不能把任意速度漂移都当作目标流场。自己更换体系时，检查模拟盒形状、加速度强度、周期边界与原生方法要求。

3 原生能量项的名称必须准确#

module load gromacs/2026.4-p100-threadmpi
printf '2CosZ*Vel-X
1/Viscosity
0
' | gmx energy   -f md.edr -o cosine-response.xvg

本版本两个能量项叫 2CosZ*Vel-X 和 1/Viscosity。不能凭习惯写成Vcos：原生程序可能提示某个名称不匹配，但因为另外一个项选中了，仍正常退出。提取完要看XVG的legend，确认两列都在。

2CosZ*Vel-X 是相应余弦速度投影，单位nm/ps；第二项是原生方法给出的逆黏度量，原生输出单位m·s/kg，即(Pa·s)⁻¹。不能把列名中的“Viscosity”直接当成已收敛、可引用的黏度数值。

4 任务结束后提取结果#

在当前算例目录执行：

module load gromacs/2026.4-p100-threadmpi
python3 ../../export_native.py
ls native-analysis

V100更换为对应模块。这个步骤读取已有轨迹和能量，输出到新建的 native-analysis，不会覆盖同名分析目录。全部参考结果及绘图脚本可在资料包根目录运行 bash view.sh 查看。

在当前条件目录画自己的新结果：

module load ase/3.29.0
python - <<'PY'
import re
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
p=Path("native-analysis/cosine-response.xvg")
labels={label:int(i)+1 for i,label in re.findall(r'@\s+s(\d+)\s+legend\s+"([^"]+)"',p.read_text())}
d=np.loadtxt(p,comments=["#","@"],ndmin=2)
plt.plot(d[:,0],d[:,labels["2CosZ*Vel-X"]])
plt.xlabel("Time (ps)");plt.ylabel("2CosZ*Vel-X (nm/ps)")
plt.tight_layout();plt.savefig("new-cosine-response.png",dpi=180)
PY

这段程序直接读取刚生成的原生数据；PNG保存在当前目录。XVG以 # 和 @ 开头的行是说明/列名，不能当作数值读入。

5 八组实际输出怎样读#

本篇读中间和右列；左列见外电场与偶极响应。全部保存时刻和重复条件保留，没有删除负值或只取看起来平稳的一段。

这些20 ps短轨迹中的速度投影仍被涨落明显影响，正向驱动的短窗均值也可能为负，反之亦然。原生逆黏度值会随时间改变甚至变号，本例不从这些数值报告材料黏度。

正确的用途是先确认任务和分析流程，再延长平衡与生产，检查受驱速度分布是否形成稳定响应，比较不同加速度、盒长、时间窗口及独立种子。不能对每个瞬时逆黏度取倒数再平均，也不能删掉负值后只平均正值。

若目标是平衡态压力相关方法，可继续GROMACS压力相关与黏度实操。两种方法的控制变量、响应和采样要求不同，需要分别检查，不能直接拼接输出。

参数说明见GROMACS非平衡MD选项。