MMPBSA的MM部分和Gromacs用g_energy算出来的Coul+LJ是一样的吗？

xiaotianzhou

请教一下老师。最近在用MMPBSA计算蛋白和配体之间的结合自由能。但是发现了一些奇怪的现象。
1. MMPBSA的MM部分对有的体系和Gromacs用g_energy算出来的Coul+LJ几乎一样。但对另外一个体系差别又很大，但是能量曲线的形状类似，只是整体上移或者下移了100-200kJ/mol。
2. MMPBSA的MM部分和Gromacs用g_energy的能量，在Coul这一项几乎相等，差别就在LJ上面。我试过相同的体系，分别单独run了2个jobs。第一个job MM和gromcas的能量几乎相等。第二个job MM比gromcas的能量小了将近200kJ/mol。
这到底是什么原因呢？到底哪个能量可以相信?期待老师的解答。万分谢谢！

zheny1379

体系带电荷  用MMPB(GB)SA方法算影响很大，看综述End-Point Binding Free Energy Calculation with MM/PBSA and MM/GBSA: Strategies and Applications in Drug Design

xiaotianzhou

比如图上的能量数据，带-1电荷的小分子的MM+PBSA算出来就是正的（蓝色数据），主要因为带-1电荷的小分子PBlig算的过大，导致总的PB部分过正，最后加上MM的部分后还是个正数。不太清楚这种情况要怎么处理，还是说MM/PBSA就不适合处理带电的体系呢？请假一下大家。谢谢！

xiaotianzhou

lyj714 发表于 2021-2-19 13:48
差别你应该弄错了。主要差异绝对在于静电相互作用，gmx_mmpbsa这个脚本和g_mmpbsa差不多，因为有个个pdie参 ...

目前大致知道差异的原因了。一方面是你说的pdie，这个影响静电的计算。还有Gromacs计算的时候设置了截断值，MMPBSA没有截断值，造成静电和VDW都有差异，但是有的体系，静电和VDW的差异刚好抵消了，所以看总的能量GROMACS和MMPBSA的MM部分一致。
另外还发现了一个问题，对同样的蛋白不同小分子的体系，发现带电荷的小分子往往算出来的结合能不很准，甚至出现正的结合能（MM+PBSA），主要是带电荷的小分子PBlig的部分算的太大了，最终导致总的MMPBSA的能量变成了正的。从MD的估计看带电小分子结合的更好，可是MMPBSA算出来结合能是正的，不知道怎么处理这种问题，小分子带电是否需要修改参数，有没有有经验的朋友赐教一二啊。

lyj714

差别你应该弄错了。主要差异绝对在于静电相互作用，gmx_mmpbsa这个脚本和g_mmpbsa差不多，因为有个个pdie参数的存在，此参数直接影响静电作用的大小，比如你设置为2，那么得到的就和gmx energy差距是两倍的关系。至于vdw，gmx_mmpbsa和g_mmpbsa和gmx energy得到的应该是相同的

xiaotianzhou

脚本是网上下载的。用的这个命令:
./gmx_mmpbsa -f ../md_out.xtc -s ../md.tpr -n ../../index.ndx -com Protein_UNK -pro Protein -lig UNK
gmx_mmpbsa里面的参数我没改，就用的默认值。下面是主要的参数。我也把完整脚本上传到附件了。
apbs=/usr/bin/apbs
export MCSH_HOME=/dev/null                              # APBS io.mc

pid=pid                         # 输出文件($$可免重复) prefix of the output files($$)
scr=_$pid.scr           # 屏幕输出文件 file to save the message from the screen
qrv=_$pid.qrv           # 电荷/半径/VDW参数文件 to save charge/radius/vdw parmeters

radType=1                       # 原子半径类型 radius of atoms (0:ff; 1:mBondi; 2:Bondi)
radLJ0=1.2                      # 用于LJ参数原子的默认半径(A, 主要为H) radius when LJ=0 (H)

