gmx MMPBSA

前言

Jerkwin老师讲的很详细，看后面之前先读这个：

GMXPBSA和g_mmPBSA的比较

GMXPBSA只能给出总体的值，g_mmPBSA可给出残基分解值；
据说使用默认参数时, GMXPBSA在精度上优于g_mmpbsa而且使用也更方便；
GMXPBSA不需要安装，解压即可使用，g_mmpbsa需编译故稍微复杂，但新版本也有编译好可直接使用的；
GMXPBSA有一些问题，比如gmxpbsa给的mdp太老会报错，需要使用时修改代码，如参照em_vac.mdp修改print_files.dat文件的Energy in gas的Mm.mdp部分；
g_mmpbsa只需要xtc、tpr和ndx文件，GMXPBSA若处理有自定义拓扑或残基的体系时需要自带top、itp和/或力场文件夹ff;
GMXPBSA的debug过程很恶心，记录错误的文件在子文件夹里，而且重新运行时需要删除原来生成文件夹。

虽然GMXPBSA的pdie默认为2，但只在计算PBSA时调用，而在计算COU时并未调用pdie，所以若GMXPBSA算出的COU为g_mmPBSA的两倍，则可能是因为计算COU时GMXPBSA使用的pdie为1而g_mmPBSA使用的pdie为2所导致。

玺洛克用GMXMMPBSA（Jerkwin老师开发的脚本）、GMXPBSA以及g_mmPBSA对g_mmPBSA中的算例1EBZ进行了测试和对比：链接

GMXPBSA

参考：GMXPBSAtool安装与使用经验教程

鉴于之前遇到的问题，今更新一下GMXPBSA的安装和使用步骤：

安装GMX和APBS

GMX安装不赘述；
APBS的安装依据参考资料，完成后在 .bashrc 里添加APBS的环境变量：
e.g.

#Path of APBS
export PATH=$PATH:/home/users/ntu/n1805727/APBS/bin
#add the APBS lib directory to my LD_LIBRARY_PATH variable
export LD_LIBRARY_PATH=/home/users/ntu/n1805727/APBS/lib:$LD_LIBRARY_PATH

解压缩GMXPBSA并设置环境变量

将GMXPBSA解压到目标文件夹；
在 .bashrc 里添加GMXPBAS的环境变量：
e.g

#Home of GMXPBSA, Position of gmxpbsa0.sh etc.
export GMXPBSAHOME=/home/users/ntu/n1805727/GMXPBSAtool

修订gmxpbsa0.sh文件

原文件调用gmx指令时缺少前缀gmx，现在添加上：
Line122(修改前):

	 pdb2gmx=pdb2gmx; trjconv=trjconv; mdrun=mdrun; grompp=grompp; editconf=editconf; tpbconv=tpbconv;

Line122(修改后):

	 pdb2gmx='gmx pdb2gmx'; trjconv='gmx trjconv'; mdrun='gmx mdrun'; grompp='gmx grompp'; editconf='gmx editconf'; tpbconv='gmx tpbconv';

修改print_files.dat文件

原文件写出的Mm.mdp缺少参数，运行时会出错，故在L41-L44添加上相应参数: e.g.

function PRINT_MINfile_n {
# echo -e "title = Energy in gas \ncpp   = /usr/bin/cpp \nconstraints = none \nnsteps = 1000 \nnstlist = 1 \nns_type = grid
#rlist = 1 \nrcoulomb = 1 \nrvdw = 1 \npbc = xyz \n;" >> Mm.mdp
echo -e "title   =     EM in Vacuum        \ndefine  =     -DFLEXIBLE          \nintegrator       = steep          \ndt               = 1E-3           \nnsteps           = 15000          \nemtol                   = 250     \nnstlist                  = 1      \nrlist                    = 1      \ncutoff-scheme            = Verlet \nns_type                  = grid   \npbc                      = xyz    \nrvdw                = 1.2         \nrcoulomb            = 1.2         \nvdw-type            = Cut-off     \ncoulombtype         = PME         \nconstraints           = none      \n" >> Mm.mdp
}

