计算化学公社

标题: g_mmpbsa能量分解过程中的异常值 [打印本页]
作者
Author: 510426762    时间: 2018-3-17 04:16
标题: g_mmpbsa能量分解过程中的异常值
本帖最后由 510426762 于 2018-3-17 06:47 编辑

模拟了蛋白质和受体在0.15M NaCl水溶液环境中的结合作用，时间长度100ns。rmsd结果稳定了。
g_mmpbsa分解能量，用的一步法分解，参数如下。
(, 下载次数 Times of downloads: 43)
  1. ;Polar calculation: "yes" or "no"
  2. polar                = yes

  4. ;=============
  5. ;PSIZE options
  6. ;=============
  7. ;Factor by which to expand molecular dimensions to get coarsegrid dimensions.
  8. cfac                 = 1.5

  10. ;The desired fine mesh spacing (in A)
  11. gridspace         = 0.5

  13. :Amount (in A) to add to molecular dimensions to get fine grid dimensions.
  14. fadd                 = 5

  16. ;Maximum memory (in MB) available per-processor for a calculation.
  17. gmemceil         = 4000

  19. ;=============================================
  20. ;APBS kwywords for polar solvation calculation
  21. ;=============================================
  22. ;Charge of positive ions
  23. pcharge         = 1

  25. ;Radius of positive charged ions
  26. prad                = 0.95

  28. ;Concentration of positive charged ions
  29. pconc           = 0.150

  31. ;Charge of negative ions
  32. ncharge         = -1

  34. ;Radius of negative charged ions
  35. nrad                = 1.81

  37. ;Concentration of negative charged ions
  38. nconc                 = 0.150

  40. ;Solute dielectric constant
  41. pdie                 = 2

  43. ;Solvent dielectric constant
  44. sdie                 = 80

  46. ;Reference or vacuum dielectric constant
  47. vdie                 = 1

  49. ;Solvent probe radius
  50. srad                 = 1.4

  52. ;Method used to map biomolecular charges on grid. chgm = spl0 or spl2 or spl4
  53. chgm            = spl4

  55. ;Model used to construct dielectric and ionic boundary. srfm = smol or spl2 or spl4
  56. srfm            = smol

  58. ;Value for cubic spline window. Only used in case of srfm = spl2 or spl4.
  59. swin                 = 0.30

  61. ;Numebr of grid point per A^2. Not used when (srad = 0.0) or (srfm = spl2 or spl4)
  62. sdens                 = 10

  64. ;Temperature in K
  65. temp                 = 300

  67. ;Type of boundary condition to solve PB equation. bcfl = zero or sdh or mdh or focus or map
  68. bcfl                 = mdh

  70. ;Non-linear (npbe) or linear (lpbe) PB equation to solve
  71. PBsolver         = lpbe


  74. ;========================================================
  75. ;APBS kwywords for Apolar/Non-polar solvation calculation
  76. ;========================================================
  77. ;Non-polar solvation calculation: "yes" or "no"
  78. apolar                = yes

  80. ;Repulsive contribution to Non-polar
  81. ;===SASA model ====

  83. ;Gamma (Surface Tension) kJ/(mol A^2)
  84. gamma           = 0.0226778

  86. ;Probe radius for SASA (A)
  87. sasrad          = 1.4

  89. ;Offset (c) kJ/mol
  90. sasaconst       = 3.84928

  92. ;===SAV model===
  93. ;Pressure kJ/(mol A^3)
  94. press           = 0

  96. ;Probe radius for SAV (A)
  97. savrad          = 0

  99. ;Offset (c) kJ/mol
  100. savconst        = 0

  102. ;Attractive contribution to Non-polar
  103. ;===WCA model ====
  104. ;using WCA method: "yes" or "no"
  105. WCA             = no

  107. ;Probe radius for WCA
  108. wcarad          = 1.20

  110. ;bulk solvent density in A^3
  111. bconc                = 0.033428

  113. ;displacment in A for surface area derivative calculation
  114. dpos                = 0.05

  116. ;Quadrature grid points per A for molecular surface or solvent accessible surface
  117. APsdens                = 20

  119. ;Quadrature grid spacing in A for volume integral calculations
  120. grid            = 0.45 0.45 0.45

  122. ;Parameter to construct solvent related surface or volume
  123. APsrfm          = sacc

  125. ;Cubic spline window in A for spline based surface definitions
  126. APswin          = 0.3

  128. ;Temperature in K
  129. APtemp          = 300
