当前位置:首页 > amber中生成小分子模板
Created a new atom named: O35 within residue: .R Created a new atom named: N34 within residue: .R
如果屏幕仅仅显示添加原子,这说明输入的PDB文件中缺失了部分重原子,leap根据模板自动补全了这些缺失的原子,这种情况
不会影响进一步的计算
Added missing heavy atom: .R.A Added missing heavy atom: .R.A Added missing heavy atom: .R.A Added missing heavy atom: .R.A
根据体系的具体情况,还可能要将成对的CYX残基中的二硫键相连,有时候还会链接其他原子,比如将糖基链接在氨基酸残基上
,用bond命令可以完成,命令用法如下: >bond comp.35.SG comp.179.SG
其中comp是刚才读入的分子名称,35和179是残基序号,SG是CYX残基模板中硫原子的名称,用comp.35.SG这样的语法就可以定位一个原子.果希望进行考虑溶剂效应的动力学模拟,可能还需要为体系加上水,加水有很多种命令,这里只列举一个:
>solvatebox comp TIP3PBOX 10.0
solvatebox命令是说要加上一个方形的周期水箱,comp指要加水的分子,TIP3PBOX是选择的水模板名称,10.0是水箱子的半径.的体系总电荷不为0,为了模拟稳定,需要加入抗衡离子,系统会自动计算体系的静电场分布,在合适的位置加上离子,命令如下:
>addions comp Na+ 0
意思是用钠离子把体系总电荷补平,还可以选择其他库里面有的离子。 完成这一步之后就可以输出拓扑文件和坐标文件了,但是为了方便起见,在运行最后一步之前要先把leap里加工好的分子单独存成一个库文件,以后可以直接调用这个库文件,免得重复上面的操作:
>saveoff comp 1taj.off
这样就生成了一个off文件在那里,下面输出拓扑文件和坐标文件 >saveamberparm comp 1t4j.prmtop 1t4j.inpcrd Checking Unit. Building topology. Building atom parameters. Building bond parameters. Building angle parameters. Building proper torsion parameters. Building improper torsion parameters. total 1 improper torsion applied Building H-Bond parameters.
Not Marking per-residue atom chain types. Marking per-residue atom chain types. (Residues lacking connect0/connect1 - these don't have chain types marked: res total affected CMET 1 )
(no restraints) >quit
现在准备好拓扑文件和坐标文件,接下来就可以开始能量优化和动力学模拟了。如果愿意的话还可以用ambpdb这个命令生成一个pdb文件,直观地看一看生成了什么东西,命令如下:
$ ambpdb -p 1t4j.prmtop <1t4j.inpcrd> kankan.pdb
| New format PARM file being parsed.
| Version = 1.000 Date = 09/08/06 Time = 16:36:09 $
可以下载之后用看图软件欣赏,如果加了溶剂盒子的话,看的时候可要小心,会恨吓人的样子。
第四步:能量优化
用amber进行分子动力学模拟需要坐标和拓扑文件,这在上一步已经完成了,分别生成了1t4j.prmtop 和1t4j.inpcrd两个文件,下面一步就是动力学模拟之前的能量优化了。
由于我们进行的起始结构来自晶体结构或者同源模建,所以在分子内部存在着一定的张力,能量优化就是在真正的动力学之前,释放这些张力,如果没有这个步骤,在动力学模拟开始之后,整个分子可能会因此散架。
能量优化由sander模块完成,运行sander至少需要三个输入文件,分别是分子的拓扑文件,坐标文件以及sander的控制文件。
现在分子的拓扑文件和坐标文件已经完成,需要建立输入文件,min_1.in Initial minimisation of our structures &cntrl
imin=1, maxcyc=4000, ncyc=2000, cut=10, ntb=1,ntr=1, restraint_wt=0.5 restraintmask=':1-283' /
文件首行是说明,说明这项任务的基本情况; &cntrl与/之间的部分是模拟的参数
其中imin=1表示任务是能量优化,maxcyc=4000表示能量优化共进行4000步,ncyc=2000表示在整个能量优化的4000步中,前 2000步采用最陡下降法,在2000步之后转换为共轭梯度法,如果模拟的时候不希望进行方法的转换,可
以再加入另一个关键词NTMIN,如果NTMIN =0则全程使用共轭梯度法,NTMIN=2则全程使用最陡下降法,此外还有=3和=4的选项,分别是xmin法和lmod法,具体情况可以看手册。
第二行的cut=10表示非键相互作用的截断值,单位是埃, ntb=1表示使用周期边界条件,这个选项要和前面生成的拓扑文件坐标文件相匹配,如果前面加溶剂时候用的是盒子水,就设置ntb=1,如果加的是层水,那就应该选择ntb=0;ntr=1表示在能量优化的过程中要约束一些原子,约束的是哪些原子呢?后面有。
第三行和第四行都是关于约束原子的信息,restraint_wt=0.5限定了约束的力常数,在这里约束原子就是把原子用一根弹簧拉在固定的位置上,一旦原子偏离固定的位置,系统就会给他施加一个回复力,偏离的越远,回复力越大,回复力就是由这个力常数决定的,单位是Kcal/(mol*A)。 restraintmask=':1-283'表示约束的是从1到283号残基,在这个分子中,1-283号残基是蛋白质上的氨基酸残基,从284号开始,就都是水了,所以这个命令的意思就是,约束整个蛋白质,自由地优化溶剂分子。因为溶剂分子是前面的tleap自动给加上的,
所以一定要额外多关注一些。 下面运行sander来执行这个优化:
$ sander -O -i min_1.in -p 1t4j.prmtop -c 1t4j.i npcrd -ref 1t4j.inpcrd -r 1t4j_min1.rst -o 1t4j_min1.out
命令中,-O表示覆盖所有同名文件,-i min_1.in表示sander的控制文件是min_1.in,-p 1t4j.prmtop表示分子的拓扑文件,-c 1t4j.inpcrd表示坐标文件,-ref 1t4j.inpcrd是参考坐标文件,只有在控制文件中出现关键词ntr=1的时候才需要给sander指定-ref文件,这是约束原子的参考坐标,- ref 1t4j.inpcrd就是说以1t4j.inpcrd中的坐标为准进行约束原子的优化。以上这四个是输入文件。-r 1t4j_min1.rst表示经过模拟之后新的原子坐标会输出到1t4j_min1.rst文件中,-o 1t4j_min1.out 则表示优化过程中的相关信息都会写入到1t4j_min1.out文件中。运行起这个命令之后,等拿到 1t4j_min1.rst文件就意味着对溶剂的优化已经差不多了,显然下面还需要对蛋白质本身进行优化,这个优化还要分两步进行,控制文件分别是min_2.in 和min_3.in
min_2.in
共分享92篇相关文档