注意，此处写出的Mm.mdp文件不要有能量组项（energygrps）（如：energygrps=Protein）。因为脚本会将蛋白和配体分别grompp做计算，计算蛋白时或许可以从蛋白中找到index中定义为Protein的索引所指向的原子而不会报错，但计算配体时必然无法找到该索引所指向的原子，故而报错。

准备模拟文件

将轨迹文件处理后命名为npt.xtc
将tpr文件命名为npt.tpr
将ndx文件命名为index.ndx，内容后面会说。
将上述三个文件放入文件夹1EBZ2中。

修改Input.dat文件

在1EZB2文件夹外面创建（或复制）Inpur.dat文件并修改:

基本设置：

root      1EBZ2
multitrj      n
root_multitrj
run                     1                               #options: integer
RecoverJobs            y                                #options: y,n

修改部分：

将 root 内容修改为包括tpr,xtc,ndx的指定文件夹

路径部分

The INPUT.dat file in each of the EXAMPLE directories should be edited to reflect the location of your Coulomb, APBS and Gromacs binaries.
Cpath path to Coulomb binary ;apbs的coulumb工具
Apath path to APBS binary ;apbs目录
Gpath path to Gromacs binary ;gromacs目录
If coulomb, in the INPUT.dat file, is set to gmx, the coulomb library is not required.

e.g.

Cpath                  /home/users/ntu/n1805727/APBS/share/tools/manip
Apath                  /home/users/ntu/n1805727/APBS/bin
Gpath		    /home/users/ntu/n1805727/GMX5.1.4_Mine/gromacs-5.1.4-cpu/bin

力场部分

use_topology            y                               #options: y,n
	itp_receptor  recept.itp
	itp_receptor  recept2.itp
	itp_ligand    ligand.itp
use_nonstd_ff           n                               #options: y,n
ffield                 6

如果包含有pdb2gmx不能识别的结构，则需要通过 use_topology 来使用自定义拓扑，注意要把top文件和其内容中包括的itp文件都一起放在根目录里；
把要分析的两个片段的拓扑文件（top或itp）提取成itp后分别命名为itp_receptor和itp_ligand，脚本会将这几个itp文件组合，若缺少某个itp文件则会提示top文件中原子个数不够；
itp_receptor可以定义多次，应用于蛋白质多条链的情况；
应该也可以通过 use_nonstd_ff 修改力场文件实现上一步，加原子类型与电荷,但是自己修改不成功；
ffield我选的amber6力场，即选择力场时的第6个选项，默认为“amber99sb-ildn”，若添加过自定义力场则序号数字可能会变动。

Gromacs变量

complex	                Complex
receptor                C1
ligand                  C2
skip                    1                               #options: integer
min                     n                               #options: y,n
double_p                n                               #options: y,n

在index.ndx文件里将两个要分析的片段分别命名为C1和C2，并将二者组成的复合物命名为Complex。

APBS变量

coulomb                  coul                     #options: coul,gmx
linearized              n                               #options: y,n
precF                   1                               #options: integer 1,2,3(..)
temp                    300
bcfl                    sdh                             #options: sdh,mdh,focus
pdie                    2                               #options: integer, usually between 1 and 20

未作修改

其他设置

      #QUEQUE VARIABLE
	
cluster                 n                             #options: y,n
Q                       ...                                     #necessary only if cluster="y"!!
mnp                     1                             #options: integer 
       #OUTPUT VARIABLE
pdf                     n                               #options: y,n
       # ALANINE SCANNING
cas                     n                               #options: y,n
	
#MUTATION        P53          112             PHE lig PHE112ALA
#MUTATION        P53          115             LEU lig LEU115ALA
#MUTATION        P53          116             TRP lig TRP116ALA
#MUTATION        P53          119             LEU lig LEU119ALA

未作修改。

完整input.dat文件

完整input.dat如下：

        #GENERAL

root      1EBZ2

multitrj      n
root_multitrj

run                     1                               #options: integer

