MMPBSA结合能计算的一些经验

casea

之前使用MCPB.py方法构建了金属蛋白体系(教程），并使用Amber进行了模拟。在成品模拟结束之后，像以往的操作步骤一样，计算配体与蛋白的结合能。
Amber和gromacs相比，它自带了计算结合能的mmpbsa.py程序。也有开发者将mmpbsa转移到了gromacs平台(介绍，中文教程）。
如果一切顺利的话，意外就要来了。
使用Amber自带的mmpbsa.py计算gbsa时所得的结合能为正值，且为指数级别。当调用pbsa模块时，可能是因为算法的问题，蛋白内部含有的非标准残基识别不出来(gmx_MMPBSA套用MMPBSA.py的算法，计算MM能量时，会从力场中读取能量)，程序直接报错。
通过网络查询，发现结合能为正的原因大抵是因为一个或者两个组分含有静电荷，气相MM部分的静电相互作用往往会非常大，从而导致总结和能绝对值过于大，这并不是什么计算错误, 而是由MMPBSA这个方法的近似导致的。在mmpbsa.py程序中可以通过引入溶质的介电常数(indi)进行校正，但是没有明确给出数值的适用范围，有文献指出indi介于2-4之间较为合适，但是亲身实验发现，某些体系下indi的改变并不会改善结合能过大的现象。
李老师提供的gmx_mmpbsa.bsh程序提供了一个新的矫正方法进行计算结合能。即使用德拜-休克尔屏蔽方法计算MM: 计算MM时使用考虑离子强度, 使用德拜特征长度对静电作用进行指数衰减. 这样处理后, 所得MM(COU)贡献就变小了, 再加上考虑熵的贡献, 最终所得的结合能总算与实验值比较一致了。同时，脚本计算PB由线性lpbe方法改为非线性npbe方法. 对于净电荷很大的体系, lpbe方法误差过大. 根据维基, 改为npbe后, 所得PB相互作用能会变小.
使用教程见李老师教程，对于gromacs产生的文件可以很容易进行操作，
下面将介绍如何使用gmx_mmpbsa.bsh计算amber产生的文件并计算结合能

• 导出amber轨迹文件为xtc格式
  
  1. cpptraj <<EOF
     
  2. parm xxx.prmtop
     
  3. trajin  xxx.crd
     
  4. trajout gmx.xtc xtc start 9500 stop 9600
     
  5. run
     
  6. quit
     
  7. EOF

  复制代码
  
  
• 使用parmed将parm和crd导出为xtc和top
  
  1. python << EOF
     
  2. import parmed as pmd
     
  3. amber=pmd.load_file("xxx.prmtop","xxx.inpcrd")
     
  4. amber.save("gmx.gro")
     
  5. amber.save("gmx.top")
     
  6. exit()
     
  7. EOF

  复制代码
  
  
• 生成需要的tpr文件
  

  gmx_mmpbsa计算MM和PBSA所需的原子参数(电荷, 半径, LJ参数等)来自GROMACS的.tpr文件, 因此, 需要调用gmx dump程序处理.tpr文件, 以便抽取所需的参数. 采用这种方法, 可以避免版本不兼容的问题, 因此也就可以支持任意版本的GROMACS. 这里只要生成tpr文件即可，也即mdp文件任意

  1. gmx grompp -f em.mdp -c gmx.gro -p gmx.top -maxwarn 5 -o gmx.tpr

  复制代码
  
  
• 制作需要的index.ndx文件
  
  1. gmx make_ndx -f gmx.gro  <<EOF
     
  2. 1|13|14|15|16|17|18|19|20
     
  3. name 26 protein
     
  4. 26|21
     
  5. name 27 complex
     
  6. name 21 lig
     
  7. q
     
  8. EOF

  复制代码
  
  
• 运行gmx_mmpbsa.bsh文件
  
  1. bash gmx_mmpbsa.bsh

  复制代码
  
  
  

tongyu

你好，我想请教下使用该方法写文章如何引用？相关文章没有找到

大村驴

“2. 使用parmed将parm和crd导出为xtc和top”应该是gro和top吧

mwx

您好，我在计算时一直提示awk: 命令行:207: (FILENAME=_pid.pdb FNR=3686023) 警告： exp：参数 -1246.16 超出范围。这个问题怎么解决呢？