复制代码
计算出来的能量，发现energy_MM里有异常值。主要集中在s4 s5 s6,Protein-MOL VdW Energy接近于0，VdW = s4 - s0 - s2
  1. @ s0 legend "Protein VdW Energy"
  2. @ s1 legend "Protein Elec. Energy"
  3. @ s2 legend "MOL VdW Energy"
  4. @ s3 legend "MOL Elec. Energy"
  5. @ s4 legend "Protein-MOL VdW Energy"
  6. @ s5 legend "Protein-MOL Elec. Energy"
  7. @ s6 legend "Protein-MOL Total Energy"
  8.           0.000          0.000          0.000          0.000          0.000         -0.001         -0.342         -0.343
  9.         100.000          0.000          0.000          0.000          0.000         -0.001         -0.437         -0.438
  10.         200.000          0.000          0.000          0.000          0.000       -278.876        -37.203       -316.079
  11.         300.000          0.000          0.000          0.000          0.000       -286.627        -43.203       -329.830
  12.         400.000          0.000          0.000          0.000          0.000         -0.001         -0.025         -0.025
  13.         500.000          0.000          0.000          0.000          0.000       -307.671          3.073       -304.598
  14.         600.000          0.000          0.000          0.000          0.000       -285.399        -26.875       -312.274
  15.         700.000          0.000          0.000          0.000          0.000         -0.001          0.014          0.013
  16.         800.000          0.000          0.000          0.000          0.000         -0.001         -0.132         -0.133
  17.         900.000          0.000          0.000          0.000          0.000         -0.001          0.037          0.036
  18.        1000.000          0.000          0.000          0.000          0.000         -0.001          0.319          0.318
  19.        1100.000          0.000          0.000          0.000          0.000         -0.001          0.156          0.155
  20.        1200.000          0.000          0.000          0.000          0.000         -0.001          0.322          0.321
  21.        1300.000          0.000          0.000          0.000          0.000       -260.126         -5.969       -266.095
  22.        1400.000          0.000          0.000          0.000          0.000         -0.001         -0.495         -0.496
  23.        1500.000          0.000          0.000          0.000          0.000         -0.001         -0.023         -0.023
  24.        1600.000          0.000          0.000          0.000          0.000         -0.001         -0.141         -0.142
  25.        1700.000          0.000          0.000          0.000          0.000         -0.001         -0.035         -0.035
  26.        1800.000          0.000          0.000          0.000          0.000         -0.001         -0.068         -0.069
  27.        1900.000          0.000          0.000          0.000          0.000         -0.001          0.088          0.088
  28.        2000.000          0.000          0.000          0.000          0.000       -264.009         -3.774       -267.782
  29.        2100.000          0.000          0.000          0.000          0.000         -0.001         -0.227         -0.227
  30.        2200.000          0.000          0.000          0.000          0.000         -0.001          0.137          0.136
  31.        2300.000          0.000          0.000          0.000          0.000         -0.001         -0.138         -0.139
  32.        2400.000          0.000          0.000          0.000          0.000         -0.001          0.105          0.105
  33.        2500.000          0.000          0.000          0.000          0.000         -0.001          0.039          0.038
复制代码
(, 下载次数 Times of downloads: 40)

这个算作异常值么？还是由于这个时刻蛋白质和配体的确没有VdW力作用了呢？


大概情况如上
作者
Author: 510426762    时间: 2018-3-17 05:10
好像可能解决了。
看到这篇帖子提示了。http://bbs.keinsci.com/thread-6822-1-1.html
g_mmpbsa里也说明了
WARNING: Trajectory should be PBC corrected and molecule should not be PBC broken. To make molecule whole in trajectory, please follow these links: PBC and trjconv.
作者
Author: HarrisLin    时间: 2018-3-26 14:55
本帖最后由 HarrisLin 于 2018-3-26 14:58 编辑

是PBC的問題沒錯 要用noPBC分析
作者
Author: 510426762    时间: 2018-5-30 11:56
510426762 发表于 2018-3-17 05:10
好像可能解决了。
看到这篇帖子提示了。http://bbs.keinsci.com/thread-6822-1-1.html
g_mmpbsa里也说明 ...

是的，之后仔细阅读，发现漏掉了这句。
用这种方式做了noPBC
echo "0"      | gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_whole.xtc -pbc whole -n index.ndx
echo "1 0"   | gmx trjconv -s md_0_1.tpr -f md_0_1_whole.xtc -o md_0_1_mol_center.xtc -center -pbc mol -n index.ndx
作者
Author: lonemen    时间: 2018-9-14 18:28
学习了，并收藏之。
作者
Author: 芮啦    时间: 2021-6-4 10:27
510426762 发表于 2018-3-17 05:10
**** 作者被禁止或删除 内容自动屏蔽 ****

我现在刚开始学习使用g_mmpbsa来计算结合能，有一点问题希望向您请教一下，请问您的mdp参数文件是采用官网上的参数还是依照自己的体系设定的呢，我的体系中不止有一种阳离子，所以在设定mdp文件中的pcharge、prad等参数时不知道怎么处理，希望向您请教一下，还请您指点一二



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