GROMACS 中 g_mmpbsa 工具使用详解
本帖最后由 solo7773 于 2014-8-14 11:20 编辑注：坛中已有前辈贡献类似帖子，介绍的比较详细专业。在此献上入门版，疏漏错讹之处敬请指正^_^
使用 g_mmpbsa 可以像在amber中一样，使用pbsa方法计算能量（自由能，结合能，结合自由能？）。
本教程实用平台，亲测 Ubuntu 12.04 LTS Server，Ubuntu Kylin，Deepin linux (Ubuntu 12.04 based)
|----------------------------------------------------------------------|
||必须环境之 GROMACS和 APBS 安装，Ubuntu系统，方便省事||
||sudo apt-get install gromacs ||
||sudo apt-get install apbs ||
|----------------------------------------------------------------------|
g_mmpbsa官方来源
http://rashmikumari.github.io/g_mmpbsa/
一 安装
参考网页
http://rashmikumari.github.io/g_mmpbsa/Download-and-Installation.html
1. 推荐使用"Pre-compiled executable program"，程序已经是编译好了的，即可运行
1.1 复制程序到 /usr/local/bin 文件夹
sudo cp g_mmpbsa /usr/local/bin
1.2 如果不能运行，添加可执行属性
sudo chmod +x /usr/local/bin/g_mmpbsa
1.3 用同样的方法按需添加其他文件（or工具，or预编译程序）
2. Installation from Source-Code，源代码安装，略麻烦，想折腾可以试试
二 使用
参考网页
http://rashmikumari.github.io/g_mmpbsa/Usage.html#top
1. 明白binding free energy 的组成，对于 g_mmpbsa中的一些表示，结合着amber中的一些解释来理解
1) "MM" terms
molecular mechanics (vdw and electrostatic) vacuum energy
2) "GB" term
polar solvation energy "GB"
3) "SA" term [可以三中方法选一：SASA, SAV and Weeks–Chandler–Andersen (WCA)]non-polar solvation energy
2. 对于1中所说的3个能量组成，g_mmpbsa是先分别算好，然后用工具 MmPbSaStat.py, calculating the average binding energy and its standard deviation/error from the output files obtained from g_mmpbsa
2.0) 轨迹处理
由于在计算能量过程中，是安每ps来计算，特别费时，所以需要用trjconv把轨迹文件，通过-skip参数设置每隔多长时间读取一次轨迹，如
my command: trjconv -f md.xtc -o trj_8-10ns -b 8000 -e 10000 -skip 100
最后得到的文件，似乎时间序列应是8100, 8200, 8300, ..., 10000，共计20个frames。但是，请注意
在我的动力学过程中是每2ps记录一次轨迹，所以实际上是，这里得到的8-10ns的xtc文件中时间是
8000, 8200, 8400, ..., 10000，共计10个frames
如果是每1ps，2ps，或4ps记录一次轨迹，实际情况发生相应变化
-------------------- Running log ---------------------
Reading frame 0 time 8000.000
Precision of ../md_noPBC.xtc is 0.001 (nm)
Using output precision of 0.001 (nm)
Last frame 1000 time 10000.000 ->frame 10 time 10000.000
---------------------- End Log ----------------------------
对于以上log的解释：(10000-8000)=2000 ps，步长2 ps，2000/2=1000 frames，由于skip为100，所以1000/100=10 frames
2.1) calculate "MM" terms
sample command:g_mmpbsa -f traj.xtc -s topol.tpr -n index.ndx -mme -mm energy_MM.xvg -decomp -mmcon contrib_MM.dat
my command: g_mmpbsa -f ../md_noPBC.xtc -s ../../md.tpr -mme -mm energy_MM.xvg
注：-mme表示计算 van der Waals and electrostatic energy
选择group，分别为 1 protein；13 LIG
2.2) calculate "PB" term, i.e. polar solvation
因为 -pbsa 打开，solvation energy will be calculated and solvation energy parameter input file -i is necessary for this calculation
sample command:g_mmpbsa -f traj.xtc -s topol.tpr -i mmpbsa.mdp -n index.ndx -nomme -pbsa -decomp -pol polar.xvg -pcon contrib_pol.dat
my command: g_mmpbsa -f ../md_noPBC.xtc -s ../../md.tpr -i mmpbsa_polar.mdp -nomme -pbsa -pol polar.xvg
注: still don't know how to set the input .mdp file
-nomme 因为2.1)已经算了gas phase的能量了，也就是"MM" terms
增加了 -pbsa，默认是 -nopbsa，所以在2.1)中没给出类似的参数
必要输入.mdp文件见附件一 选择gruop，分别为 1 protein；13 LIG
2.3) calculate "SA" term, i.e. nonpolar solvation
sample command: g_mmpbsa -f traj.xtc -s topol.tpr -i mmpbsa.mdp -n index.ndx -nomme -pbsa -decomp -apol apolar.xvg -apcon contrib_apol.dat
my command: g_mmpbsa -f ../md_noPBC.xtc -s ../../md.tpr -i mmpbsa_nonpolar_sasa.mdp -nomme -pbsa -apol apolar.xvg
注：与2.2)的细微差别 选择gruop，分别为 1 protein；13 LIG
必要输入.mdp文件见附件二
注：对于2.1), 2.2)和2.3)如果出现运行错误，提示group之类的问题，请加上 -n 参数，提供index文件，使用如下类似命令
g_mmpbsa -f ../md_noPBC.xtc -s ../../md.tpr -n ../../index.ndx -mme -mm energy_MM.xvg
2.4) 运行统计工具 MmPbSaStat.py
由于缺少相应的model，需要先安装 numpy, scipy
2.4.1) 比如，安装numpy
方法一，推荐，参考http://kelly-1122.blog.163.com/b ... 072010711114123703/
Python Numpy Scipy搭建过程详解
在ubuntu里可以是用apt-get命令非常简单的就完成安装
sudo apt-get install python-numpy
sudo apt-get install python-scipy
方法二
go to this website to download correct package:
NumPy Home Page
http://numpy.scipy.org/
http://sourceforge.net/projects/numpy/files/
unzip the package
go to the document
use this command to install numpy: python setup.py install
方法三，参考http://blog.163.com/yang_jianli/ ... 006201162152724339/
setuptools工具安装
a. 使用ez_setup.py安装setuptools
这是 setuptools 自豪的一种安装方式，只需要一个大约 8K 作为的脚本ez_setup.py，就能自动为用户安装包括 setuptools 自身在内的许多 Python 包。使用这种方式，用户只需要下载 ez_setup.py 并运行，就可以自动下载和安装适合用户当前 Python 版本的适当的 setuptools egg 文件(当然，用户需要 Python 2.3.5 以上的版本，64 位操作系统的用户则需要 Python 2.4 以上的版本)。此外，这段脚本还会将可执行的 easy_install 脚本安装到用户所有的操作系统 Python 可执行脚本正常应该安装的位置(例如，Windows 用户会安装到 Python 安装目录下的 Scripts 目录中)。关于这种安装方法的更详细说明和注意事项，请参考其官方说明（见扩展阅读）。简单的安装命令如下：wget -q ez_setup。py下载地址（见扩展阅读） 安装完后，最好确保
b. 使用完整的安装包安装setuptools
当然，用户也可以直接使用 setuptools发布版本来安装。对于使用 Windows 的用户，这也是挺方便的方法，许多 Linux 发行版的官方包管理仓库都包含 setuptools 的某个版本。例如，如果你跟我一样使用 Ubuntu ，那安装 setuptools 只是简单的进行如下操作：
# apt-get install python-setuptools
安装模块
easy_install package-name，比如 easy_install pylab
模块卸载
easy_install -m package-name， 比如easy_install -m pylab
easy_install -m 包名，可以卸载软件包，但是卸载后还要手动删除遗留文件。
setuptools它可以自动的安装模块，只需要你提供给它一个模块名字就可以，并且自动帮你解决模块的依赖问题。一般情况下用setuptools给安装的模块会自动放到一个后缀是.egg的目录里。
在Windows里，easy_install这个命令在python安装目录下的scripts里面，所以需要把scripts加到环境变量的PATH里，这样用起来就更方便，linux下不需要注意这个问题。
2.4.2) 运行命令
my command: MmPbSaStat.py -m energy_MM.xvg -p polar.xvg -a apolar.xvg -of full_energy.dat -os summary_energy.dat
注：此处我把这个MmPbSaStat.py已经放到了/usr/local/bin目录中，所有可以直接运行这个python脚本，否则需要执行类似于如下命令
example command: python /此脚本的完整路径/MmPbSaStat.py -m energy_MM.xvg -p polar.xvg -a apolar.xvg -of full_energy.dat -os summary_energy.dat
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ 附件一 mmpbsa_polar.mdp///////////////////////////////////
;Polar calculation: "yes" or "no"
polar = yes
;=============
;PSIZE options
;=============
;Factor by which to expand molecular dimensions to get coarsegrid dimensions.
cfac = 1.5
;The desired fine mesh spacing (in A)
gridspace = 0.5
:Amount (in A) to add to molecular dimensions to get fine grid dimensions.
fadd = 5
;Maximum memory (in MB) available per-processor for a calculation.
gmemceil = 4000
;=============================================
;APBS kwywords for polar solvation calculation
;=============================================
;Charge of positive ions
pcharge = 1
;Radius of positive charged ions
prad = 0.95
;Concentration of positive charged ions
pconc = 0.150
;Charge of negative ions
ncharge = -1
;Radius of negative charged ions
nrad = 1.81
;Concentration of negative charged ions
nconc = 0.150
;Solute dielectric constant
pdie = 2
;Solvent dielectric constant
sdie = 80
;Reference or vacuum dielectric constant
vdie = 1
;Solvent probe radius
srad = 1.4
;Method used to map biomolecular charges on grid. chgm = spl0 or spl2 or spl4
chgm = spl4
;Model used to construct dielectric and ionic boundary. srfm = smol or spl2 or spl4
srfm = smol
;Value for cubic spline window. Only used in case of srfm = spl2 or spl4.
swin = 0.30
;Numebr of grid point per A^2. Not used when (srad = 0.0) or (srfm = spl2 or spl4)
sdens = 10
;Temperature in K
temp = 300
;Type of boundary condition to solve PB equation. bcfl = zero or sdh or mdh or focus or map
bcfl = mdh
;Non-linear (npbe) or linear (lpbe) PB equation to solve
PBsolver = lpbe
;========================================================
;APBS kwywords for Apolar/Non-polar solvation calculation
;========================================================
;Non-polar solvation calculation: "yes" or "no"
apolar = no
;Repulsive contribution to Non-polar
;===SASA model ====
;Gamma (Surface Tension) kJ/(mol A^2)
gamma = 0.0226778
;Probe radius for SASA (A)
sasrad = 1.4
;Offset (c) kJ/mol
sasaconst = 3.84928
;===SAV model===
;Pressure kJ/(mol A^3)
press = 0
;Probe radius for SAV (A)
savrad = 0
;Offset (c) kJ/mol
savconst = 0
;Attractive contribution to Non-polar
;===WCA model ====
;using WCA method: "yes" or "no"
WCA = no
;Probe radius for WCA
wcarad = 1.20
;bulk solvent density in A^3
bconc = 0.033428
;displacment in A for surface area derivative calculation
dpos = 0.05
;Quadrature grid points per A for molecular surface or solvent accessible surface
APsdens = 20
;Quadrature grid spacing in A for volume integral calculations
grid = 0.45 0.45 0.45
;Parameter to construct solvent related surface or volume
APsrfm = sacc
;Cubic spline window in A for spline based surface definitions
APswin = 0.3
;Temperature in K
APtemp = 300
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ 附件一 END, version 20140807, N. Zhou, SCU ///////////////////////////////////
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ 附件二 mmpbsa_nonpolar_sasa.mdp ///////////////////////////////////
;Polar calculation: "yes" or "no"
polar = no
;=============
;PSIZE options
;=============
;Factor by which to expand molecular dimensions to get coarsegrid dimensions.
cfac = 1.5
;The desired fine mesh spacing (in A)
gridspace = 0.5
:Amount (in A) to add to molecular dimensions to get fine grid dimensions.
fadd = 5
;Maximum memory (in MB) available per-processor for a calculation.
gmemceil = 4000
;=============================================
;APBS kwywords for polar solvation calculation
;=============================================
;Charge of positive ions
pcharge = 1
;Radius of positive charged ions
prad = 0.95
;Concentration of positive charged ions
pconc = 0.150
;Charge of negative ions
ncharge = -1
;Radius of negative charged ions
nrad = 1.81
;Concentration of negative charged ions
nconc = 0.150
;Solute dielectric constant
pdie = 2
;Solvent dielectric constant
sdie = 80
;Reference or vacuum dielectric constant
vdie = 1
;Solvent probe radius
srad = 1.4
;Method used to map biomolecular charges on grid. chgm = spl0 or spl2 or spl4
chgm = spl4
;Model used to construct dielectric and ionic boundary. srfm = smol or spl2 or spl4
srfm = smol
;Value for cubic spline window. Only used in case of srfm = spl2 or spl4.
swin = 0.30
;Numebr of grid point per A^2. Not used when (srad = 0.0) or (srfm = spl2 or spl4)
sdens = 10
;Temperature in K
temp = 300
;Type of boundary condition to solve PB equation. bcfl = zero or sdh or mdh or focus or map
bcfl = mdh
;Non-linear (npbe) or linear (lpbe) PB equation to solve
PBsolver = lpbe
;========================================================
;APBS kwywords for Apolar/Non-polar solvation calculation
;========================================================
;Non-polar solvation calculation: "yes" or "no"
apolar = yes
;Repulsive contribution to Non-polar
;===SASA model ====
;Gamma (Surface Tension) kJ/(mol A^2)
gamma = 0.0226778
;Probe radius for SASA (A)
sasrad = 1.4
;Offset (c) kJ/mol
sasaconst = 3.84928
;===SAV model===
;Pressure kJ/(mol A^3)
press = 0
;Probe radius for SAV (A)
savrad = 0
;Offset (c) kJ/mol
savconst = 0
;Attractive contribution to Non-polar
;===WCA model ====
;using WCA method: "yes" or "no"
WCA = no
;Probe radius for WCA
wcarad = 1.20
;bulk solvent density in A^3
bconc = 0.033428
;displacment in A for surface area derivative calculation
dpos = 0.05
;Quadrature grid points per A for molecular surface or solvent accessible surface
APsdens = 20
;Quadrature grid spacing in A for volume integral calculations
grid = 0.45 0.45 0.45
;Parameter to construct solvent related surface or volume
APsrfm = sacc
;Cubic spline window in A for spline based surface definitions
APswin = 0.3
;Temperature in K
APtemp = 300
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ 附件二 END, version 20140807, N. Zhou, SCU ///////////////////////////////////
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ 附件三 my example result ///////////////////////////////////
#Complex Number: 1
===============
SUMMARY
===============
van der Waal energy = -101.565 +/- 11.709 kJ/mol
Electrostattic energy = -329.161 +/- 36.608 kJ/mol
Polar solvation energy = 322.298 +/- 45.109 kJ/mol
SASA energy = -13.990 +/- 1.021 kJ/mol
SAV energy = 0.000 +/- 0.000 kJ/mol
WCA energy = 0.000 +/- 0.000 kJ/mol
Binding energy = -122.418 +/- 24.559 kJ/mol
===============
END
===============
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ 附件三 END, version 20140814, N. Zhou, SCU /////////////////////////////////// 本帖最后由 zhuimeng0812 于 2014-10-17 20:48 编辑
我有问题请教楼主：
我用g_mmpbsa计算蛋白和小分子的结合能量，
算non-polar用的是SASA模型，但是算出来数值偏大好多，是不是哪里的细节问题没有注意到，
求指教
用到的命令：g_mmpbsa -f 7-8ns.xtc -s md.tpr -i nonpolar-wca.mdp -n mm.ndx -nomme -pbsa -decomp -apolapolar.xvg -apcon apolar.dat
计算结果：（一组）van der Waal energy = -193.035 +/- 9.768 kJ/mol
Electrostattic energy = -117.392 +/- 14.038 kJ/ mol
Polar solvation energy = 169.845 +/- 8.646 kJ/mol
SASA energy = 19.360 +/- 12.209 kJ/mol
SAV energy = 0.000 +/- 0.000 kJ/mol
WCA energy = 0.000 +/- 0.000 kJ/mol
Binding energy = -121.223 +/- 17.473 kJ/mol
（二组）van der Waal energy = -280.684 +/- 9.631 kJ/mol
Electrostattic energy = -136.934 +/- 15.977 kJ/mol
Polar solvation energy = 203.455 +/- 8.377 kJ/mol
SASA energy = -245.935 +/- 18.634 kJ/mol
SAV energy = 0.000 +/- 0.000 kJ/mol
WCA energy = 0.000 +/- 0.000 kJ/mol
Binding energy = -460.098 +/- 20.665 kJ/mol
我测试了两个配体。结果都是SASA energy 这一项，出了问题。
我是不是哪里的小细节，没注意到？
请赐教！
楼主你好，小弟是gromacs初学者，在使用g_mmpbsa计算结合自由能时出现了下面的问题，想向你请教：1、我的蛋白是个复合物（由3条链组成，A链是抗原，H链与L链是一个2聚体），现在我想看看A链与H链之间的结合能，行吗？
2、我在运行在使用mmpbsa分析的时候（附件1）出现这个错误（附件2），该如何解决？3、之后发现ndex.ndx 是通过（附件3）得到，运行这里（附件4）改输入什么或者是选择什 么？
4、假若我前面2、3个问题的做法是错误的，那么第1个问题，该如何解决，谢谢。
附件：
solo7773 发表于 2014-8-14 11:19
你好，迟复为歉。你的这个问题，确实也困扰着我。我测试了用amber算和g_mmpbsa来算的差别，比如amber算出 ...
楼上正解！！amber中的能量单位是kcal/mol；gromcas的能量单位是kj/mol；1kcal/mol=4.18kj/mol。另外这两个软件的长度单位，有的地方是Å，有的地方是nm，也要注意；1Å=0.1nm。 我靠，这个贴子必须顶啊！！！！！！！！！！！！！！！！ nice。。 您好，看到您发的g_mmpbsa的帖子。我有一个小问题想请教您：tutorial中1EBZ的ΔGbind=-210KJ/mol。而实验中1EBZ的Ki=0.4nM,根据ΔGbind = RT ln IC50 计算得到的ΔGbind=-53.973KJ/mol。这两个值也差的太多了。本人刚学此软件，还望大神指导:loveliness: 哆啦A梦 发表于 2014-8-13 08:15 static/image/common/back.gif
您好，看到您发的g_mmpbsa的帖子。我有一个小问题想请教您：tutorial中1EBZ的ΔGbind=-210KJ/mol。而实验中 ...
你好，迟复为歉。你的这个问题，确实也困扰着我。我测试了用amber算和g_mmpbsa来算的差别，比如amber算出来只有-54，而g_mmpbsa算出来却有-200左右。我也不是很懂，我个人认为，可能依据你的计算方式，amber里面的pbsa会更准确一些。而对于此g_mmpbsa而言，作者提供了一个在gromacs里面计算pbsa的工具，在一个相对的水平上是可以认为是正确的，但是不一定像amber一样优秀。 谢谢解答，:)。计算方面的话，我也是按照g_mmpbsa教程中计算的。这个软件确实使gromacs使用者受益良多，恩，从相对水平的话，可能靠谱些。总之，谢谢您 zhuimeng0812 发表于 2014-10-16 16:48
我有问题请教楼主：
我用g_mmpbsa计算蛋白和小分子的结合能量，
算non-polar用的是SASA模型，但是算出来数 ...
你好，请问数值偏大很多是什么意思？和谁比偏大？ 和amber的mmpbsa计算值比较 zhuimeng0812 发表于 2014-10-30 08:35
和amber的mmpbsa计算值比较
单位不一样哈。一个是Kcal，一个是Kj/mol
页:
[1]
2