|
|
马上注册，结交更多好友，下载更多分子模拟资源。
您需要 登录 才可以下载或查看，没有帐号？我想注册
x
本帖最后由 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 /////////////////////////////////// |
评分
-
查看全部评分
|
窥视卡
雷达卡
发表于 2014-8-11 10:56:01
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
千斤顶
显身卡
楼主
。计算方面的话，我也是按照g_mmpbsa教程中计算的。这个软件确实使gromacs使用者受益良多，恩，从相对水平的话，可能靠谱些。总之，谢谢您
