异。从 AMBER 7 开始，程序的波恩半径与此前的版本有了明显的不同，而且似乎我们模拟所获得的折叠结果与文献报道的有所差异，后面的分析会告诉我们更多。
模拟的起始阶段，能量上升的很快，这是由于我们对体系进行加热的缘故，而后能量开始下降，这表明蛋白开始折叠。总体来看我们发现体系在50ns的时间尺度里下降了40 kcal/mol 这与文献上所报道的是相符合的。
下一个步骤是定位低能构像这个低能构像就是我们所认为的蛋白质折叠的状态，我们还要以这个构像为基准来做出文献中 Figure 1 B的图。我们可以通过查阅脚本process_mdout.perl生成的
summary.EPTOT文件来找到低能构像。需要主意的是，我们首先需要剥离模拟过程的前50ps，那是升温过程，不算数的。 (summary_EPTOT_stripped.dat.gz). 下面的awk脚本可以迅速找出能量最低的构像。
>gunzip summary_EPTOT_stripped.dat.gz
>cat summary_EPTOT_stripped.dat | awk '{if($2 51.000 -453.0616 52.000 -458.9413 53.000 -466.4404 57.000 -471.1558 61.000 -478.0281 71.000 -479.0216 72.000 -482.4780 74.000 -488.0607 122.000 -494.2667 188.000 -506.0259 224.000 -511.8178 510.000 -519.0695 1278.000 -521.5923 1292.000 -526.4179 3903.000 -528.1458 14896.000 -534.8027 16909.000 -537.3546 17745.000 -539.2624 17759.000 -539.3218 19120.000 -539.6930 19719.000 -539.9273 21316.000 -545.5322 26834.000 -545.9479 42505.000 -547.4913 注意：如果你是用你自己运行的模拟轨迹来运行上述命令，你所获得的输出结果会不同，这是由于你使用的机器不同，所以搜索到的是相空间的不同区域。 如结果所示能量最低的结构出现在42.505 ns其能量为-547.4913 kcal/mol。我们可以使用grep命令来找到这个结构： >grep 42505.000 *.out amber 动力学过程(1)： 1、选择力场 tleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99 2、导入晶体结构 model=loadpdb \ 保存crd和top文件 saveamberparm model polyAT_vac.top polyAT_vac.crd 此时注意电荷是否平衡： 如果缺正电荷 addions model Na+ 0 负离子就用Cl- 选择水箱 solvateoct model TIP3PBOX 8.0 3、保存crd和top文件 saveamberparm model model_wat.top model_wat.crd 4、退出tleap quit 5、保存新的pdb ambpdb -p model_wat.top < model_wat.crd > model2.pdb 6、溶剂环境能量最优化。这一步保持溶质（蛋白）不变，去除溶剂中能量不正常的范德华相互作用。 该步骤的配置文件min1.in如下： --------------------------------------------------------------------- oxytocin: initial minimisation solvent + ions ##说明信息###### &cntrl ##模拟参数起始 imin = 1, ##任务是优化，0 是分子动力学 cut = 10 ##非键相互作用的截断值为10 挨 ntb = 1, ##周期边界条件 0 不采用；1 定容 ；2 定压 maxcyc = 4000, ##优化步数 ntr = 1, ##优化时需要一些约束原子 -ref ncyc = 2000, ##前2000最陡下降，后面步骤共轭梯度 / Hold the protein fixed ##约束说明 500.0 ##作用在肽键上的力 kcal/mol RES 1 9 ##限制的残基序号 同restrain=’:1-9’ END END ------------------------------------------------------------------------------ 任务命令：如果 sander -O -i min1.in -p model_wat.top -c model_wat.crd -o min1.out -r min1.rst –ref model_wat.crd & 7、对蛋白进行优化，min2.in文件将min1.in中的限制原子修改，限制水的位置。 也可以考虑利用restrainmask=’:1-9@CA,N,C’约束蛋白主链上的原子。 sander -O -i min2.in -p model_wat.top -c min1.rst -o min2.out -r min2.rst& 8、整体的优化,去掉限制条件 sander -O -i min3.in -p model_wat.top -c min2.rst -o min3.out -r min3.rst & 9、有限制的分子动力学 第一步分子动力学保持蛋白分子位置不变，但是不是完全固定每个原子，同时缓解蛋白分子周围的水分子，是溶剂环境能量优化。在这个步骤中，我们将主要目的是对特定的原子使用作用力使其能量优化。 Eq1.in 如下： fix protein ,relax H2O &cntrl nstlim=25000, dt=0.002, ntx=1, irest=0, ntpr=500, ntwr=500,ntwx=500, tempi=0.0,temp0=300,ntt=3, gamma_ln=1.0, ntb=1, ntp=0, nrespa=1, cut = 10, ntc=2,ntf=2, NTR=1, / fix protein and HEM 10 RES 1 284 END END ----------------------------- nstlim = #：＃表示计算的步数。dt = 0.002：表示步长，单位为ps，0.002表示2fs。ntx=1 irest=0 默认 ntb = 1：表示分子动力学过程保持体积固定。 imin = 0：表示模拟过程为分子动力学，不是能量最优化。 temp0 = 300：表示最后系统到达并保持的温度，单位为K。 tempi = 100：系统开始时的温度。 ntc=2,ntf=2 忽略氢键 gamma_ln = 1：表示当ntt＝3时的碰撞频率，单位为ps-1（请参考AMBER手册） ntt = 3：温度转变控制，3表示使用兰格氏动力学。 sander -O -i eq1.in -p model_wat.top -c min3.rst -o eq1.out -r eq1.rst -ref min3.rst -x eq1.mdcrd 11整系统分子动力学模拟: eq2 ------------------------------------------------- f2:500ps MD &cntrl imin = 0, irest=1, ntx=5, ntb=2, pres0 = 1.0, ntp=1, taup = 2.0, ntc=2, ntf=2, cut = 10, ntr = 0, ntt = 3, gamma_ln = 1.0, tempi = 300.0 , temp0 = 300.0 nstlim=500000, dt=0.002, ntpr=500, ntwr=500, ntwx=500 / ---------------------------------------------------------- ntb=2：表示分子动力学过程的压力常数。Pres0＝1：引用1个单位的压强。 ntp＝1：表示系统动力学过程各向同性。taup = 2.0：压力缓解时间，单位为ps。 使用以下命令进行MD： sander -O -i eq2.in -p model_wat.top -c eq1.rst -o eq2.out -r eq2.rst -x eq2.mdcrd -ref eq1.rst & 结果处理： 在动力学过程中，将会产生 .rst 和.mdrcd 文件（需要的），由于crd文件很多，可以将其压缩成.gz 文件：gzip filename，将产生一个filename.gz 文件 做出已跑动力学的rmd图，判断是否平衡： ptraj xxx.top < xxx.in xxx.in 内容： trajin eq2.mdcrd.gz trajin eq3.mdcrd.gz trajin eq4.mdcrd.gz trajin eq5.mdcrd.gz rms first mass out xin.rms.dat1 :1-284@CA,C,N time 0.1 #产生一个xin.rms.dat1的文件，整体1-284骨架上的C与N原子 rms first mass out xin.rms.dat2 :88-91,172,201-204,230@CA,C,N time 0.1 #产生一个xin.rms.dat2的文件，保守残基骨架上的C与N原子 ---------------------------------------------------------------------------- 利用xmgrace 作图： xmgrace xin.rms.dat1 xin.rms.dat2 如果需要取随即的点的构型： ptraj xxx.top < xxx.in xxx.in内容： trajin eq7.mdcrd.gz 117 117 #eq7中的117ps 的构型 strip :WAT #冒号前有空格，后没有，注意wat与top的水的大小写一致 strip :Na+ trajout eq7.pdb pdb nobox 产生一个eq7.pdb.117的文件 ------------------------------------------- traj AR.top << EOF > trajin AR.md7.crd 300 500 > center :1-250 > image center familiar > rms first mass out rms.dat :1-250 > average average.rst restrt > average average.pdb pdb > EOF amber 动力学过程(2) 1、利用gaussian * antechamber -i 49.mol2 -fi mol2 -o 49.in -fo gzmat 这样可以生成49.in文件，下载到windows环境，运行gaussian计算这个文件，如果发现计算时间过长或者内存不足计算中断，可以修改文件选择小一些的基组。 you have Gassian output file.out (don’t forget to put some keyword in the Gassian input for the special format used by amber) * -antechamber –i file.out –fi gout –o filr.prep –fo prepi –c resp (this step prepare the “file.prep” for your drug -c bcc all right ) * -parmchk –i file.prep –f prepi –o file.frcmod )this step prepare the file “file.frcmod”for the drug)