Amber软件中动力学模拟的步骤.docx
《Amber软件中动力学模拟的步骤.docx》由会员分享，可在线阅读，更多相关《Amber软件中动力学模拟的步骤.docx（13页珍藏版）》请在冰点文库上搜索。
Amber软件中动力学模拟的步骤
MolecularDynamicssimulation
——从能量最小化到实际模拟
1基本流程图
1）概述
前面我们已经得到了Amber用来动力学模拟的prmtop和inpcrd文件，它们分别是参数文件和坐标文件。
我们先从一条命令说起来解释Amber是如何做动力学模拟的：
sander–O–imdin–omdout–pprmtop–cinpcrd–rrst–xmdcrd
动力学过程是一个连续地解牛顿运动方程的过程：
上一个牛顿方程结束时，蛋白质中各原子的位置和速度保留给下一个牛顿方程，惟一改变的是原子的加速度，它会根据各种势能函数重新计算（势能随原子坐标改变：
E=f（r,…））。
只不过每个牛顿方程的时间很短，短到fs（10-15s）级，Amber软件提供的sander主程序可以用来自动地做这样的数值计算。
它需要参数文件（prmtop）、坐标文件（inpcrd）、sander程序配置文件（mdin）来启动运行，我们已经有了前两种文件，本节内容最主要的就是讲解如何配置我们需要的动力学模拟。
sander程序运行过程中会输出临时文件（rst）保存坐标和速度，还有轨迹文件（mdcrd）。
2）动力学过程
从基本流程图可以知道，一般的动力学过程也就可以分为三步：
能量最小化（minimization）、体系平衡（equilibrium）、实际动力学模拟。
由于我们进行的初始结构来自晶体结构或同源建模，所以在分子内部存在着一定的结构张力，能量最小化就是真正的动力学之前释放这些张力，如果没有这个步骤，在动力学模拟开始之后，整个体系可能会因此变形、散架。
另外，由于动力学模拟的是真实的生物体环境，因此必须使研究对象升温升压到临界值，体系达到平衡，才能做实际的动力学模拟。
2各流程输入文件
要通过Amber软件做动力学模拟，需要明白如何去配置上述过程中的每一步。
一般来说就是指定一些键/值对。
1）Minimization
第一行是标题，&cntrl是起始符，”/”是结束符，中间的键/值对就是参数配置。
上述的参数配置可以归纳如下：
imin=1
Chooseaminimizationrun，指定做能量最小化的动力学。
ntx=1
ReadcoordinatesbutnotvelocitiesfromASCIIformattedinpcrd coordinatefile
ntx指定如何获得坐标信息，这里直接从inpcrd中读取坐标信息。
maxcyc=2000
Maximumminimizationcycles
能量最小化的算法涉及循环迭代，这里指定迭代次数。
ncyc=1000
Thesteepestdescentalgorithmforthefirst0-ncyc cycles,thenswitchestheconjugategradientalgorithmfor ncyc-maxcyccycles
如上所说，循环迭代的算法不同，此处指定到哪一步第一种迭代算法结束。
ntpr=100
PrinttotheAmber mdout outputfileevery ntpr cycles
这个表明采集计算信息的频率，输出到mdout文件中。
cut=8.0
NonbondedcutoffdistanceinAngstroms
由于计算能量时需要有一个截断距离，这个参数指定截断距离。
其余两个参数都是能量最小化时不可用的参数，下面会具体说明。
2）Heating（Equilibriumstep1）
上述的参数配置可以归纳如下：
imin=0
Chooseamoleculardynamics（MD）run[nominimization]
开始做动力学，而不是能量最小化。
nstlim=10000
NumberofMDstepsinrun（nstlim*dt=runlengthinps）
动力学过程的步长，前面说过动力学的原理，它是连续地解牛顿运动方程。
动力学的时长会等于步长（nstlim）乘以每一步的时间间隔（如下，dt）。
dt=0.002
Timestepinpicoseconds（ps）.ThetimelengthofeachMDstep
每一步的时间间隔，单位是ps。
ntf=2
SettingtonotcalculateforceforSHAKEconstrainedbonds
ntc=2
EnableSHAKEtoconstrainallbondsinvolvinghydrogen
以上两个参数针对氢原子做shake限制，主要由于氢原子的振动频率过高，氢原子的运动还存在量子效应。
tempi=0.0
InitialthermostattemperatureinK（see NMROPT section）
temp0=300.0
FinalthermostattemperatureinK（see NMROPT section）
以上两个参数指定初始温度和升温后的温度。
ntwx=1000
WriteAmbertrajectoryfile mdcrd every ntwx steps
将坐标信息写入轨迹文件的频率。
ntb=1
Periodicboundariesforconstantvolume
采用周期性边界的恒容条件，表明这是一个NVT系综，周期性边界在附录中说明。
ntp=0
Nopressurecontrol
暂时不考虑控制压力。
ntt=3
TemperaturecontrolwithLangevinthermostat
gamma_ln=2.0
Langevinthermostatcollisionfrequency
以上两个参数指定如何控制温度。
nmropt=1
NMRrestraintsandweightchangesread（see NMROPTsection）
这个参数原本用于指定nmr限制，现在被机智地借来用于控制升温过程。
升温过程的控制：
这三行配置表示，在10000步的升温中，在最初的9000步中，温度将从0K升到300K，最后9001步到10000步，温度会保持在300K。
3）Equilibriumstep2、Production
注意：
这里配置的动力学时间只有60ps，主要是为了演示；一般来说视体系的大小和研究的需要，模拟的时间一般会更长。
第二步平衡（体系升压）与实际动力学的参数配置是一样的，因为实际动力学便是在恒温恒压下进行的。
上述的参数配置可以归纳如下：
ntx=5
Readcoordinatesandvelocitiesfromunformatted inpcrdcoordinatefile
ntx等于5表示将读取另一种格式的inpcrd文件，它往往配合irest=1（如下）使用。
irest=1
RestartpreviousMDrun[Thismeansvelocitiesareexpectedintheinpcrdfileandwillbeusedtoprovideinitialatomvelocities]
irest等于1表明将重启上一次的动力学，也就意味着希望能从输入文件中读取到上一次动力学结束时对应的（除了坐标信息外的）原子速度信息（储存在rst文件中）。
temp0=300.0
Thermostattemperature.Runat300K
动力学过程中需要控制保持的温度，一般在300K左右
ntb=2
Useperiodicboundaryconditionswithconstantpressure
采用周期性边界的恒压条件。
一般动力学即采用此NPT系综。
ntp=1
UsetheBerendsenbarostatforconstantpressuresimulation
设定控压的算法。
3.进入Amber实际操作
1）min1：
限制蛋白，溶剂部分能量最小化
终端运行命令：
$AMBERHOME/bin/sander-O-imin1.in-omin1.out–pNAME.prmtop-cNAME.inpcrd–rmin1.rst-refNAME.inpcrd
$AMBERHOME指定Amber软件安装包的目录信息。
NAME.prmtop、NAME.inpcrd即表示prmtop文件盒inpcrd文件。
-ref表示按照给定的坐标信息来限制蛋白，这里按照蛋白的inpcrd内的坐标信息。
2）min2：
放松蛋白，整个体系能量最小化
终端运行命令：
$AMBERHOME/bin/sander-O-imin2.in-omin2.out–pNAME.prmtop-cmin1.rst-rmin2.rst
输入的坐标信息来源于上一步输出的rst文件，其储存了坐标信息。
下面的步骤是同理的，上一步输出的rst文件将作为下一步输入的坐标文件。
3）Equilibrationstep1：
（heating：
系统在约束蛋白下升温）
给系统加热，从0K到300K，运行20ps动力学（nstlim的值可以适当提高，视体系而定）：
终端运行命令：
$AMBERHOME/bin/sander-O-imd1.in-omd1.out–pNAME.prmtop-cmin2.rst-rmd1.rst-xmd1.mdcrd-refmin2.rst
这一步开始会收集轨迹文件（-x指定的输出），用于后续可能的分析。
4）Equilibrationstep2：
（体系升压）
为了平衡体系，将运行100ps的恒温恒压的动力学计算：
终端运行命令：
$AMBERHOME/bin/sander-O-imd2.in-omd2.out-pNAME.prmtop–cmd1.rst-rmd2.rst–xmd2.mdcrd
5）ProductionRun：
实际模拟
获取轨迹，将运行10ns的分子动力学计算。
一般来说，实际模拟中的配置文件与系统平衡时的配置文件是一致的：
这里只是”nstlim”这个配置参数的值改变了。
终端运行命令：
$AMBERHOME/bin/sander–O–iprod.in–oprod.out–pNAME.prmtop–cmd2.rst–rprod.rst–xprod.mdcrd
最终我们将对prod.out、prod.mdcrd进行分析、整合，这个将再下一节中讲述。
4.不同的sander运行命令：
集群限制
本阶段可以实现多种运行方式，如串行、并行、多集群并行和GPU计算等等。
整理执行命令如下。
1）串行：
$AMBERHOME/bin/sander–O–imdin–omdout–pprmtop–c–rrst–xmdcrd
$AMBERHOME/bin/pmemd–O–imdin–omdout–pprmtop–c–rrst–xmdcrd
Amber中有两个模块程序支持动力学模拟：
sander和pmemd。
2）并行：
mpirun–npN$AMBERHOME/bin/sander.MPI–O–imdin–omdout–pprmtop–c–rrst–xmdcrd
mpirun–npN$AMBERHOME/bin/pmemd.MPI–O–imdin–omdout–pprmtop–c–rrst–xmdcrd
N值为需要加载的服务器的核数，注意，N值最大值为运行服务器的CPU个数，例如8、12等。
3）多集群并行：
mpirun–hostfilehosts–npN$AMBERHOME/bin/sander.MPI–O–imdin–omdout–pprmtop–c–rrst–xmdcrd
mpirun–hostfilehosts–npN$AMBERHOME/bin/pmemd.MPI–O–imdin–omdout–pprmtop–c–rrst–xmdcrd
hosts用于指定并行计算的节点，这种方法多用于跨节点的并行计算中，一般会提供专门用于提交任务的脚本。
4）单块GPU计算
$AMBERHOME/bin/pmemd.cuda–O–imdin–omdout–pprmtop–c–rrst–xmdcrd
只有pmemd模块支持GPU。
5.附录：
1）周期性边界条件简介
一个二维盒子排列，分子1将从中心盒子（有阴影的盒子）运动到盒子B中，这时为了保持中心盒子的粒子守恒，盒子镜像F中有一个相应的粒子会进入到中心盒子中来，就好像分子1从上缘出去，而从下缘进入。
周期性边界条件使得粒子包裹在无穷多的溶液分子中，这样在使用相对少量的粒子的情况下模拟了真实的生物体系，使动力学模拟的成本大大减少。
假设一个有粒子在内部运动的立方体盒子从各个方向被复制，形成了周期性排列的盒子，就像上图的二维盒子展示的一样。
在二维的情况下，每个盒子周围都有8个相邻的盒子，而在三维的情况下，每个盒子周围则会有26个相邻的盒子。
对于镜像盒子，相当于中心盒子朝各个方向平移了，因此镜像盒子中粒子的坐标也就等于中心盒子的坐标加上或减去平移向量：
r + ix + jy + kz （i,j,k ->-inf,+inf）。
如果一个粒子在模拟过程中从某一面离开了盒子，那么相应地就会有一个镜像粒子从另一面进入盒子，始终保持盒子中粒子数目不变。
立方体是最简单的可用于周期性排列的几何形状，视图查看和程序化都很方便；然而，有些情况下另一种盒子形状更适合模拟体系，例如被水分子包围的单个大分子或复合物，往往需要研究的对象也就是被包裹在中心的大分子、复合物，因此周围溶剂（水分子）的计算量越少越好，也就是说尽可能减少水分子的数目。
在这种情况下，立方体盒子事实上是最不节省的盒子。
原则上，任何形状的盒子都可以用于周期性模拟，只要通过平移操作，能够填满整个三维空间（也就是中心盒子与镜像盒子间没有空隙）。
因此大概有5类盒子满足条件：
立方体（同样的，长方体也可以）、六方柱、截八面体、菱形十二面体。
另外，在选取合适的盒子形状后，模拟中盒子的大小也是非常重要的，既要盒子足够大，保证模拟的精确，又要尽量减小盒子的大小降低计算量。
体系的大小一般要大于体系中相互作用能的计算范围的两倍。
6.Reference
http:
//ambermd.org/tutorials/basic/tutorial1/section5.htm
http:
//www.charmm-gui.org/?
doc=input/mdsetup
升级会员 