GROMACS 本征动力学约束与定向采样

本征动力学监测只记录轨迹投影。本篇实际施加约束，让集体坐标固定、按步长移动、只接受某方向的运动，或朝目标收缩。

本例是三原子的柔性TIP3P类型教学水分子，适合快速掌握参数和输出。PCA来自真实40 fs轨迹，随后每个条件再运行40 fs；它不是收敛的蛋白质PCA，也不能用于比较蛋白构象采样效率。

1 下载并重画#

下载本征动力学约束与flooding完整实例包，约2.7 MiB，含28个实际计算条件的完整TRR、TPR、EDR、EDI与绘图脚本。

unzip gromacs-ED-constraints-flooding-wave233*.zip -d ED-example
cd ED-example
sha256sum -c MANIFEST.sha256
bash view.sh
ls analysis/figures

view.sh重画参考图和CSV，不启动新模拟。分析依赖MDAnalysis2.9.0、NumPy和Matplotlib；可用 ANALYSIS_PYTHON指定自己的Python。输入在 inputs，参考结果在 reference-results，新计算在 cases。

为比较键长、键角，图中分别将O原子置于原点、第一条O–H转到横轴。这是显示坐标变换，原始轨迹保留完整坐标。

2 选择双精度CPU入口#

module purge
module load gromacs/2026.4-cpu-double
gmx_d --version

程序使用固定v2026.4源码构建，显示 2026.4-dev，GNU13.3、thread-MPI/OpenMP、内部FFTPACK/BLAS/LAPACK。本例在CPU上运行。P100、V10016G模板在相应节点分配一个CPU核心；GPU队列要求同时申请匹配GPU，模板按此保留申请，GPU不参与计算。这些任务用于操作示范，不用于比较GPU速度；已有可用CPU任务资源时可用同一双精度入口。

3 从实际轨迹建立PCA方向#

python3 make_basis.py
ls basis-new
cat basis-new/pca-eigenvalues.xvg

脚本读取201帧 source-unbiased40fs.trr，生成带显式物理质量的分析TPR，再调用非质量加权PCA。核心命令在 basis-new目录中执行：

gmx_d grompp -f nm.mdp -c minimum.g96 -p topol.top -o basis.tpr
gmx_d covar -s basis.tpr -f source-unbiased40fs.trr -nomwa -o pca-eigenvalues.xvg -v pca-vectors.trr -av pca-average.pdb -ascii pca-covariance.dat
# 本例两次选System，即0、0

这里没有重新计算Hessian。-nomwa表示非质量加权，拟合与后续ED保持同一约定。本例System只有3个原子，换成蛋白质应选择匹配的拟合组和分析组，不把含溶剂System直接套给蛋白主成分。

本篇用PCA的PC1–PC3；正常模教程的水分子内部Hessian模是7–9号，两套编号与向量约定不同。本例ED约束限定采用这里的非质量加权PCA。Hessian按原正常模与监测教程使用，不直接替换进本篇约束输入。附件按TRR的step读取模式编号，time不能当作模式编号。

4 准备并提交#

bash prepare.sh
ls cases
sbatch job-p100.slurm
squeue -u "$USER"

也可选择 job-v100.slurm。一个任务顺序运行28个小条件，使用初速度种子13579、24680，包括监测、五类约束、方向/零步长对照和偏置流程。每条轨迹200步，步长0.0002 ps，即 0.2 fs、总计40 fs、201帧。

完整拓扑如下，单位为nm、kJ/mol和u：

; Explicit flexible TIP3P-type teaching monomer; nm, kJ/mol, e units.
[ defaults ]
1 2 yes 0.5 0.8333333333
[ atomtypes ]
OW 8 15.9994 0 A 0.315061 0.636386
HW 1 1.008 0 A 0 0
[ moleculetype ]
SOL 2
[ atoms ]
1 OW 1 SOL OW 1 -0.834 15.9994
2 HW 1 SOL HW1 1 0.417 1.008
3 HW 1 SOL HW2 1 0.417 1.008
[ bonds ]
1 2 1 0.09572 502416
1 3 1 0.09572 502416
[ angles ]
2 1 3 1 104.52 628.02
[ system ]
Flexible water monomer
[ molecules ]
SOL 1

实际MDP：

