计算化学公社

标题: 求助PCA分析后画了自由能景观图，但没观察到低谷 [打印本页]
作者
Author: HanLuyao    时间: 2025-6-23 16:11
标题: 求助PCA分析后画了自由能景观图，但没观察到低谷
1、我对模拟之后的轨迹计算了PCA，然后画了自由能景观图，但是景观图很奇怪，没有观察到能量低谷，这种情况下，通常是什么原因呢？
2、还想请教一下老师，这个自由能景观图计算的是哪些能量呢？
(, 下载次数 Times of downloads: 40) (, 下载次数 Times of downloads: 42) (, 下载次数 Times of downloads: 42) (, 下载次数 Times of downloads: 43)
作者
Author: student0618    时间: 2025-6-23 18:15
http://sobereva.com/73
作者
Author: sobereva    时间: 2025-6-23 19:12
如置顶的新社员必读贴、论坛首页的公告栏、版头的红色大字非常明确所示，求助帖必须在帖子标题明确体现出此帖内容是求助或提问，并清楚、准确反映出帖子具体内容，避免有任何歧义和含糊性，仔细看http://bbs.keinsci.com/thread-9348-1-1.html。我已把你的不恰当标题 “自由能景观图分析” 改了，以后务必注意
作者
Author: sobereva    时间: 2025-6-23 19:15
没有模拟的任何细节、对象的说明，没法回答。只能说可能操作不当，或者取的PC不当，或者当前动力学模拟情况不适合做PCA分析。并且画图取的格点间距偏大

图片体现的是相对自由能
作者
Author: 13277552957    时间: 2025-6-23 20:02
猜测是没有做好周期性矫正，多亚基蛋白跨越盒子，RMSD是否波动跳跃
作者
Author: HanLuyao    时间: 2025-6-23 21:42
sobereva 发表于 2025-6-23 19:15
没有模拟的任何细节、对象的说明，没法回答。只能说可能操作不当，或者取的PC不当，或者当前动力学模拟情况 ...

谢谢老师。我的模拟体系是同源多聚体和小分子，模拟的mdp文件为
  1. integrator          = md
  2. nsteps              = 50000000       ; 100 ns（50000000 步，每步 2 fs）
  3. dt                  = 0.002

  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

  13. tcoupl              = V-rescale
  14. tc-grps             = Protein Non-Protein
  15. tau_t               = 0.1 0.1
  16. ref_t               = 298.15 298.15

  18. pcoupl              = Parrinello-Rahman
  19. pcoupltype          = isotropic
  20. tau_p               = 2.0
  21. ref_p               = 1.0
  22. compressibility     = 4.5e-5

  24. constraints         = h-bonds
  25. pbc                 = xyz

  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

  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

  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

  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

  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

  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

  42. # Step 6: 使用 gmx sham 生成自由能图
  43. echo ">>> Step 6: 生成自由能图..."
  44. gmx sham -f gsham_input.xvg -ls FES.xpm
复制代码
作者
Author: HanLuyao    时间: 2025-6-23 21:43
13277552957 发表于 2025-6-23 20:02
猜测是没有做好周期性矫正，多亚基蛋白跨越盒子，RMSD是否波动跳跃

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



欢迎光临 计算化学公社 (http://ccc.keinsci.com/) Powered by Discuz! X3.3