研究案例——一种稳定蛋白质的全部原子结构预测和折叠模拟
这段教程展示的是一个研究实例，像您演示如何重现下述文章中的研究工作：
Simmerling, C., Strockbine, B., Roitberg, A.E., J. Am. Chem. Soc., 2002, 124, 11258-11259 (http://dx.doi.org/10.1021/ja0273851)
我们建议您在开始本教程前首先阅读上述文章，获得该蛋白的氨基酸序列及其他有用信息。
警告1: 本教程中的一些计算耗时很长，我使用了由16个1.3GHz cup的SGI Altix进行了27小时计算才完成整个工作，因此如果您没有足够的计算能力，我强烈建议您在重复本教程的过程中使用我为您提供的out文件，以使得您能够流畅地完成整个教程。
警告2: 如果您重复本教程，我们并不能保证您能够精确地重现我的计算结果，在计算过程中，不同结构的计算机会产生不同的近似误差，从而使得计算过程搜索的是相空间的不同部位，但是模拟的平均结果是大致相同的。另外，尽管您完全重复了本教程也有可能无法获得论文中给出的结果，而且即便是我们自己也无法保证论文中的结果能够重现，这可能是因为我模拟的时间不够长，获取的仅仅是一个局部最小点，但是尽管如此，本教程的工作还是展示了蛋白折叠中一些有趣的行为。 背景
这篇论文应用AMBER FF99力场和经典的全原子动力学对一个肽的折叠过程进行了模拟。模拟的对象\是一个由20个氨基酸构成的小肽，华盛顿大学的 Andersen已经对这个蛋白做过了结构优化，它是现在已知最小的能够显示两种不同折叠状态的蛋白，而且这个蛋白在室温下可以稳定存在。该蛋白的小身量使得它成为模拟蛋白质折叠的绝嘉对象。当最早的关于这个蛋白的折叠的计算结果出炉时，对这个蛋白结构的实验测定还没有完成，所以整个模拟过程是在没有实验数据作为指导的情况下完成的。当蛋白的结构经由实验手段测定之后，人们惊喜地发现，计算机模拟的结果与实验测定的数值之间的RMSD值仅为1.4A。考虑到整个模拟过程是从蛋白的一级结构开始并且完全没有同源蛋白作为参考，这样的一个计算结果是非常精确的。
本教程中，我们试图重复论文中的结果，计算的设定都与论文非常接近，只是由于计算能力的限制，在教程中我们只进行一个50ns级的模拟。这已经足够重见蛋白质折叠的结果了。在这里必须提醒的是，由于模拟过程的长度所限，在不同的计算机，或在处理器数量不同的情况下，计算的结果将会是不同的。这是由分子动力学模拟的方法决定的，实施过程的细微变化或者浮点计算中舍入的变化都意味着由不同的计算机进行采样的动力学轨迹会随着时间的流逝发生不可预知的分化。这并非误差或者程序的bug，也并不意味着某一个模拟过程比其他的过程更合理。这仅仅意味着不同的模拟过程搜索的是相空间的不同区域，如果我们平均一下模拟的结果，或者运行更长时间的动力学过程，我们会在不同的机器上得到完全相同的结果，他们之间仅仅在过程上有所不同。因而我们说在本教程中我们很难精确的再现论文中的结果，但是我们试图重新创造那个重要的结果，即用AMBER程序来预测一个20氨基酸的小蛋白的空间结构是可以完成的。 那么记住这一点，让我们开始吧 第一步：构建起始结构
在以往的教程中，我们要么有一个可用的晶体结构，要么可以通过程序生成一个已经初步优化的结构。而在这个教程中我们要用的结构太复杂，没法通过手画的办法完成，同时我们也没有一个可用的PDB结构，因此我们就需要构建一个线形的肽链，非常幸运的是，在LEAP中有一个命令可以完成这个工作，就是 sequence。
蛋白的一级结构序列在所列论文中可以查到，如下所示: NLYIQWLKDGGPSSGRPPPS
这是用单字母符号显示的蛋白质一级结构序列，在Leap中使用之前我们需要将其转换成标准的三字母表示下面的表格给出了单字母表示和三字母表示之间的转换关系：
单字母与三字母的转换 conversion 甘氨酸 (Gly) 脯氨酸 (Pro) G P 丙氨酸(Ala) A 缬氨酸(Val) V 亮氨酸 (Leu) L 异亮氨酸 (Ile) I 蛋氨酸 (Met) M 半胱氨酸 (Cys) C 苯丙氨酸 (Phe) F 酪氨酸 (Tyr) Y W 色氨酸 (Trp) H 组氨酸 (His) K 赖氨酸 (Lys) R 精氨酸 (Arg) Q 谷氨酸盐 (Glu) N 天冬酰氨 (Asn) E 谷氨酸 (Glu) D S 天冬氨酸 (Asp) T 丝氨酸 (Ser) 苏氨酸 (Thr) 那么上述序列可以转写为： ASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO SER 但是这还没有结束，LEaP 不能自动识别序列的两端，所以我们必须手工为这个序列标定N末端和C末端，标定的方法就是在N末端氨基酸前方加上N，C末端氨基酸前方加上字母C。最终在 LEaP中使用的序列如下： NASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO CSER 下面启动xleap并调用ff99力场: $AMBERHOME/exe/xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99 (使用xLeap的时候一定要记住要关闭Num Lock键！否则工具栏会无法使用) 下面使用 sequence 命令来建立蛋白的起始结构 (如需了解 sequence 命令的详细情况可以在Leap中键入: help sequence). 注意： 为了版面设计的需要，下面将命名分为三行显示，实际上您必须将所有内容在一行内输入，其间不能回车。 >TC5b = sequence { NASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO CSER } 我们需要的起始结构就放在 TC5b中我们可以使用 edit命令来观察这个结构。 >edit TC5b 现在我们获得了一个线形的蛋白质序列作为起始结构，但是在这个起始结构中很多原子是相互抵触的，所以在进行分子动力学模拟之前我们要对这个结构首先进行短时间的优化。我们暂时将Unit中的这个结构存成一个.lib文件，这样在之后的操作中，我们只要调用这个lib就可以简单地取出起始结构，同时我们还要将这个结构存成一个PDB文件，以便直观地进行观察。
>saveoff TC5b TC5b_linear.lib
>savepdb TC5b TC5b_linear.pdb
(TC5b_linear.lib, TC5b_linear.pdb) 第二步：创建prmtop和inpcrd文件
我们已经有了起始结构，下一步的工作是创建prmtop以及inpcrd文件。在进行这一步之前我们需要首先确认我们使用的参数和文献中报道的是一样的，在论文的第三段讲到：
We initiated our simulations using only the trpcage TC5b2 amino acid sequence
(N20LYIQWLKDGGPSSGRPPPS39), with an extended initial conformation built by the LEaP module of AMBER version 6.0.4 All molecular dynamics (MD) simulations were fully unrestrained and carried out in the canonical ensemble using the SANDER module, which we modified to improve performance on the
5
Linux/Intel PC cluster that was used for all calculations. The ff99 force field was employed, with the
6
exception of [phi/psi] dihedral parameters which were refit (see Supporting Information) to improve
7
agreement with ab initio relative energies of alanine tetrapeptide conformations. Parameters were not fit
8
to data for the trpcage. Solvation effects were incorporated using the Generalized Born model, as
9
implemented in AMBER.
文献显示，他们首先建立了一个线形的起始结构，这一步我们已经完成了，之后他们运行了没有限制的恒温动力学模拟过程（即正则系综中的模拟）在动力学过程中他们使用了广义波恩近似来模拟溶剂效应的影响。AMBER程序可以支持很多不同的广义波恩模型，在这些模型中最先进的是由A. Onufriev, D. Bashford 和 D.A. Case等人开发的改良GB模型，这个GB模型使用了模型II的半径
(IGB=5) 具体可以参考AMBER用户手册的GB模型一章。在论文中，他们没有使用特殊的GB模型，这是因为在那时AMBER程序中只有IGB=1这个设定可用。为了使我们的教程尽可能接近文献报道，我们也使用IGB=1的设定。由于Leap默认的设定就是IGB=1，所以我们无需专门对此作出设定。 论文中还声明他们使用了FF99力场，这与我们之前设定的是一样的，但是他们的立场有改进的 phi/psi二面角参数，这是对FF99立场中phi/psi二面角参数的一种校正，可以更好的模拟蛋白质中 alpha螺旋的结构。为了更好地重复文献中的工作，我们需要建立一个包含上述修正的参数文件。但是比较麻烦的是，文献中并没有明确给出那些参数做了如何的改变，仅仅给出了一个修正后的
parm99.dat文件。出现这种情况的原因我认为可能是 AMBER6本身不带FF99力场，在当时存在很多不同的版本，文献的作者为了让人们了解他们使用的是官方版本的FF99力场所以在论文中展示了 parm99.dat文件。但很不幸，ACS以PDF文件格式给出了这个文件，这使得我们很难直接使用这个文件。幸运的是，在AMBER8版本中，给出了这个修正的力场，位于下述路径：
$AMBERHOME/dat/leap/parm/ as frcmod.mod_phipsi.1. 下面我也列出了文件的内容，以备不时：
frcmod.mod_phipsi.1 from Simmerling, Strockbine, Roitberg, JACS 124:11258, 2002. Modifies parm99. MASS BOND ANGL DIHEDRAL N -CT-C -N 1 0.700 180.000 -1. N -CT-C -N 1 1.100 180.000 2. C -N -CT-C 1 1.000 0.000 1. NONB 如你所见，只有三个二面角参数发生了变化，所以我们只需要打开xLEaP，读取这个文件，其中的参数就会自动覆盖原有的参数。如果现在你已经关闭了xLEaP，可以重新打开并调入蛋白质结构： $AMBERHOME/exe/xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99 >loadoff TC5b_linear.lib 然后调入修正的二面角参数： >loadamberparams frcmod.mod_phipsi.1 下面就可以存储我们的 prmtop和 inpcrd文件: >saveamberparm TC5b TC5b.prmtop TC5b.inpcrd 下面是生成的输入文件 TC5b.prmtop, TC5b.inpcrd 第三步：预优化蛋白质. 在运行分子动力学模拟之前我们必须对起始结构进行短时间的优化，这样的话我们的体系就不会因为局部能量的聚集而在动力学过程中出现问题。结构优化的过程会整理整个分子结构，重新修整氢的位置，这样的过程之后我们的动力学模拟会比较稳定 。 下面是结构优化的输入文件: min1.in Stage 1 - minimisation of TC5b &cntrl imin=1, maxcyc=1000, ncyc=500, cut=999., rgbmax=999.,igb=1, ntb=0, ntpr=100 / 我们总共运行1000步优化过程，其中500步为最陡下降法(ncyc=500)，然后紧跟500步共轭梯度法(maxcyc-ncyc)。这样的设置已经足够充分地释放聚集在起始结构中的能量。需要提醒的是我在输入文件中设置了非常大的截断值(cut=999. angstroms)，这样设置是因为我们使用了非周期模拟(ntb=0) ，故而我们没有使用PME方法，也就不会出现长程的静电相互作用。如果使用了 PME，推荐的截断值是8埃，在这样一个范围内实际上模拟的主要是范德华相互作用。 但是如果不使用PME而设置截断值的话，范德华相互作用和静电相互作用都在截断值的范围内被截断了，所以在没有使用PME方法的状况下，最好要将截断值尽可能设大。需要提醒的是，模拟的耗时是与截断值的平方成正比的，（参见 教程一 的 5.1.2节）所幸我们模拟的体系非常小，足够承受没有截断值(cut=999)的计算。基于同样的原理我们将rgbmax也设置为 999埃，这个参数控制了在计算非键相互作用过程中列用于计算有效波恩半径的粒子对的最远间距。这个值设定的越大，计算的结果就越好，当然也就需要花费越多的计算时间。考虑到我们面对的体系只有20个氨基酸残基我们可以把所有的粒子都纳入到有效波恩半径中来，所以我们设定的rgbmax远远大于计算的尺度。. 下面开始运行优化过程: $AMBERHOME/exe/sander -O -i min1.in -o min1.out -p TC5b.prmtop -c TC5b.inpcrd -r min1.rst
Input Files: TC5b.prmtop, TC5b.inpcrd, min1.in Output Files: min1.out, min1.rst
在16个1.3GHzCPU的 SGI Altix上这个过程需要3.5秒完成 为了直观的比较优化前后的结构，我们生成一个pdb文件:
$AMBERHOME/exe/ambpdb -p TC5b.prmtop < min1.rst > min1.pdb
将优化前后的两个文件打开(min1.pdb and TC5b_linear.pdb)你可以选择任何可用的显示软件，比如VMD
起始结构用蓝色显示，优化后的结构用黄色显示。如你所见，优化过程并未造成主链结构太大的变化但是色氨酸和酪氨酸残基发生了比较明显的移动，这些能量热点集中的区域有可能在我们开始分子动力学模拟之后带来麻烦，如果你不相信，可以用未经优化的结构跑一个动力学过程看看，肯定飞！ 第四步：体系加热.
接下来我们要在这个体系中正式开始分子动力学模拟，首先我们要分7步花费50ps时间对体系进行升温模拟。将升温过程分为7步完成可以在每一步升温之后维持一段时间，以免一次升温造成体系能量聚集最终跑飞，另外一种可行的方法是对升温过程加一个权重限制。您可以参阅AMBER用户手册以获取更多信息。一般而言我们升温的最终目标是室温即300K但是为了重复文献的运算我们选择325K：
MD simulations of 100 ns were performed at 300 K, but all were kinetically trapped on this time scale, showing strong dependence on initial conditions and failing to converge to similar conformational ensembles. We therefore increased the temperature to 325 K.
文献认为必须将体系加温到325K进行模拟，否则有可能使模拟的结果最终落入局部最小点，所以我们也做同样的设定。但是你必须牢记更高的模拟温度会导致体系中各化学键发生更加显著的振动，这意味着如果你打算做一个600K，以2fs为步长的动力学模拟，你就要考虑一个应用了shaken的300k效果会与之相同，但600K的模拟却要临步长过大的问题，过大的步长会导致体系不稳定。还好325K不算太高，还比较接近常用的300K，2fs的步长可以处理含氢的键的振动。可是假如我们要在400K的条件下运行动力学模拟的话，那模拟的步长就要缩减到1.5fs。
我们的起始结构是手工搭建的，不是通常常见的来自实验的晶体结构，所以我们的体系在模拟的开始阶段要面临不如晶体结构稳定的问题。为了让我们的体系能够在可控制的状况下来释放能量，在模拟起始的升温阶段我们选择步长为0.5fs，进入相对稳定的生成相之后，我们再选择常规的2fs步长0.5fs的步长确实有些矫枉过正，但是保证体系的安全毕竟还是最重要的。 我们进行升温模拟的方案如下：
第一阶段 - 10,000 步, 步长 0.5fs (共5ps), 起始温度 0.0K, 结束温度 50.0K, 温度耦合系数 1.0ps 第二阶段 - 10,000 步, 步长 0.5fs (共5ps), 结束温度 100.0K, 温度耦合系数 1.0ps 第三阶段 - 10,000 步, 步长 0.5fs (共5ps), 结束温度 150.0K, 温度耦合系数 1.0ps 第四阶段 - 10,000 步, 步长 0.5fs (共5ps), 结束温度 200.0K, 温度耦合系数 1.0ps 第五阶段 - 10,000 步, 步长 0.5fs (共5ps), 结束温度 250.0K, 温度耦合系数 1.0ps 第六阶段 - 10,000 步, 步长 0.5fs (共5ps), 结束温度 300.0K, 温度耦合系数 1.0ps 第七阶段 - 40,000 步, 步长 0.5fs (共20ps),结束温度 325.0K, 温度耦合系数 1.0ps
下面是第一阶段的输入文件：
heat1.in
Stage 1 heating of TC5b 0 to 50K &cntrl
imin=0, irest=0, ntx=1, nstlim=10000, dt=0.0005, ntc=2, ntf=2,
ntt=1, tautp=1.0,
tempi=0.0, temp0=50.0, ntpr=50, ntwx=50, ntb=0, igb=1,
cut=999.,rgbmax=999. /
其他六个阶段的输入文件与之非常接近，只需要改变一下相应的温度就可以了，可以从此处下载现成的输入文件： (heat2.in, heat3.in, heat4.in, heat5.in, heat6.in, heat7.in).
下面是一个运行升温模拟的PBS脚本，你也可以根据你的系统自己写一个脚本。
#PBS -l ncpus=16
#PBS -l walltime=500:00:00 #PBS -l cput=2000:00:00 #PBS -j oe
setenv AMBERHOME /usr/people/rcw/amber9 cd ~rcw/initial_heating
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat1.in -p TC5b.prmtop -c min1.rst -r heat1.rst -o heat1.out -x heat1.mdcrd
gzip -9 heat1.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat2.in -p TC5b.prmtop -c heat1.rst -r heat2.rst -o heat2.out -x heat2.mdcrd gzip -9 heat2.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat3.in -p TC5b.prmtop -c heat2.rst -r heat3.rst -o heat3.out -x heat3.mdcrd gzip -9 heat3.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat4.in -p TC5b.prmtop -c heat3.rst -r heat4.rst -o heat4.out -x heat4.mdcrd gzip -9 heat4.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat5.in -p TC5b.prmtop -c heat4.rst -r heat5.rst -o heat5.out -x heat5.mdcrd gzip -9 heat5.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat6.in -p TC5b.prmtop -c heat5.rst -r heat6.rst -o heat6.out -x heat6.mdcrd gzip -9 heat6.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat7.in -p TC5b.prmtop -c heat6.rst -r heat7.rst -o heat7.out -x heat7.mdcrd
gzip -9 heat7.mdcrd echo \
译者提供的bash脚本如下：
#!/bin/bash #heating
sander -O -i heat1.in -p TC5b.prmtop -c min1.rst -r heat1.rst -o heat1.out -x heat1.mdcrd gzip -9 heat1.mdcrd
sander -O -i heat2.in -p TC5b.prmtop -c heat1.rst -r heat2.rst -o heat2.out -x heat2.mdcrd gzip -9 heat2.mdcrd
sander -O -i heat3.in -p TC5b.prmtop -c heat2.rst -r heat3.rst -o heat3.out -x heat3.mdcrd gzip -9 heat3.mdcrd
sander -O -i heat4.in -p TC5b.prmtop -c heat3.rst -r heat4.rst -o heat4.out -x heat4.mdcrd gzip -9 heat4.mdcrd
sander -O -i heat5.in -p TC5b.prmtop -c heat4.rst -r heat5.rst -o heat5.out -x heat5.mdcrd gzip -9 heat5.mdcrd
sander -O -i heat6.in -p TC5b.prmtop -c heat5.rst -r heat6.rst -o heat6.out -x heat6.mdcrd gzip -9 heat6.mdcrd
sander -O -i heat7.in -p TC5b.prmtop -c heat6.rst -r heat7.rst -o heat7.out -x heat7.mdcrd gzip -9 heat7.mdcrd mkdir initial_heating
cp heat1.out initial_heating cp heat2.out initial_heating cp heat3.out initial_heating cp heat4.out initial_heating cp heat5.out initial_heating cp heat6.out initial_heating cp heat7.out initial_heating
cp heat1.mdcrd.gz initial_heating cp heat2.mdcrd.gz initial_heating cp heat3.mdcrd.gz initial_heating cp heat4.mdcrd.gz initial_heating cp heat5.mdcrd.gz initial_heating cp heat6.mdcrd.gz initial_heating cp heat7.mdcrd.gz initial_heating echo \
在16个1.3GHz CPU的SGI Altix上，7个步骤全部完成共耗时7分钟。下面提供了全部过程的输出文件，你可以分别下载，也可以下载最终的一个压缩文件。
Heating Stage Stage 1 Stage 2 Stage 3 Stage 4 Stage 5 Stage 6 Stage 7 Complete file set Output File heat1.out heat2.out heat3.out heat4.out heat5.out heat6.out heat7.out Restrt File heat1.rst heat2.rst heat3.rst heat4.rst heat5.rst heat6.rst heat7.rst Mdcrd File heat1.mdcrd.gz heat2.mdcrd.gz heat3.mdcrd.gz heat4.mdcrd.gz heat5.mdcrd.gz heat6.mdcrd.gz heat7.mdcrd.gz heating.tar.gz (5.2 Mb) 将轨迹文件用VMD打开就可以看到在升温过程中究竟发生了什么。你可以看到体系随着温度升高开始折叠，我们对这一阶段的轨迹并不关心，观看升温过程主要的目的在于确认整个升温过程一切OK，没有发生什么意外。 下图显示了升温过程结束后肽链的结构:
从这个结构我们看出，一些alpha螺旋已经形成了，但是这个结构距离最终的稳定折叠构像还有很长的路要走。
第五步：生产相动力学模拟
本教程分子模拟部分的最后一步是在325K条件下进行一个时间非常长的动力学模拟。在文献中他们做了50ns的动力学模拟，但是实际上我们看到的结果在模拟进行了20ns之后就已经呈现在人们面前了，之后继续进行的30ns模拟的意义仅仅在于说明之前的计算获得的就是一个稳定的结果。 尽管文献的作者发现在模拟的最初5~20ns中蛋白就已经充分折叠，我们在本教程中仍然进行50ns的动力学计算，以重复文献的报道。
\ns.\
我们将这个总时间长度为50ns的模拟分为10个阶段，每段5ns，这样做是为了一旦系统崩溃，我们不会损失已经进行的所有工作。而且这样分开处理还可以保证每个输出文件和轨迹文件的大小都适合处理。这10个阶段的模拟我们会使用相同的输入文件，文件如下所示:
equil.in Stage 2 equilibration 1 0-5ns &cntrl
imin=0, irest=1, ntx=5, nstlim=2500000, dt=0.002, ntc=2, ntf=2,
ntt=1, tautp=0.5,
tempi=325.0, temp0=325.0, ntpr=500, ntwx=500, ntb=0, igb=1,
cut=999.,rgbmax=999. /
每一个阶段的模拟都会进行250,000步(由nstlim的取值决定)，步长为2fs(由dt的取值决定) 即总共进行 5 ns的模拟。请注意我们在模拟全过程中使用了SHAKE(ntc=2, ntf=2)，我们使用了Berendsen 恒温器来保证模拟过程中体系温度恒定(ntt=1)，另外由于我们的体系已经完成升温，并且保持稳定，所以我们选择了耦合地更加紧密的恒温器(tautp=0.5)，这个恒温器的作用在于保持我们模拟的对象始终保持325K的温度。如同文献中所列，我们将模拟的温度设定在325K，每隔500步书写一次输出文件和轨迹文件. 过于频繁地写入这些文件会产生非常巨大的文件，这是不容易处理的。按照我们进行的设定，每5ns的模拟会产生一个 35 Mb的轨迹文件，对于我们的研究来说，500步记录一次已经足够了。
下面是进行生产相模拟的PBS 脚本，您可以根据您系统的状况修改脚本。
#PBS -l ncpus=16
#PBS -l walltime=500:00:00 #PBS -l cput=8000:00:00 #PBS -j oe
setenv AMBERHOME /usr/people/rcw/amber9 cd ~rcw/production
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c heat7.rst -r equil1.rst -o equil1.out -x equil1.mdcrd gzip -9 equil1.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil1.rst -r equil2.rst -o equil2.out -x equil2.mdcrd gzip -9 equil2.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil2.rst -r equil3.rst -o equil3.out -x equil3.mdcrd gzip -9 equil3.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil3.rst -r equil4.rst -o equil4.out -x equil4.mdcrd gzip -9 equil4.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil4.rst -r equil5.rst -o equil5.out -x equil5.mdcrd gzip -9 equil5.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil5.rst -r equil6.rst -o equil6.out -x equil6.mdcrd gzip -9 equil6.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil6.rst -r equil7.rst -o equil7.out -x equil7.mdcrd gzip -9 equil7.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil7.rst -r equil8.rst -o equil8.out -x equil8.mdcrd gzip -9 equil8.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil8.rst -r equil9.rst -o equil9.out -x equil9.mdcrd gzip -9 equil9.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil9.rst -r equil10.rst -o equil10.out -x equil10.mdcrd gzip -9 equil10.mdcrd echo \
在16个1.3GHz CPU的SGI Altix上全部的十段模拟共耗时27小时。下面是每一步模拟的结果输出，您可以下载这些结果，但是您需要注意每一个轨迹文件在压缩之后都有13Mb 左右
你可以将轨迹文件载入到VMD等看图软件中，提醒一下，你要首先解压缩这个轨迹文件，然后载入它，整个过程有可能需要花费一天时间，然后你就可以欣赏整个动力学过程了，如果你用飘带模型来显示的话会非常好看。下面是一个抓图： 第六步：分析模拟结果
下面我们要开始分析模拟的结果，我们要绘制模拟过程的温度、总能量、动能、势能变化曲线，这些信息可以从输出文件中获取，检查这些数据也可以让我们知道在动力学模拟过程中有没有发生什么异常状况。我们可以看到，温度的曲线非常平滑，并且维持在325K。动能和势能的曲线也一样，会在平衡位置比较平滑。这些曲线任何的突刺都暗示了我们的模拟过程发生了一些意外，我们就需要看看是不是用了不好的起始结构、太长的动力学步长或者不合适的参数等等。
为了从输出文件中提取我们需要的数据，我们可以使用下述 perl脚本： process_mdout.perl。用下面的命令我们可以把10个输出文件中的信息提取到一个文件中：
mkdir analysis cd analysis
../process_mdout.perl ../heat1.out ../heat2.out ../heat3.out ../heat4.out ../heat5.out ../heat6.out ../heat7.out ../equil1.out ../equil2.out ../equil3.out ../equil4.out ../equil5.out ../equil6.out ../equil7.out ../equil8.out ../equil9.out ../equil10.out
下面这个文件是运行process_mdout.perl提取出来的信息： analysis.tar.gz (2.3 Mb)
下面我们可以用一些图形化的数据处理软件如xmgr来分析动力学数据。xmgr我们在教程1中使用过，当时也是用来绘制温度和能量数据的曲线。 首先处理温度：
温度 升温相的温度 结果正如我们所愿，温度非常平缓地在325K附近波动，升温相也显示为平缓的上升，没有突然的跃升。 下面我们来看能量的状况，我们最感兴趣的是势能，我们可以找到势能最低的构像，把它作为标准的折叠模式，然后又可以得到一条RMSd的曲线 总能量(黑色)，势能 (红色)， 动能 (绿色) 看起来还不错，动能的曲线非常平滑，这从此前的温度曲线可以预测出来，因为温度与动能是直接相关的，所以温度稳定，动能也一定稳定。最初总能量和势能呈现减少的趋势，最终也趋于稳定，这表明我们的体系从最初的相对不稳定的状态折叠成一个稳定的状态，从图上我们能够比较明显地看到势能的衰减。我们也可以做一个每10ps取平均值的曲线，那样的话，势能衰减的趋势更容易看出来。要注意的是在 xmgr 中你可以通过Data->Transformations->Running averages来获得这样的一个曲线 势能 (红色曲线是10Ps平均的势能变化曲线) 上图是直接摘取自文献中的，我们教程中所使用的起始结构其能量要高于文献。这有可能是由于他们使用的是Amber6中的力场，而我们用的是Amber8，这导致两次模拟在广义波恩方程一项上产生差异。从 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) 2、 直接用mol2 转 也可以直接把sybyl 的mol2 做参数 * -antechamber –i file.mol2 –fi mol2 –o filr.prep –fo prepi –c resp / -c bcc * -parmchk –i file.prep –f prepi –o file.frcmod now you have 2 files: file.prep and file.frcmod for the drug . 注意修改file.prep 中文件名。要与pdb中的配体一致。（vi prep。假设是abc） 对应pdb与prep文件中的原子名称,，原子顺序要一致。 prep 在xleap中打开 using tleap for the protein or drug or/and complex -xleap –s –f leaprc.ff99 -mod = loadamberparams file.frcmod -loadamberprep file.prep edit abc pdb 可以用sybyl打开。用vi 修改pdb -source leaprc.gaff -RL = loadpdb file.pdb -Solvateoct RL TIP3PBOX 10 -Addions RL Cl- 0 (0:zero for neutral charge,Cl- can be replaced by Na+) -Saveamberparm RL fileWT.top fileWT.crd (save topology and coord. With water) amber动力学，博大精深！做参数还需好好学习，记录一下过程，以免忘记。 一般来说，小分子的参数是最麻烦的，首先具备的是mol2文件，做参数主要是电荷问题，一般有bcc和resp电荷，个人觉得选用resp电荷可能好一点，曾经用bcc发现优化做不下去。 1、用resp电荷前，需用gaussian计算电荷分布，其中gaussian计算的关键词为： %mem=400Mb %nproc=1 #P HF 6-31G* pop=mk iop(6/33=2) iop(6/42=6) 2、生成aaa.out文件后，利用antechamber制作prep文件 antechamber -i aaa.out -fi gout -o aaa.prep -fo prepi -c resp #输入文件名、文件类型、输出文件名、输出文件类型，电荷类型 3、产生prep文件后，就是检查参数， parmchk -i aaa.prep -fi prepi -o aaa.frcmod 在aaa.frcmod文件中是一些非标准的键长，键角，二面角等信息。 4、修改prep文件中，小分子的名字，默认的都是MOL，改成pdb中复合物中的小分子一致的名字（设为ABC )。 5、利用xleap修改原子名称。 由于mol2文件转成的prep文件的原子名称往往与pdb中小分子的名称不一致，这是可以用xleap中的edit ABC。（假设pdb中小分子的名字叫ABC). 在图形中选择显示原子名称，这样这就是amber将会读的原子名称，是prep文件中的。可以用sybyl打开pdb，现在原子id，对应位置，利用vi编辑器，打开pdb，修改pdb中的原子名称与xleap中的一致。 6、修改工作一般很细小，一定要小心，完成后保存退出。为了以防万一，备份一下pdb。然后vi打开pdb ，将原子键连的信息删掉。小分子与蛋白间有必要可以加上TER，断开各个小分子，这是pdb格式中的一个通用分段方式，具体可以查看pdb的文件格式信息。 7、完成后就可以sus=loadpdb abc.pdb导入pdb，这样一般不会有新增加原子了。 在做动力学的时候有时候出现以下错误信息： Coordinate resetting (SHAKE) cannot be accomplished, deviation is too large 这时候可以用ptraj中的checkoverlap检查原子间是否有很近的原子 ptraj << EOF YourParm.top trajin YourCRD.crd checkoverlap go EOF