GROMACS PCA、聚类与接触图
PCA、聚类与接触图，包含终端操作、输入文件、实际数据和适用范围。
阅读本篇目录
本文以鸢算科技集群 GROMACS 2026.4 实际结果为例。下载输入、参考轨迹、分析与绘图资料包，资料包约 529 MB，包含参考轨迹。解压后进入 gromacs-examples。参考输出在 cases；新计算由 new-case.py 建立独立目录。
本篇任务与操作要点#
PCA 前先拟合去掉整体平移和转动。本例 covar 的拟合组和分析组都是 C-alpha，输出特征值和特征向量；anaeig 再将同一选择投影到 PC1、PC2。改变原子组后不能继续复用旧特征向量。
聚类采用 GROMOS 方法、0.12 nm 阈值、每 2 ps 一帧，输出簇大小和代表结构。阈值与原子选择都改变结果，不能脱离这些条件比较“簇的数量”。mdmat 的平均残基距离图用于观察接触，不等同于按阈值统计的接触概率。
这条短轨迹只显示局部涨落。二维 PCA 散点不能直接叫作充分采样的自由能面，未采到的区域也不是无限高势垒。
1 准备与提交#
先用参考轨迹练习；下列命令只复制分析所需的轨迹与输入，不复制已有分析结果。
python new-analysis.py ubiquitin
cd work-analysis/ubiquitin/analysis
sbatch run.sh
若已有自己刚跑完的同名模型，把附件 cases/对应目录/analysis/analyze.sh 和 run.sh 放到新结果目录的 analysis 子目录，并保证上一级包含所需的 TPR、XTC、TRR、EDR 与 GRO。电势分析另需复制 System 索引或自行生成。
| 算例目录 | 本次内容 |
|---|---|
ubiquitin |
完整输入、前置文件与真实计算输出 |
2 实际终端命令#
以下是在对应算例目录中实际执行并生成结果的命令。分析命令在该算例的 analysis 子目录执行；交互回答已用 printf 写入，组名需随自己的体系核对。
ubiquitin/analysis
printf %s 'C-alpha
C-alpha
' | gmx covar -s ../md.tpr -f center.xtc -o eigenvalues.xvg -v eigenvectors.trr -av average.pdb
printf %s 'C-alpha
C-alpha
' | gmx anaeig -s ../md.tpr -f center.xtc -v eigenvectors.trr -eig eigenvalues.xvg -first 1 -last 2 -proj projection.xvg -2d pc12.xvg
printf %s 'C-alpha
Protein
' | gmx cluster -s ../md.tpr -f center.xtc -method gromos -cutoff 0.12 -dt 2 -g cluster.log -o clusters.xpm -sz cluster-size.xvg -cl centers.pdb
printf %s 'Protein
' | gmx mdmat -s ../md.tpr -f center.xtc -mean contacts.xpm -dt 2
3 输入文件与关键设置#
本篇读取计算产生的 TPR/轨迹/能量文件；完整来源随资料包提供。先核对 gmx check -f md.xtc 的时间和帧数，选择表达式中的原子名称可从 md.gro 查看。
4 实际结果与画图#
下表为本次 100 ps 生产轨迹后半段的简单均值，不是收敛物性或误差分析。
| 实例 | 原子数 | 温度 K | 密度 kg/m³ |
|---|---|---|---|
| ubiquitin | 26742 | 299.40 | 994.22 |
用附件中的 Python 脚本重画所有参考图：
# 回到 gromacs-examples 根目录
module load ase/3.29.0
python plot.py --root .
plot.py 使用 NumPy 和 Matplotlib，直接读取 XVG/COLVAR/坐标数值，图片保存在 figures。对单项新结果，可照着脚本对应段替换文件路径；不要把不同模型或单位的数据直接拼在一起。
参数查阅#
以本集群 gmx 命令名 -h 为准；GROMACS 2026.4 命令手册和MDP 参数手册可查完整选项。