GROMACS flooding 偏置与恒力采样
以真实轨迹比较高斯、反向高斯与恒力，说明强度、作用区域、实际偏置能与投影力的终端提取。
阅读本篇目录
本篇在ED约束教程相同的三原子模型与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。
module load gromacs/2026.4-cpu-double
bash prepare.sh
sbatch job-p100.slurm
任务顺序运行28个小条件，包含两个种子、五类约束和三种偏置。每条轨迹40 fs、201帧。CPU双精度版本、完整拓扑与MDP见约束教程。
2 固定强度高斯#
gmx_d make_edi -f pca-vectors.trr -eig pca-eigenvalues.xvg -s minimum.g96 -flood 1-3 -tau 0 -Eflnull 2.5 -T 300 -outfrq 1 -o sampling.edi
-tau 0让强度保持常数，-Eflnull指定强度。本例比较0.02与2.5 kJ/mol；它们是教学模型对照值，实际蛋白应依据自身基底与采样分布选择。
本例高斯形式为：
V(q)=Efl·exp[-KT/(2·Efl·alpha²)·Σ wi·(qi−qi,ref)²]
wi是PCA特征值的倒数，由原生程序写入EDI；默认中心是PCA平均结构。这里不使用质量加权，投影长度为nm。查看实际参数与向量：
head -35 sampling.edi
rg -A 15 'COMPONENTS GROUP 7' sampling.edi
弱强度同时让作用区域变窄，离开中心后可能几乎不起作用，应该看实际势和力，而非仅看输入数值。
3 反向高斯和恒力#
反向高斯：
gmx_d make_edi -f pca-vectors.trr -eig pca-eigenvalues.xvg -s minimum.g96 -flood 1-3 -tau 0 -Eflnull 2.5 -T 300 -restrain -outfrq 1 -o sampling.edi
-restrain反转偏置方向，程序按约定处理内部强度和alpha²符号。直接使用这个选项，不单独改变高斯指数符号。附件保留正向/反向两档强度和两个种子的全部数据。
恒力：
gmx_d make_edi -f pca-vectors.trr -eig pca-eigenvalues.xvg -s minimum.g96 -flood 1 -constF "0.02" -outfrq 1 -o sampling.edi
沿PC1直接施加0.02 kJ/(mol·nm)子空间力。原生Vfl列在恒力模式为0，力列仍为0.02；Vfl=0不能解释为没有偏置。
4 实际运行并提取#
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 -50 sampling.xvg
gmx_d energy -f md.edr -o energy.xvg
# 选择Potential、Kinetic-En.、Total-Energy、Temperature，最后0
XVG中每个EV依次给出投影、Vfl和该方向的力。高斯Vfl是全局偏置势，在多个EV旁重复显示，不要把相同的Vfl列相加。常规能量文件与偏置输出按各自定义读取。
| 条件 | 种子 | 最大偏置能绝对值/(kJ/mol) | 最大子空间力绝对值/(kJ/mol/nm) |
|---|---|---|---|
| flood-constant-s13579 | 13579 | 0 | 0.02 |
| flood-constant-s24680 | 24680 | 0 | 0.02 |
| flood-gaussian-E2p5-s13579 | 13579 | 1.79794 | 5081.47 |
| flood-gaussian-E2p5-s24680 | 24680 | 2.0082 | 6502.88 |
| flood-gaussian-s13579 | 13579 | 2.09186e-06 | 0.0748428 |
| flood-gaussian-s24680 | 24680 | 1.26086e-08 | 0.00273756 |
| flood-restrain-E2p5-s13579 | 13579 | 2.48164 | 6510.38 |
| flood-restrain-E2p5-s24680 | 24680 | 1.84979 | 5992.15 |
| flood-restrain-s13579 | 13579 | 2.09186e-06 | 0.0748429 |
| flood-restrain-s24680 | 24680 | 1.26086e-08 | 0.00273756 |
最大子空间力取所有所选方向与时刻的绝对值最大值，不是矢量模，也不是某个原子的总受力。两档强度和全部种子都保留。
5 分析约定与新结果绘图#
- flooding输出与同一时刻TRR配对；ED约束/监测则按更新后的步序配对，附件分别处理。
- 原生q、V、F列保留六位有效数字，高斯会放大投影舍入误差。脚本按实际末位区间传播误差，并独立核对坐标投影。
- 本例按300 K生成初速后没有恒温器；偏置改变运动，40 fs不足以形成平衡分布。这些结果用于学习设置和提取，不能直接当作自由能面。
重画新任务：
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
head -8 analysis-new/flood-gaussian-E2p5-s13579/flooding.csv
三种偏置均采用同一组非质量加权PCA。自适应强度、谐波或移动中心另有参数含义，不能只改一个开关后沿用本篇能量解释。参数见GROMACS make_edi说明。