Amber 教程

$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