RecoverJobs            y                                #options: y,n

Cpath                  /home/users/ntu/n1805727/APBS/share/tools/manip	
Apath                  /home/users/ntu/n1805727/APBS/bin
Gpath		    /home/users/ntu/n1805727/GMX5.1.4_Mine/gromacs-5.1.4-cpu/bin


	#FFIELD

use_topology            y                               #options: y,n
	itp_receptor  recept.itp
	itp_ligand    ligand.itp
use_nonstd_ff           n                               #options: y,n
ffield                 6


        #GROMACS VARIABLE

complex	                System
receptor                C1
ligand                  C2

skip                    1                               #options: integer
min                     n                               #options: y,n

double_p                n                               #options: y,n


        #APBS VARIABLE

coulomb                  coul                     #options: coul,gmx

linearized              n                               #options: y,n
precF                   1                               #options: integer 1,2,3(..)
temp                    300
bcfl                    sdh                             #options: sdh,mdh,focus
pdie                    2                               #options: integer, usually between 1 and 20


       #QUEQUE VARIABLE
	
cluster                 n                             #options: y,n
Q                       ...                                     #necessary only if cluster="y"!!

mnp                     1                             #options: integer 


       #OUTPUT VARIABLE
pdf                     n                               #options: y,n


       # ALANINE SCANNING

cas                     n                               #options: y,n
	
#MUTATION        P53          112             PHE lig PHE112ALA
#MUTATION        P53          115             LEU lig LEU115ALA
#MUTATION        P53          116             TRP lig TRP116ALA
#MUTATION        P53          119             LEU lig LEU119ALA

使用

使用步骤依据

cd EXAMPLE1 which include Input.dat and the folder including xtc, tpr and ndx
$GMXPBSAHOME/gmxpbsa0.sh and wait few seconds to finish the calculations
when the step 2. is done, type $GMXPBSAHOME/gmxpbsa1.sh and have a cup of coffe :)
when the step 3. is done, type $GMXPBSAHOME/gmxpbsa2.sh and wait few seconds to finish the calculations

“Input.dat”文件参数表Table

玺洛克修订版的GMXPBSA参见：链接

g_MMPBSA

参考资料

入门教程
详细教程
[g_mmpbsa之简单教程]{http://kangsgo.com/18.html}
[g_mmpbsa之参数注释]{http://kangsgo.com/17.html}
g_mmpbsa之参数注释-生信杂谈

若不想编译g_mmpbsa，则在安装好GMX的情况下可直接使用“Includes APBS functionality”的“Pre-compiled executable program”版，不用安装apbs。链接
若使用“To use with external APBS program”的“Pre-compiled executable program”版则可能在计算polar时报错 Segmentation fault。

将“Pre-compiled executable program”版压缩包里的g_mmpbsa和energy2bfac两个文件的所在路径添加到环境变量里即可调用。

使用步骤参照上面的参考资料即可。

否则需要在安装GMX之外安装APBS，然后编译g_mmpbsa.

使用方法

mdp文件及算例包

计算MM：g_mmpbsa -f traj.xtc -s topol.tpr -n index.ndx -mme -mm energy_MM.xvg -decomp -mmcon contrib_MM.dat 依次选择两个group；
计算PB: g_mmpbsa -f traj.xtc -s topol.tpr -n index.ndx -i mmpbsa.mdp -nomme -pbsa -pol polar.xvg -decomp -pcon contrib_pol.dat；
计算SA: g_mmpbsa -f traj.xtc -s topol.tpr -n index.ndx -i mmpbsa.mdp -nomme -pbsa -apol apolar.xvg -decomp -apcon contrib_apol.dat；

g_mmpbsa参数选择：

In addition to Hiqmet’s answer, you must also read the following papers which have exhaustively tested various parameters for MMPBSA (note these have been tested using AMBER but the results should be valid for Gromacs as well,) 1) http://pubs.acs.org/doi/abs/10.1021/ci100275a 2) http://pubs.acs.org/doi/abs/10.1021/ci300385h

前言