AMBER分子动力学简例（三）
分子动力学（2）
水环境中的分子动力学模拟
溶剂环境中的分子动力学模拟分为以下四步进行：
1、溶剂环境能量最优化。这一步保持溶质（蛋白）不变，去除溶剂中能量不正常的范德华相 互作用。
2、整系统能量最优化。去除整个系统中能量不正常的相互作用。
3、有限制的分子动力学。保持蛋白质不动，溶解溶剂的不同层，同时逐渐将系统温度从0K提升到300K。
4、整系统分子动力学模拟。在一个大气压，300K的环境下整个系统分子动力学模拟。可以得到成果的分子动力学模拟。
#############################
1、溶剂环境能量最优化。
该步骤的配置文件min1.in如下：
--------------------------------------------------------------------------------
oxytocin: initial minimisation solvent + ions
&cntrl
imin = 1,
maxcyc = 1000,
ncyc = 500,
ntb = 1,
ntr = 1,
cut = 10
/
Hold the protein fixed
500.0
RES 1 9
END
END
--------------------------------------------------------------------------------
该过程保持肽链不动，其中500.0单位是kcal/mol，表示作用在肽链上使其不动的力。“RES 1 9”表示肽链残基数目，因为我们学习使用的oxytocin有9个残基。
模拟命令如下：
sander –O –i min1.in –o min1.out –p oxy.top –c oxy.crd –r oxy_min1.rst –ref oxy.crd &
2、整系统能量最优化。
配置文件min2.in如下：
--------------------------------------------------------------------------------
oxytocin: initial minimisation whole system
&cntrl
imin = 1,
maxcyc = 2500,
ncyc = 1000,
ntb = 1,
ntr = 0,
cut = 10
/
--------------------------------------------------------------------------------
命令如下：
sander –O –i min2.in –o min2.out –p oxy.top –c oxy_min1.rst –r oxy_min2.rst &
3、有限制的分子动力学。
第一步分子动力学保持蛋白分子位置不变，但是不是完全固定每个原子，同时缓解蛋白分子周围的水分子，是溶剂环境能量优化。在这个步骤中，我们将主要目的是对特定的原子使用作用力使其能量优化。我们要优化溶剂环境，至少需要10ps，我们将使用20ps用来优化我们上两步制作的分子系统的周期性边界的溶剂环境。
命令配置文件md1.in如下：
--------------------------------------------------------------------------------
oxytocin: 20ps MD with res on protein
&cntrl
imin = 0,
irest = 0,
ntx = 1,
ntb = 1,
cut = 10,
ntr = 1,
ntc = 2,
ntf = 2,
tempi = 0.0,
temp0 = 300.0,
ntt = 3,
gamma_ln = 1.0,
nstlim = 10000, dt = 0.002,
ntpr = 100, ntwx = 500, ntwr = 1000
/
Keep protein fixed with weak restraints
10.0
RES 1 9
END
END
--------------------------------------------------------------------------------
上述参数解释如下：
ntb = 1：表示分子动力学过程保持体积固定。
imin = 0：表示模拟过程为分子动力学，不是能量最优化。
nstlim = #：＃表示计算的步数。
dt = 0.002：表示步长，单位为ps，0.002表示2fs。
temp0 = 300：表示最后系统到达并保持的温度，单位为K。
tempi = 100：系统开始时的温度。
gamma_ln = 1：表示当ntt＝3时的碰撞频率，单位为ps-1（请参考AMBER手册）
ntt = 3：温度转变控制，3表示使用兰格氏动力学。
tautp = 0.1：热浴时间常数，缺省为1.0。小的时间常数可以得到较好的耦联。
vlimit = 20.0：保持分子动力学稳定性速度极限。20.0为缺省值，当动力学模拟中原子速度大于极限值时，程序将其速度降低到极限值以下。
comp = 44.6： 溶剂可压缩单位。
ntc = 2：Shake算法使用标志位。1表示不实用使用，2表示氢键将被计算，3表示所有键都将被计算在内。
tol = #.#####：坐标位置重新设置的几何位置相对容忍度。
我们将使用一个较小的作用力，10kcal/mol。在分子动力学中，当ntr＝1时，作用力只需要5-10kcal/mol（我们需要引用一个坐标文件做分子动力学过程的比较，我们需要使用"-ref"参数）。太大的作用力同时使用Shake算法和2fs步长将使整个系统变得不稳定，因为大的作用力使系统中的原子产生大频率的振动，模拟过程并步需要。
运行命令如下：
sander –O –i md1.in –o md1.out –p oxy.top –c oxy_min2.rst –r oxy_md1.rst –x oxy_md1.mdcrd –ref oxy_min2.rst –inf md1.info&
进行MD的运行时间一般较长，可以使用程序的并行版本提交集群计算。在主频位2.0GHz的P4单机上，大概需要一个钟头。可以随时查看md1.in文件的程序输出。
4、整系统分子动力学模拟。
这一步中，我们将进行整个系统的分子动力学模拟，而不对某些特定原子位置进行限制。因为知识一个小例子，我们将只进行250ps的MD计算。配置文件md2.in如下：
--------------------------------------------------------------------------------
oxytocin: 250ps MD
&cntrl
imin = 0, irest = 1, ntx = 7,
ntb = 2, pres0 = 1.0, ntp = 1,
taup = 2.0,
cut = 10, ntr = 0,
ntc = 2, ntf = 2,
tempi = 300.0, temp0 = 300.0,
ntt = 3, gamma_ln = 1.0,
nstlim = 125000, dt = 0.002,
ntpr = 100, ntwx = 500, ntwr = 1000
/
--------------------------------------------------------------------------------
上面一些参数解释如下:
ntb=2：表示分子动力学过程的压力常数。
ntp＝1：表示系统动力学过程各向同性。
taup = 2.0：压力缓解时间，单位为ps。
pres＝1：引用1个单位的压强。
使用以下命令进行MD：
sander –O –i md2.in –o md2.out –p oxy.top –c oxy_md1.rst – r oxy_md2.rst –x oxy_md2.mdcrd –ref oxy_md1.rst –inf md2.info &
模拟时间较长，在P4单机2.0GHz的上需要7.5个小时。
到此，模拟全部完成，接下来要对得到的数据进行分析。主要数据文件.out文件，包含系统能量、温度，压力等等；.mdcrd文件，是分子动力学轨迹文件，可以求系统蛋白的RMSD，回转半径等等。数据分析根据不同的研究目的不同而不同，我们将在后文中进行一些简单的分析。
