我们的起始结构是手工搭建的，不是通常常见的来自实验的晶体结构，所以我们的体系在模拟的开始阶段要面临不如晶体结构稳定的问题。为了让我们的体系能够在可控制的状况下来释放能量，在模拟起始的升温阶段我们选择步长为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，这导致两次模拟在广义波恩方程一项上产生差