meshType=0                      # 网格大小 mesh (0:global  1:local)
gridType=1                      # 格点大小 grid (0:GMXPBSA 1:psize)

cfac=3                          # 分子尺寸到粗略格点的放大系数
fadd=10                         # 分子尺寸到细密格点的增加值(A)
df=.5                           # 细密格点间距(A) The desired fine mesh spacing (A)

# 极性计算设置(Polar)
PBEset='
  temp  300         # 温度
  pdie  2           # 溶质介电常数
  sdie  78.54       # 溶剂介电常数, 真空1, 水78.54

  lpbe              # PB方程求解方法, lpbe(线性), npbe(非线性), smbpe(大小修正)
  bcfl  mdh         # 粗略格点PB方程的边界条件, zero, sdh/mdh(single/multiple Debye-Huckel), focus, map
  srfm  smol        # 构建介质和离子边界的模型, mol(分子表面), smol(平滑分子表面), spl2/4(三次样条/7阶多项式)
  chgm  spl4        # 电荷映射到格点的方法, spl0/2/4, 三线性插值, 立方/四次B样条离散
  swin  0.3         # 立方样条的窗口值, 仅用于 srfm=spl2/4

  srad  1.4         # 溶剂探测半径
  sdens 10          # 表面密度, 每A^2的格点数, (srad=0)或(srfm=spl2/4)时不使用

  ion  1 0.15 0.95  # 阳离子的电荷, 浓度, 半径
  ion -1 0.15 1.81  # 阴离子

  calcforce  no
  calcenergy comps'
# 非极性计算设置(Apolar/Non-polar)
PBAset='
  temp  300  # 温度
  srfm  sacc # 构建溶剂相关表面或体积的模型
  swin  0.3  # 立方样条窗口(A), 用于定义样条表面

  # SASA
  srad  1.4    # 探测半径(A)
  gamma 1      # 表面张力(kJ/mol-A^2)

  #gamma const 0.027     0        # 表面张力, 常数
  #gamma const 0.0226778 3.84928  # 表面张力, 常数

  press  0     # 压力(kJ/mol-A^3)
  bconc  0     # 溶剂本体密度(A^3)
  sdens 10
  dpos  0.2
  grid  0.1 0.1 0.1

  # SAV
  #srad  1.29      # SAV探测半径(A)
  #press 0.234304  # 压力(kJ/mol-A^3)

  # WCA
  #srad   1.25           # 探测半径(A)
  #sdens  200            # 表面的格点密度(1/A)
  #dpos   0.05           # 表面积导数的计算步长
  #bconc  0.033428       # 溶剂本体密度(A^3)
  #grid   0.45 0.45 0.45 # 算体积分时的格点间距(A)

  calcforce no
  calcenergy total'

chemcomp

@xiaotianzhou 社长的意思是你用什么工具、什么参数计算的mmpbsa。知道这些才能进一步分析差异的来源。

sobereva

xiaotianzhou 发表于 2021-2-18 20:47
这里的能量都没有考虑溶剂化。
老师按照您的说法，那MMPBSA里MM的部分应该和Gromacs用g_energy抽出来的 ...

你仍然没说清楚MMPBSA你是具体怎么实现的

xiaotianzhou

sobereva 发表于 2021-2-18 19:38
得说清楚MMPBSA具体怎么实现的
本身MMPBSA的MM部分的静电作用和范德华作用原理上直接就对应于gmx的蛋白和 ...

这里的能量都没有考虑溶剂化。
老师按照您的说法，那MMPBSA里MM的部分应该和Gromacs用g_energy抽出来的一样才对啊。怎么有的体系差的那么多啊？

sobereva

得说清楚MMPBSA具体怎么实现的
本身MMPBSA的MM部分的静电作用和范德华作用原理上直接就对应于gmx的蛋白和配体组之间的静电作用和范德华作用。

nianbin

mmpbsa是考虑了溶剂化自由能的，gmx energy没有考虑这部分