《Amber 教程.docx》由会员分享,可在线阅读,更多相关《Amber 教程.docx(14页珍藏版)》请在课桌文档上搜索。
1、(IGB=5)具体可以参考AMBER用户手册的GB模型一举。在论文中,他们没有运用特别的GB模型,这是因为在那时AMBER程序中只右IGB=I这个设定可用。为了使我们的教程尽可能接近文献报道,我们也运用IGB=I的设定。由于1.eap默认的设定就是IGB=I,所以我们无需特地对此作出设定。论文中还声明他们运用了FF99力场,这与我们之前设定的是一样的,但是他们的立场有改进的phi/psi二面角参数,这是对FF99立场中phi/psi二面角参数的一-种校正,可以更好的模拟蛋白质中alpha螺旋的结构。为了更好地重复文献中的工作,我们须要建立一个包含上述修正的参数文件。但是比较麻烦的是,文献中并没
2、有明确给出那些参数做了如何的变更,仅仅给出了-个修正后的parm99.dat文件.出现这种状况的缘由我认为可能是AMBER6本身不带FF99力场,在当时存在许多不同的版本,文献的作者为了让人们了解他们运用的是官方版本的FF99力场所以在论文中展示了Parm99.dat文件。但很不幸,ACS以PDF文件格式给出了这个文件,这使得我们很难干脆运用这个文件。幸运的是,在AMBER8版本中,给出了这个修正的力场,位于下述路径:SAMBERHOMEdatleapparmasfrcmod.mod_phipsi.1.下面我也列出了文件的内容,以备时常:fromSimmcrling,Strockbine,Ro
3、itberg,JACS124:11258,2002.Modifiesparm99.MASSBONDANG1.DIHEDRA1.N-CT-C-N10.700180.000-1.N-CT-C-N11.100180.0002.C-N-CT-C11.0000.0001.NONB如你所见,只有三个二面角参数发生了变更,所以我们只须要打开X1.EaP,读取这个文件,其中的参数就会自动覆盖原有的参数假如现在你已经关闭了周内被截断了,所以在没有运用PME方法的状况K最好要将截断值尽可能设大。须要提示的是,模拟的耗时是与截断值的平方成正比的,(参见教程一的节)所幸我们模拟的体系特别小,足够承受没有截断值(CUt
4、=999)的计算。基于同样的原理我们将rgbmax也设置为999埃,这个参数限制了在计算非犍相互作用过程中列用于计算有效波恩半径的粒子对的最远间距。这个值设定的越大,计算的结果就越好,当然也就须要花费越多的计算时间。考虑到我们面对的体系只有2。个氟基酸残基我们可以把全部的粒子都纳入到有效波恩半径中来,所以我们设定的rgbmax远远大F计算的尺度0.卜面起先运行优化过程:SAMBERHOMEexesander-O-imini.in-omin!.out-pTC5b.prmto-cTC5b.inpcrd-rmini.rstInputFiles:TC5b.prmtop,TC5b.inpcrd,mini
5、.inOutputFiles:mini.out,minl.rst在16个1.3GHzCPU的SGIAltix上这个过程须要3.5秒完成为了直观的比较优化前后的结构,我们生成一个Pdb文件:SAMBERHOME/exe/ambpdb-pTC5b.prmtopminl.pdb将优化前后的两个文件打开(minipdbandTC5binear.pdb)你可以选择任何可用的显示软件,比如VMD起始结构用蓝色显示,优化后的结构用黄色显示。如你所见,优化过程并未造成主链结构太大的变更但是色氨酸和酪氨酸残基发生了比较明显的移动,这些能玷热点其次阶段-10,000步,步长0.5fs(共5ps),结束温度100.
6、0K,温度糊合系数1.0ps第三阶段-10,000步,步长05fs(共5ps),结束温度150.0K,温度期合系数1.0ps第四阶段10,000步,步长0.5fs(共5ps),结束温度200.0K,温度耦合系数1.0ps第五阶段-10,000步,步长05fs供5ps),结束温度250.0K,温度期合系数1.0ps第六阶段-IOQOO步,步长05fs供5ps),结束温度300.0K,温度摘合系数1.0ps第七阶段-40,000步,步长0.5fs供20ps),结束温度325.0K,温度描合系数1.0ps下面是第一阶段的输入文件:heatl.inStage1heatingofTC5b0to50K&c
7、ntrlimin=0,irest=O,ntx=l,nstlim三10000,dt-0.0005,ntc=2,ntf=2,ntt-l,tautp-1.0,tempi=0.0,temp0=50.0,ntpr-50fntwx-50fntb=O,igb=1,cut-999.,rgbmax-999./其他六个阶段的输入文件与之特别接近,只须要变更二下相应的温度就可以了,可以从今处IC载现成的输入文件:(heat2.in,heat3.in,heat4.in,heat5.in,heat6.in,heat7.in).下面是一个运行升温模拟的PBS脚本,你也可以依据你的系统臼己写一个脚本。#PBS-1ncpus
8、=16#PBS-1WalItime=500:00:00sander-O-iheat2.in-pTC5b.prmtop-cheatl.rst-rheat2.rst-oheat2.out-xheat2.mdcrdgzip-9heat2.mdcrdsander-O-iheat3.in-pTC5b.prmtop-cheat2.rst-rheat3.rst-oheat3.out-xheat3.mdcrdgzip-9heat3.mdcrdsander-O-iheat4.in-pTC5b.prmtop-cheat3.rst-rheat4.rst-oheat4.out-xheat4.mdcrdgzip-9he
9、at4.mdcrdsander-O-iheat5.in-pTC5b.prmtop-cheat4.rst-rheat5.rst-oheat5.out-xheat5.mdcrdgzip-9heat5.mdcrdsander-O-iheat6.in-pTC5b.prmtop-cheat5.rst-rheat6.rst-oheat6.out-xheat6.mdcrdgzip-9heat6.mdcrdsander-O-iheat7.in-pTC5b.prmtop-cheat6.rst-rheat7.rst-oheat7.out-xheat7.mdcrdgzip-9heat7.mdcrdmkdirinit
10、ial-heatingcpheat1.outinitial_heatingcheat2,outinitialeatingcpheat3.outinitial-heatingcpheat4,outinitial_heatingcpheat5.outinitial-heatingcpheat6.outinitia1.heatingcpheat7.outinitial-heatingcpheatl.mdcrd.gzinitial_heatingcpheat2.mdcrd.gzinitial_heatingcpheat3,mdcrd.gzinitial_heatingcpheat4.mdcrd.gzi
11、nitial_heatingcpheat5,mdcrd.gzinitial_heatingcpheat6.mdcrd.gzinitial-heatingcpheat7.mdcrd.gzinitial_heatingsetenvAMBERHOME/usrpeople/rcwamber9cd-rcwproductionmpirun-np16SAMBERHOME/exe/sander-O-iequil.in-pTC5b.prmtop- cheat7.rst-requill.rst-oequill.out-xequill.mdcrdgzip-9cquill.mdcrdmpirun-np16SAMBER
12、HOME/exe/sander-O-iequil.in-pTC5b.prmtop- ccquill.rst-rcquil2.rst-ocquil2.out-xequil2.mdcrdgzip-9equil2.mdcrdmpirun-np16SAMBERHOME/exe/sander-O-iequil.in-pTC5b.prmtop- cequi12.rst-requil3.rst-oequil3.out-xequil3.mdcrdgzip-9cquil3.mdcrdmpirun-np16SAMBERHOME/exe/sander-O-iequil.in-pTC5b.prmtop- ccquil
13、3.rst-requil4.rst-oequil4.out-xequil4.mdcrdgzip-9equil4.mdcrdmpirun-np16SAMBERHOMEcxesander-O-iequil.in-pTC5b.prmtop- cequi14.rst-requil5.rst-oequil5.out-xequil5.mdcrdgzip-9cquil5.mdcrdmpirun-np16SAMBERHOME/exe/sander-O-iequil.in-pTC5b.prmtop- cequil5.rst-requil6.rst-oequil6.oul-xequil6.mdcrdgzip-9e
14、qui!6.mdcrdmpirun-n16SAMBERHOME/exe/sander-O-iequil.in-pTC5b.prmtop- cequil6.rst-requil7.rst-oequil7.out-xequil7.mdcrdgzip-9equil7.mdcrdmpirun-np16SAMBERHOME/exe/sander-O-iequil.in-pTC5b.prmtop- cequil7.rst-requil8.rst-oeTransformations-Runningaverages来获得这样的一个曲线势能(红色曲线是WPs平均的势能变更曲线)上图是干脆摘取自文献中的,我们教程
15、中所运用的起始结构其能班要高于文献。这有可能是由于他们运用的是Amber6中的力场,而我们用的是Amber8,这导致两次模拟在广义波恩方程一项上产生差异。从AMBER7起先,程序的波恩半径与此前的版本有了明显的不同,而且好像我们模拟所获得的折强结果与文献报道的有所差异,后面的分析会告知我们更多。模拟的起始阶段,能最上升的很快,这是由于我们对体系进行加热的原因,而后能盘起先卜.降,这表明蛋白起先折狡。总体来看我们发觉体系在50ns的时间尺度里下降了40kcal/mol这与文献上所报道的是相符合的。下一个步骤是定位低能构像这个低能构像就是我们所认为的蛋白质折登的状态,我们还要以这个构像为基准来做出
16、文献中Figure1B的图。我们可以通过查阅脚本process_mdout.perl生成的SUmmaryEPTOT文件来找到低能构像。须要办法的是,我们首先须要剥离模拟过程的前50ps,那是升温过程,不答数的(summary_EPTOT_stripped.dat.gz).下面的awk脚本可以快速找出能R最低的构像。grep42505.000*.outamber动力学过程(1):1、选择力场tleap-S-f$AMBERHOME/dat/leap/cmd/leaprc.ff992、导入晶体结构model=Ioadpdb-sbp_lin.pdb保存crd和top文件saveamberparmmod
17、elpolyAT_vac.toppolyAT_vac.crd此时留意电荷是否平衡:假如缺正电荷addionsmodelNa+0负离子就用Cl-选择水箱SolvateoctmodelTIP3PBOX8.03、保存crd和top文件Saveamberparmmodelmodel_wat.topmodel_wat.crd4、退出tleapquit5、保存新的pdbambpdb-pmodel_wat.topmodel2.pdb6、溶剂环境能量最优化。这一步保持溶质(蛋白)不变,去除溶剂中能量不正常的范德华相互作用。该步骤的配置文件mini.in如下:ifccntrlimin=1,cut=10ntb=1
18、,maxeye=4000,ntr=1,oxytocin:initialminimisationsolvent+ions#说明信息#模拟参数起始# #任务是优化,0是分子动力学# #非键相互作用的截断值为10挨# #周期边界条件0不采纳;1定容;2定压#优化步数#优化时须要一些约束原子-refncyc=2000,#前2000最陡F降,后面步骤共辗梯度Holdtheproteinfixed#约束说明500.0#作用在肽键上的力kcal/molRES19#限制的残基序号同restrain-,:1-9ENDEND任务吩咐:假如sander-O-iminl.in-pmodel_wat.top-cmode
19、l-wat.crd-omin!.out-rminl.rst-refmodel_wat.crd&7、对蛋白进行优化,min2in文件将minl.in中的限制原子修改,限制水的位置。也可以考虑利用restrainmask=:1-9CA,N,C约束蛋白主链上的原子。sander-O-imin2.in-pmode1.Wat.top-cminl.rst-omin2.out-rmin2.rst&8、整体的优化,去掉限制条件sander-O-imin3.in-pmodel_wat.top-cmin2.rst-omin3.out-rmin3.rst&9、有限制的分子动力学第一步分子动力学保持蛋白分子位置不变,
20、但是不是完全固定每个原子,同时缓解蛋白分了四冏的水分子,是溶剂环境能量优化.在这个步腺中,我们将主要目的是对特定的原子运用作用力使其能Ia优化。Eq1.in如下:fixprotein,relaxH2O&cntrlnstlim=25000,dt=0.002,ntx=l,irest=O,ntpr=500,ntwr=500,ntwx=500,temi=0.0,temp0=300fntt=3fgammajn=l.,ntb=1,ntp=O,nrespa=1,cut=10,ntc=2,ntf=2,NTR=1,/fixproteinandHEM10RES1284ENDENDnstlim=#:#表示计算的步数
21、。dt=0.002:表示步长,总位为ps,0.002表示2fs。ntxlirestO默认ntb-1:表示分子动力学过程保持体积固定。imin=0:表示模拟过程为分子动力学,不是能量最优化。temp=300:表示最终系统到达并保持的温度,胞位为K。tempi=100:系统起先时的温度。ntc=2,ntf=2忽视氢健gammajn=1:表示当ntt=3时的碰撞频率,单位为PS-I(请参考AMBER手册)ntt=3:温度转变限制,3表示运用兰格氏动力学。sander-O-ieql.in-pmodel_wat.top-cmin3.rst-oeq1.out-reql.rst-refmin3.rst-xe
22、ql.mdcrd11整系统分子动力学模拟:eq2f2:500psMD85cntrlimin=0,irest=1,ntx=5,ntb=2,pres=1.0,ntp=l,taup=2.0,ntc=2,ntf=2,cut-10,ntr-0,ntt-3,gamma_ln-1.0,tempi-300.0,temp-300.0nstlim=500000,dt=0.002,ntpr=500,ntwr=500,ntwx=50011tb=2:表示分子动力学过程的压力常数。PresO=I:引用1个帽位的压强。ntp=l:表示系统动力学过程各向同性。taup-2.0:压力缓解时间,单位为pso运用以下吩咐进行MD:
23、sander-O-ieq2.in-pmodel_wat.top-ceql.rst-oeq2.out-req2.rst-xeq2,mdcrd-refeql.rst&结果处理:在动力学过程中,将会产生rst和.mdrcd文件(须要的),由于Crd文件许多,可以将其压缩成gz文件:gzip,将产生一个文件做出已跑动力学的md图,推断是否平衡:ptrajXXX.topxxx.inxxx.in内容:rmsfirstmassoutxin.rms.dat1:1-284CA,C,Ntime0.1#产生一个的文件,整体1-284骨架上的C与N原子rmsfirstmassoutxin.rms.dat2:88-91
24、,172,201-204,230CA,C,Ntime0.1#产生一个的文件,保守残基骨架上的C与N原子利用xmgrace作图:假如须要取随即的点的构型:ptrajXXX.topxxx.inxxx.in内容:trajineq7.mdcrd.gz117117#eq7中的117ps的构型strip:WAT#冒号前有空格,后没有,留意Wat与top的水的大小写一样strip:Na+trajouteq7.pdbpdbnobox产生一个的文件trajAR.top trajinAR.md7.crd300500 center:1-250 imagecenterfamiliar rmsfirstmassoutr
25、ms.dat:1-250 averageaverage.rstrestrt averageaverage.pdbdb EOFamber动力学过程(2)1、利用gaussian* antechamber-i49.mol2-fimol2-o49.in-fogzmat这样可以生成49.in文件,卜.载到WindQWS环境,运行gaussian计算这个文件,假如发觉计算时间过长或者内存不足计算中断,可以修改文件选择小一些的基组。youhaveGassianoutput(don,tforgettoputsomekeywordintheGassianinputforthespecialformatused
26、byamber)* -antechamber-i-figout-ofilr.prcp-foprepi-cresp(thissteppreparetheunforyourdrug-cbccallright)* -parmchk-i-fprepi-o)thissteppreparethefileforthedrug)2、干脆用mol2转也可以干脆把sybyl的mol2做参数* -antechamber-i-fimol2-ofilr.prep-foprepi-cresp/-cbcc* -parmchk-i-fprepi-onowyouhave2files:andforthedrug.留意修改中文件名
27、。要与Pdb中的配体一样。(viprep=假设是abc)对应pdb与prep文件中的原子名称,,原子依次要一样。prep在XIeaP中打开usingtleapfortheproteinordrugor/andcomplex* xleap-s-fIeaprc.ff99* mod=Ioadambcrparams-Ioadamberprcpeditabcpdb可以用sybyl打开。用vi修改pdb-sourceleaprc.gaff* R1.-Ioadpdb* SolvateoctR1.TIP3PB0X10-AddionsR1.Cl-0(O:zeroforneutralcharge,Cl-canbe
28、replacedbyNa+)-SaveamberparmR1.(savetopologyandcoord.Withwater)amber动力学,博大精深!做参数还需好好学习,记录一下过程,以免遗忘。一般来说,小分子的参数是最麻烦的,首先具备的是mol2文件,做参数主要是电荷问题,-般有bcc和resp电荷,个人觉得选用resp电荷可能好一点,曾经用bcc发觉优化做不下去。1、用resp电荷前,需用gaussian计算电荷分布,其中gaussian计算的关键词为:%mem=400Mb%nroc=1#PHF6-31G*pop=mkiop(633=2)iop(642=6)2、生成aaa.out文件后
29、,利用antechamber制作prep文件antechamber-iaaa.out-figout-oaaa.prcp-foprcpi-cresp#输入文件名、文件类型、输出文件名、输出文件类型,电荷类型3、产生PreP文件后,就是检杳参数,armchk-iaaa.prep-fiprepi-oaaa.frcmod在aaa.frcmod文件中是一些非标准的犍长,键角,二面角等信息。4、修改PreP文件中,小分子的名字,默认的都是MO1.,改成Pdb中曳合物中的小分子一样的名字(设为ABC5、利用XIeaP修改原子名称C由于mol2文件转成的prep文件的原子名称往往与pdb中小分子的名称不一样,
30、这是可以用XleaP中的editABC0(假设Pdb中小分子的名字叫ABC).在图形中选择显示原子名称,这样这就是amber将会读的原子名称,是PreP文件中的。可以用Sybyl打开Pdb,现在原子id,对应位置,利用Vi编辑器,打开Pdb,修改pdb中的原子名称与Xleap中的一样。6、修改工作一般很细小,确定要当心,完成后保存退出。为了以防万一,备份一下Pdbo然后Vi打开Pdb,将原子键连的信息删掉。小分子与蛋白间有必要可以加上TER,断开各个小分子,这是Pdb格式中的一个通用分段方式,具体可以查看Pdb的文件格式信息。7、完成后就可以SUS=IOadPdbabC.pdb导入pdb,这样一般不会有新增加原子了。在做动力学的时候有时候出现以卜错误信息:Coordinateresetting(SHAKE)cannotbeaccomplished,deviationistoolarge这时候可以用ptraj中的checkoverlap检查原子间是否有很近的原子ptrajEOFYourParm.toptrajinYourCRD.crdcheckoverlapgoEOF