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）。MinimizationEquilibriumProductionRun2）动力学过程从基本流程图可以知道，一般的动力学过程也就可以分为三步：能量最小化（minimization）、体系平衡（equilibrium）、实际动力学模拟。由于我们进行的初始结构来自晶体结构或同源建模，所以在分子内部存在着一定的结构张力，能量最小化就是真正的动力学之前释放这些张力，如果没有这个步骤，在动力学模拟开始之后，整个体系可能会因此变形、散架。另外，由于动力学模拟的是真实的生物体环境，因此必须使研究对象升温升压到临界值，体系达到平衡，才能做实际的动力学模拟。2各流程输入文件要通过Amber软件做动力学模拟，需要明白如何去配置上述过程中的每一步。一般来说就是指定一些键/值对。1）Minimization第一行是标题，&cntrl是起始符，”/”是结束符，中间的键/值对就是参数配置。上述的参数配置可以归纳如下：imin=1Chooseaminimizationrun，指定做能量最小化的动力学。ntx=1ReadcoordinatesbutnotvelocitiesfromASCIIformattedinpcrdcoordinatefilentx指定如何获得坐标信息，这里直接从inpcrd中读取坐标信息。maxcyc=2000Maximumminimizationcycles能量最小化的算法涉及循环迭代，这里指定迭代次数。ncyc=1000Thesteepestdescentalgorithmforthefirst0-ncyccycles,thenswitchestheconjugategradientalgorithmforncyc-maxcyccycles如上所说，循环迭代的算法不同，此处指定到哪一步第一种迭代算法结束。ntpr=100PrinttotheAmbermdoutoutputfileeveryntprcycles这个表明采集计算信息的频率，输出到mdout文件中。cut=8.0NonbondedcutoffdistanceinAngstroms由于计算能量时需要有一个截断距离，这个参数指定截断距离。其余两个参数都是能量最小化时不可用的参数，下面会具体说明。2）Heating（Equilibriumstep1）上述的参数配置可以归纳如下：imin=0Chooseamoleculardynamics(MD)run[nominimization]开始做动力学，而不是能量最小化。nstlim=10000NumberofMDstepsinrun(nstlim*dt=runlengthinps)动力学过程的步长，前面说过动力学的原理，它是连续地解牛顿运动方程。动力学的时长会等于步长（nstlim）乘以每一步的时间间隔（如下，dt）。dt=0.002Timestepinpicoseconds(ps).ThetimelengthofeachMDstep每一步的时间间隔，单位是ps。ntf=2SettingtonotcalculateforceforSHAKEconstrainedbondsntc=2EnableSHAKEtoconstrainallbondsinvolvinghydrogen以上两个参数针对氢原子做shake限制，主要由于氢原子的振动频率过高，氢原子的运动还存在量子效应。tempi=0.0InitialthermostattemperatureinK(seeNMROPTsection)temp0=300.0FinalthermostattemperatureinK(seeNMROPTsection)以上两个参数指定初始温度和升温后的温度。ntwx=1000WriteAmbertrajectoryfilemdcrdeveryntwxsteps将坐标信息写入轨迹文件的频率。ntb=1Periodicboundariesforconstantvolume采用周期性边界的恒容条件，表明这是一个NVT系综，周期性边界在附录中说明。ntp=0Nopressurecontrol暂时不考虑控制压力。ntt=3TemperaturecontrolwithLangevinthermostatgamma_ln=2.0Langevinthermostatcollisionfrequency以上两个参数指定如何控制温度。nmropt=1NMRrestraintsandweightchangesread(seeNMROPTsection)这个参数原本用于指定nmr限制，现在被机智地借来用于控制升温过程。升温过程的控制：这三行配置表示，在10000步的升温中，在最初的9000步中，温度将从0K升到300K，最后9001步到10000步，温度会保持在300K。3）Equilibriumstep2、Production注意：这里配置的动力学时间只有60ps，主要是为了演示；一般来说视体系的大小和研究的需要，模拟的时间一般会更长。第二步平衡（体系升压）与实际动力学的参数配置是一样的，因为实际动力学便是在恒温恒压下进行的。上述的参数配置可以归纳如下：ntx=5Readcoordinatesandvelocitiesfromunformattedinpcrdcoordinatefilentx等于5表示将读取另一种格式的inpcrd文件，它往往配合irest=1（如下）使用。irest=1RestartpreviousMDrun[Thismeansvelocitiesareexpectedintheinpcrdfileandwillbeusedtoprovideinitialatomvelocities]irest等于1表明将重启上一次的动力学，也就意味着希望能从输入文件中读取到上一次动力学结束时对应的（除了坐标信息外的）原子速度信息（储存在rst文件中）。temp0=300.0Thermostattemperature.Runat300K动力学过程中需要控制保持的温度，一般在300K左右ntb=2Useperiodicboundaryconditionswithconstantpressure采用周期性边界的恒压条件。一般动力学即采用此NPT系综。ntp=1UsetheBerendsenbarostatforconstantpressuresimulation设定控压的算法。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.mdcrd5）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–xmdcrdAmber中有两个模块程序支持动力学模拟：sander和pmemd。2）并行：mpirun–npN$AMBERHOME/bin/sander.MPI–O–imdin–omdout–pprmtop–c–rrst–xmdcrdmpirun–npN$AMBERHOME/bin/pmemd.MPI–O–imdin–omdout–pprmtop–c–rrst–xmdcrdN值为需要加载的服务器的核数，注意，N值最大值为运行服务器的CPU个数，例如8、12等。3）多集群并行：mpirun–hostfilehosts–npN$AMBERHOME/bin/sander.MPI–O–imdin–omdout–pprmtop–c–rrst–xmdcrdmpirun–hostfilehosts–npN$AMBERHOME/bin/pmemd.MPI–O–imdin–omdout–pprmtop–c–rrst–xmdcrdhosts用于指定并行计算的节点，这种方法多用于跨节点的并行计算中，一般会提供专门用于提交任务的脚本。4）单块GPU计算$AMBERHOME/bin/pmemd.cuda–O–imdin–omdout–pprmtop–c–rrst–xmdcrd只有pmemd模块支持GPU。5.附录：1）周期性边界条件简介一个二维盒子排列，分子1将从中心盒子（有阴影的盒子）运动到盒子B中，这时为了保持中心盒子的粒子守恒，盒子镜像F中有一个相应的粒子会进入到中心盒子中来，就好像分子1从上缘出去，而从下缘