本示例来自gromacs网站上的示例:gromacs tutorial for drug-enzyme complex.有些步骤按其上面所述服务器无法运行,自己稍微做了修改。
1.从http://www./pdb/home/home.do下载1az8蛋白的pdb文件。
2.用UltraEdit打开1az8.pdb文件,将HETATM一段(除去其中的水)拷入到The Dundee PRODRG2 Server,在“Chirality”,“Full charges”,“Energy Minimization”三个选项分别选为“Yes”,“Yes”,“No”,点击“Run PRODRG”生成抑制剂文件,解压,将DRGGMX.ITP文件重命名为drg.itp文件;
3. 编辑蛋白质坐标文件trp_minv.pdb导入工作站;
4. 在进行任何分子动力学模拟之前,必须建立分子的拓扑文件。Gromacs的分子拓扑文件是用pdb2gmx命令生成,文件后缀名为 .top。pdb2gmx的输入唯一文件是分子的PDB文件。把pdb分子文件转化成gromacs独特的gro分子结构文件类型,同时产生分子拓扑文件,使用命令:
pdb2gmx -ignh -ff gmx -f trp_minv.pdb -o trp2.pdb -p trp.top -water spce
解释:
-f : 指定你的坐标文件,可以是pdb、gro、tpr等等包含有分子坐标的文件; -o : 输出文件,也就是处理过的分子坐标文件,同样可以是pdb、gro、g96等文件类型; -p :输出拓扑文件。pdb2gmx读入力场文件,根据坐标文件建立分子系统的拓扑; -water :指定使用的水模型,使用pdb2gmx的时候最好加这个参数,不然后面会吃苦头。它会提前在拓扑文件中添加水分子模型文件; -ff :指定力场文件(下文讨论),也可以不用这个参数,再自行选择; -ignh : 舍弃分子文件中的H原子,因为H原子命名规则多,有的力场不认; 5. 打开生成的trp2.pdb文件,将DRGPOH.PDB文件复制到trp2.pdb文件的结尾,将残基数列改为224,将原子序列也作相应调整,保存(原文说接在trp.gro文件的ASN223后面,但现在还没有生成这个文件,感觉应该是trp2.pdb);
6. 编辑tro.top文件:在力场部分添加 #include “drg.itp”,在分子部分添加 DRG 1;
7. 真空条件下进行分子模拟输出结果误差较大,所以模拟之前,必须给分子添加水环境。使用genbox命令为gromacs模拟分子添加水环境。
使用editconf为体系添加一个盒子:
editconf -bt cubic -f trp2.pdb -o trp2.pdb -d 0.9
解释:
-f : 指定你的坐标文件。 -o : 输出文件,即放进盒子里面的分子系统。 -bt : 盒子类型,有正方型,长方形,八面型等等,看个人需要跟癖好啦。 -d : 分子离盒子表面的最短距离。这个跟-bt一起使用,基本就足够了;如果蛋白在模拟过程尺寸变化很大,那就用-box吧。 使用genbox向盒子里添加水分子:
genbox -cp trp2.pdb -cs spc216.gro -o trp_b4ion.pdb -p trp.top
解释:
-cp :带盒子参数的分子坐标文件,也就是editconf的输出文件; -cs :添加的水分子模型,如spc216、spce、tip3p、tip4p等,关于各个模型的区别,请参考scholar google; -o :输出坐标文件,就是添加水分子之后的分子坐标文件,默认是.gro文件,但是也可以输出其他文件格式,如pdb; -p :系统拓扑文件,genbox会往里面写入添加水分子的个数。
8. 编辑MDP文件。在gromacs里,控制分子动力学模拟的参数文件称mdp(Molecular Dynamic Parameters)文件。该文件中定义分子动力学模拟的各种行为,内容如原文。每一项的内容为:
title:文件头,你想定义你模拟文件的名字,可以随意定义。
cpp:预处理器,于C/C++的预处理器一样,默认为(/lib/cpp)
define;预定义,设有默认预定义。可以使用如何模拟的预定义方法控制模拟进程,gromacs目前手册可供选择的有以下两种(其实还有其他,可以查看邮件列表): " -DFLEX_SPC ":即让gromacs为SPC水模型引入软性性质。这样可以在进行能量最优话是使共轭梯度法产生作用,同时也可以更进一步进行最速下降法进行能量最优话。 " -DPOSRE ":让gromacs把posre.itp文件包含到拓扑文件中,这样做是为了进行坐标限制性动力学模拟。
dt: 模拟步长,默认为0.001ps(只对md,sd和bd有效)。
nsteps: 整合算法的最大模拟步数。
其余的选项包括控制时间、精度、温度、坐标什么的。
9. 在给蛋白质添加水环境之后,系统中蛋白质本身已经带了静电量,需要要给系统加几个带相反电量的金属离子,使系统处于电中性。使用命令:
grompp -f em.mdp -c trp_b4ion.pdb -p trp.top -o trp_b4ion.tpr
形成参数,溶液呈电中性,需要三个负电荷中和。使用命令:
genion -s trp_b4ion.tpr -o trp_b4em.pdb -nname Cl -nn 8 -g trp_ion.log
为体系添加金属离子,选择12(SOL)。重新编辑trp.top文件。在“molecules”一栏,添加“CL 8”,然后在SOL数量中减8。重新保存;
解释:
-s: 指定系统tpr文件。 -p: 指定系统拓扑文件,在往系统中添加金属离子时,genion会往拓扑文件最后的分子类型中写入添加的离子数,并修改拓扑文件中系统原子数。 -o: 指定输出文件,genion的输出是pdb文件或者gro等结构文件。也就是说你产生这个文件之后,还要再用这个文件产生tpr文件。 -np/-nn: 带正/负电金属粒子的数目。 -pn/-nn 指定正负金属离子的名字,比如" NA+ "或者" CL- "。
10. 添加分子系统过程中,系统有很多原子距离太近,局部能量太高。这些相互距离太近的原子多是由于genbox程序产生的,溶剂分子与蛋白分子之间存在不稳定的高能量。如果现在开始模拟计算,而不进行能量最优化,系统将可能很不稳定。去除这些局部高能量的办法是对系统进行能量最优化。能量最优化过程是改变系统中局部高能量的原子的位置,降低这些点的能量。在进行能量最优化之前,先用GROMACS的预处理程序grompp处理所有输入文件。grompp预处理拓扑文件(.top),坐标文件(.gro)和一个参数文件(.mdp),然后输出一个二进制拓扑文件(.tpr)。这个二进制文件包含所有模拟需要的信息,利用这个文件即可以进行能量最优化和动力学模拟。使用命令:
grompp -f em.mdp -c trp_b4em.pdb -p trp.top -o em.tpr
11.使用mdrun命令进行能量最优化:
mdrun -v -deffnm em
在进行能量最优化过程mdrun程序的输出文中,从左到右第一个数字是模拟步数,第二个数字是计算步长,第三个数字是系统能量。可以看到系统能量从一个很高的值迅速降低,最后稳定在一个大负值。
12. 能量最优化之后,首先保持蛋白质不动,对蛋白质周围水环境进行动力学模拟,该过程称为位置限制性分子动力学(position restrained MD)。
位置限制分子动力学保持蛋白质位置不变,对溶剂分子进行平衡计算,可以使溶剂分子填补空间空洞。这个可能存在的空洞是genbox程序产生的。
首先对输入文件进行预处理得到二进制拓扑文件,这些输入文件就是能量最优化得到的输出文件。然后再写一个参数配置文件和索引文件。
默认情况下,程序对模拟系统是分部分的。我们利用两个部分进行位置限制性模拟:Protein和 SOL部分,分别表示蛋白质和溶剂。在这个过程中,即保持protein位置不变。
参数文件(.mdp)包含了位置限制性模拟的所有参数,包括步长,步数,温度等等。同时参数文件告诉GROMACS模拟的类型,如能量最优化、位置限制性或者是分子动力学模拟。需对原文提供的文件做一下修改如下:
tau_t = 0.1 0.1 0.1 0.1 tc-grps = protein sol DRG Cl ref_t = 300 300 300 300
后面的就完全跟文献上的步骤一样了,
13.将轨迹文件压缩成一个可用ngmx分析的*.xtc文件中,命令为:
trjconv -f md.trr -o md.xtc
14.使用轨迹观察器绘制分子的3-D结构。
后面的都是分析了,现在还不太懂,以后再续,呵呵
|