求助PCA分析后画了自由能景观图，但没观察到低谷

HanLuyao

1、我对模拟之后的轨迹计算了PCA，然后画了自由能景观图，但是景观图很奇怪，没有观察到能量低谷，这种情况下，通常是什么原因呢？
2、还想请教一下老师，这个自由能景观图计算的是哪些能量呢？

HanLuyao

13277552957 发表于 2025-6-23 20:02
猜测是没有做好周期性矫正，多亚基蛋白跨越盒子，RMSD是否波动跳跃

我做了周期性矫正，之前计算的RMSD也很正常

HanLuyao

sobereva 发表于 2025-6-23 19:15
没有模拟的任何细节、对象的说明，没法回答。只能说可能操作不当，或者取的PC不当，或者当前动力学模拟情况 ...

谢谢老师。我的模拟体系是同源多聚体和小分子，模拟的mdp文件为

1. integrator          = md
   
2. nsteps              = 50000000       ; 100 ns（50000000 步，每步 2 fs）
   
3. dt                  = 0.002
   
4. 
   
5. cutoff-scheme       = Verlet
   
6. nstlist             = 10
   
7. rlist               = 1.0
   
8. vdwtype             = cutoff
   
9. rvdw                = 1.0
   
10. coulombtype         = PME
    
11. rcoulomb            = 1.0
    
12. 
    
13. tcoupl              = V-rescale
    
14. tc-grps             = Protein Non-Protein
    
15. tau_t               = 0.1 0.1
    
16. ref_t               = 298.15 298.15
    
17. 
    
18. pcoupl              = Parrinello-Rahman
    
19. pcoupltype          = isotropic
    
20. tau_p               = 2.0
    
21. ref_p               = 1.0
    
22. compressibility     = 4.5e-5
    
23. 
    
24. constraints         = h-bonds
    
25. pbc                 = xyz
    
26. 
    
27. nstxout             = 5000           ; 每 10 ps 保存一次坐标
    
28. nstvout             = 5000           ; 每 10 ps 保存一次速度
    
29. nstenergy           = 5000           ; 每 10 ps 记录能量
    
30. nstlog              = 5000           ; 每 10 ps 记录日志

复制代码

下面是我计算PC的脚本

1. # Step 1: 去除PBC并居中轨迹
   
2. echo ">>> Step 1: 去除PBC并居中..."
   
3. cd "$INPUT_DIR"
   
4. gmx trjconv -s md.tpr -f md.trr -o md_noPBC.trr -center -pbc mol -ur compact << EOF
   
5. 1
   
6. 1
   
7. EOF
   
8. 
   
9. # Step 2: 对齐轨迹
   
10. echo ">>> Step 2: 对齐轨迹..."
    
11. gmx trjconv -s md.tpr -f md_noPBC.trr -o md_fit.trr -fit rot+trans << EOF
    
12. 1
    
13. 1
    
14. EOF
    
15. 
    
16. # Step 3: 计算协方差矩阵
    
17. echo ">>> Step 3: PCA - 计算协方差矩阵..."
    
18. cd "$PCA_DIR"
    
19. gmx covar -s "$INPUT_DIR/md.tpr" -f "$INPUT_DIR/md_fit.trr" -o eigenvalues.xvg -v eigenvectors.trr -xpma covar_matrix.xpm << EOF
    
20. 3
    
21. 3
    
22. EOF
    
23. 
    
24. # Step 4: 投影 PC1 和 PC2
    
25. echo ">>> Step 4: 计算主成分投影..."
    
26. gmx anaeig -s "$INPUT_DIR/md.tpr" -f "$INPUT_DIR/md_fit.trr" -v eigenvectors.trr -first 1 -last 1 -proj pc1.xvg << EOF
    
27. 3
    
28. 3
    
29. EOF
    
30. 
    
31. gmx anaeig -s "$INPUT_DIR/md.tpr" -f "$INPUT_DIR/md_fit.trr" -v eigenvectors.trr -first 2 -last 2 -proj pc2.xvg << EOF
    
32. 3
    
33. 3
    
34. EOF
    
35. 
    
36. # Step 5: 合并 PC1 与 PC2
    
37. echo ">>> Step 5: 合并PC1和PC2数据..."
    
38. awk '/^[^@#]/ {print $2}' pc2.xvg > pc2_values.txt
    
39. awk '/^[^@#]/ {print $1, $2}' pc1.xvg > pc1_clean.xvg
    
40. paste pc1_clean.xvg pc2_values.txt > gsham_input.xvg
    
41. 
    
42. # Step 6: 使用 gmx sham 生成自由能图
    
43. echo ">>> Step 6: 生成自由能图..."
    
44. gmx sham -f gsham_input.xvg -ls FES.xpm

复制代码

13277552957

猜测是没有做好周期性矫正，多亚基蛋白跨越盒子，RMSD是否波动跳跃

sobereva

没有模拟的任何细节、对象的说明，没法回答。只能说可能操作不当，或者取的PC不当，或者当前动力学模拟情况不适合做PCA分析。并且画图取的格点间距偏大

图片体现的是相对自由能

sobereva

如置顶的新社员必读贴、论坛首页的公告栏、版头的红色大字非常明确所示，求助帖必须在帖子标题明确体现出此帖内容是求助或提问，并清楚、准确反映出帖子具体内容，避免有任何歧义和含糊性，仔细看http://bbs.keinsci.com/thread-9348-1-1.html。我已把你的不恰当标题 “自由能景观图分析” 改了，以后务必注意

student0618

http://sobereva.com/73