cutoff-scheme = Verlet
nstlist = 1
rlist = 1.0
coulombtype = Cut-off
rcoulomb = 1.0
vdwtype = Cut-off
rvdw = 1.0
pbc = xyz
constraints = none
verlet-buffer-tolerance = -1
nstenergy = 1
nstcalcenergy = 1
nstlog = 20
integrator = md
dt = 0.0002
nsteps = 200
nstxout = 1
nstvout = 1
comm-mode = Linear
nstcomm = 1
gen-vel = yes
gen-temp = 300
gen-seed = 13579

按300 K生成初速后，本例没有恒温器。柔性O–H振动使用0.2 fs步长；ED改变运动，不能把它解释为无偏平衡采样。md-processed.mdp保存程序实际读取的设置。

5 五类约束怎样设置#

下列命令使用同一组PCA文件与匹配结构，在准备好的案例目录执行。自动任务完成后，另建目录进行手动练习。按提示选择本例System，即0。

LINFIX：固定增量或固定投影#

gmx_d make_edi -f pca-vectors.trr -eig pca-eigenvalues.xvg -s minimum.g96 -linfix 1 -linstep "0.00001" -outfrq 1 -o sampling.edi

-linstep单位为nm/步。本例PC1每步增加0.00001 nm，200步增加0.002 nm；负值沿反方向移动，0则固定该投影，其他自由度仍能运动。附件有零、正、负三组。

LINACC：只接受某方向#

gmx_d make_edi -f pca-vectors.trr -eig pca-eigenvalues.xvg -s minimum.g96 -linacc 1 -accdir "1" -outfrq 1 -o sampling.edi

-accdir只使用正负号。正号接受投影增加的运动，负号相反；它没有固定增量，不应期待匀速曲线。附件保留两个方向。

RADFIX：固定半径增量#

gmx_d make_edi -f pca-vectors.trr -eig pca-eigenvalues.xvg -s minimum.g96 -radfix 1-3 -radstep 0.00001 -ori target.g96 -maxedsteps 0 -outfrq 1 -o sampling.edi

半径是PC1–PC3子空间到指定原点的距离，不是O–H键长。这里显式指定最小化结构为扩张原点，初始结构有小位移；-maxedsteps 0关闭周期性重设原点，从非零初始半径开始扩张。

RADACC：只接受半径增加#

gmx_d make_edi -f pca-vectors.trr -eig pca-eigenvalues.xvg -s minimum.g96 -radacc 1-3 -maxedsteps 1000000 -slope 0 -outfrq 1 -o sampling.edi

接受远离起点的运动，拒绝半径减小，增量不是常数。附件按原生步序读取首个扩张起点，与同种子、同输入的监测组配对；改动构型、种子或步长时，对照组同步生成。真实问题也可用 -ori明确指定原点。

RADCON：朝目标收缩#

gmx_d make_edi -f pca-vectors.trr -eig pca-eigenvalues.xvg -s minimum.g96 -radcon 1-3 -tar target.g96 -outfrq 1 -o sampling.edi

目标是最小化结构，接受向目标子空间位置接近的运动。目标文件须与向量的原子顺序、数量及拟合定义匹配。三类半径可能采用不同原点，应先看定义再比较数值。

半径来自保存的TRR坐标，各自使用实际扩张原点或收缩目标。

6 实际运行、读取与绘图#

gmx_d grompp -f md.mdp -c initial.g96 -p topol.top -o md.tpr
gmx_d mdrun -deffnm md -ei sampling.edi -eo sampling.xvg -ntmpi 1 -ntomp 1 -nb cpu -bonded cpu -update cpu
head -45 sampling.xvg
head -30 sampling.edi

每个向量在同一份EDI里只选一次，监测和约束方向不能随意重复填写。EDI记录模式、步长和原点，XVG记录实际投影与半径。

ED约束/监测输出发生在积分更新后，对应下一帧TRR坐标。脚本按步序比较200对数据；最后一条ED输出没有另存的更新后坐标，不需要用户手动平移曲线。

处理新任务时，先为绘图解释器取短变量名，再运行：

analysis_py=/public/toolkit/miniconda3/envs/qvasp-2.26/bin/python
"$analysis_py" analyze.py --cases cases --output analysis-new
"$analysis_py" plot_geometry.py --cases cases --output analysis-new/figures
ls analysis-new/figures

依赖与 view.sh相同，输出逐帧投影、半径、键长和键角CSV及图片。继续偏置与恒力教程。参数见GROMACS make_edi说明。