AMBER分子动力学简例

AMBER分子动力学简例AMBER分子动力学简例（一）概述以下是使用AMBER包的简单教程，希望对开始学习分子动力学的同学有用处。申明一下，以下教程原版来自网上，是最最基本的教程，同时也非常实用，有非常好的借鉴意义。AMBER分子动力学程序包是加州圣弗兰西斯科大学（UniversityofCaliforniaSanFrancisco,UCSF)的PeterAKollman和其同事编的，程序很全，现在已经发展到版本9.0。AMBER功能涵盖种类非常多的生物分子，包括蛋白、核算以及药物小分子。软件详细情况请浏览软件包中四个主要的大程序：Leap：用于准备分子系统坐标和参数文件，有两个程序：xleap：X－windows版本的leap，带GUI图形界面。tleap：文本界面的Leap。Antechamber：用于生成少见小分子力学参数文件的。有的时候一些小分子Leap程序不认识，需要加载其力学参数，这些力学参数文件就要antechamber生成。Sander：MD数据产生程序，即MD模拟程序，被称做AMBER的大脑程序。Ptraj：MD模拟轨迹分析程序。学习项目本教程研究的题目是脑下垂体荷尔蒙之一的oxytocin，需要X光衍射晶体结构文件1NPO.PDB。该文件包含了该荷尔蒙和其运载蛋白的复合物，可以从蛋白数据库下载。PDB文件是不包含氢原子的，Leap程序会自动的加上PDB文件缺少的东西。当第一次使用PDB文件的时候，要十分留意文件包含的信息，所以PDB文件缺少的残疾、侧链或者添加的变异都在这个地方记录。可以用文本阅读程序阅读PDB文件头部的信息，即以REMARKS开头的信息文本行。PDB一个重要的信息是SSBOND记录，该记录说明结构中二硫键的位置，这样的信息在使用Leap程序建立分子结构的时候需要。在本教程中，我们将比较oxytocin在真空和溶液中分子动力学的差异，如果没有二硫键，将影响整个结果。整个过程在Linux系统下完成，如对Linux不熟悉，请翻阅有关Linux书籍。建立项目目录，并进入该目录：mkdirmyprojcdmyproj下载1NPO.PDB文件到该目录下，使用文本阅读器阅读该PDB文件。从在文件的开头部分可以得到该文件是一个二聚物的PDB文件，删除文件中A、C、和D链对应的文本行，剩下的就是oxytocin的晶体机构了，保存为oxyt.pdb文件。使用SwissPDBviewer分析以下oxyt.pdb文件，可以得知该pdb文件缺少了一个侧链。如果使用DeepView，它会自动给文件加上侧链。重新保存pdb文件为oxyt.pdb文件。xleap和tleap程序功能是一样的，都是准备分子结构的坐标文件和拓扑文件。xleap启动一个X界面，比较慢。tleap纯文本，要快一些，我们选择tleap：tleap-s-fleaprc.ff03其中leaprc.ff03是AMBER的2003力场文件。力场是一个很重要的文件，定义了分子、原子和残疾等的信息，请查阅有关资料。oxy=loadpdboxyt.pdb该命令读入oxyt.pdb文件到oxy变量中。bondoxy.1.SGoxy.6.SG连接oxy变量的第一和第六个残疾的SG原子，即是连接二硫键。可以使用check命令检查oxy是否完好，没有特殊情况的话，最后应该是“OK”，表示分子系统状态是好的，可以保存。checkoxy接下来就是保存分子系统了：saveamberparmoxyoxy_vac.topoxy_vac.crd其中oxy_vac.top和oxy_vac.crd分别是分子系统真空状态下的拓扑和坐标文件。接下来要给分子系统添加水环境，solvateoctoxyTIP3PBOX9.0该命令让Leap程序使用TIP3PBOX水模型,水环境距离分子为9纳米。也可以用solvatebox命令，但是solvatebox添加的水环境是一个立方体，不是八面体。一般要求水表面到蛋白距离要达到8.5纳米，避免MD过程中蛋白冲出水环境，但是水环境越大计算时间将越长。如果MD过程内含PME，那么水箱长度要大于2倍截矩（cutoff），截矩在MD配置文件中定义。加溶剂之后，要计算系统是否为电中性，使用命令：chargeoxy如果现实不是电中性，要使用addions添加反性电荷，如Cl－或者Na＋等。最后保存溶剂环境的系统拓扑结构和坐标文件，并退出Leap程序：saveamberparmoxyoxy.topoxy.crdquit也可以把命令集中到一个，让tleap读取。如将上面的命令集中到以下文件中oxy.leaprc--------------------------------------------------------------------------------sourceleaprc.ff03oxy=loadpdboxyt.pdbbondoxy.1.SGoxy.6.SGcheckoxysaveamberparmoxyoxy_vac.topoxy_vac.crdsolvateoctoxysaveamberparmoxyoxy.topoxy.crdquit将以上两横线之间的命令保存为oxy.leaprc，然后使用以下命令一步搞定：tleap-s-foxy.leaprc到此，分子系统准备已经完成，目录中会多了一个leap.log文件，这是Leap程序的记录文件，它包含了leap程序所做的一切，包括自动的和手动命令。现在已经有了可以进行分子动力学的基本文件了，再编写sander的控制文件，就可以进行模拟了。这些将在后文中介绍。AMBER分子动力学简例（二）分子动力学（1）真空模式真空模式分子动力学模拟将使用NVT系宗分两步进行，即系统能量最优化和分子动力学过程。1、系统能量最优化。我们将使用淬火能量最优化解除系统内原子之间的不正常相互作用，这些原子之间的高能量相互作用如果不消除，可能影响后续的分子动力学过程。因为动力学过程是能量梯度变化的，太高的能量壁垒可能让MD局限在某一个能量局部最小化位置中。Sander程序是分子动力学模拟程序，它的主要功能是能量最优化，动力学模拟和NMR优化计算。我们必须给Sander一个运行的配置文件，使其按照我们的要求进行计算。能量最优化的配置文件如下：min_vac.inoxytocin:initialminimizationpriortoMD&cntrlimin=1,maxcyc=500,ncyc=250,ntb=0,igb=0,cut=12/--------------------------------------------------------------------------------Sander程序的基本命令参数如下：sander–O–iin–oout–pprmtop–cinpcrd–rrestrt[-refrefc–xmdcrd–vmdvel–emden–infmdinfo]所以我们的命令可以如下：sander–O–imin_vac.in–omin_vac.out–poxy_vac.top–coxy_vac.crd–roxy_vacmin.rst&可以使用more命令或者tail命令查看min_vac.out的输出内容，那是能量最优化的记录。2、分子动力学模拟。Sander程序的配置文件为：md_vac.inoxytocinMDin-vacuo,12angstromcutoff,250ps&cntrlimin=0,ntb=0,igb=0,ntpr=100,ntwx=500,ntt=3,gamma_ln=1.0,tempi=300.0,temp0=300.0,nstlim=125000,dt=0.002,cut=12.0/--------------------------------------------------------------------------------命令为：sander–O–imd_vac.in–omd_vac.out–poxy_vac.top–coxy_vacmin.rst–roxy_vacmd.rst–xoxy_vacmd.mdcrd–refoxy_vacmin.rst–infmdvac.info&MD的计算过程一般比较久，真空相对与溶剂中要快以下，250ps的模拟大概在一个主频为2。0GHz的Linux单机上运行5分钟。以上配置文件中用到很多参数，这些参数这是sander程序参数的一小部分，以下将相应解释。如果要了解sander的其他参数，请阅读AMBER用户指南。&ctrl和/：sander的参数一般要求出现在这两个标识符号之间，参数以及这两个标识符被称做控制模块。cutoff：以纳米为单位的截矩。即超出截矩范围的非键连接相互作用将不计。ntr：原子位置能量抑制位，1表示抑制，0表示不抑制。imin：能量最优化标志位。1表示sander将进行能量最优化，0表示让sander进行分子动力学模拟。macyc：能量最优化次数。ncyc：便是经过多少次能量优化以后，能量优化从淬火过程变为梯度变化过程。ntmin：能量优化方法标志位。0表示前10个能量最优化为淬火过程，然后进行梯度能量优化；1表示ncyc次淬火过程，然后进行梯度能量优化，为默认值；2表示只进行淬火过程。dx0：表示启动模拟步长。dxm：最大优化步数。drms：梯度能量优化标准，默认值为1.0E-4kcal/mol.A。更多参数将在后文中解释。AMBER分子动力学简例（三）分子动力学（2）水环境中的分子动力学模拟溶剂环境中的分子动力学模拟分为以下四步进行：1、溶剂环境能量最优化。这一步保持溶质（蛋白）不变，去除溶剂中能量不正常的范德华相互作用。2、整系统能量最优化。去除整个系统中能量不正常的相互作用。3、有限制的分子动力学。保持蛋白质不动，溶解溶剂的不同层，同时逐渐将系统温度从0K提升到300K。4、整系统分子动力学模拟。在一个大气压，300K的环境下整个系统分子动力学模拟。可以得到成果的分子动力学模拟。#############################1、溶剂环境能量最优化。该步骤的配置文件min1.in如下：--------------------------------------------------------------------------------oxytocin:initialminimisationsolvent+ions&cntrlimin=1,maxcyc=1000,ncyc=500,ntb=1,ntr=1,cut=10/Holdtheproteinfixed500.0RES19ENDEND该过程保持肽链不动，其中500.0单位是kcal/mol，表示作用在肽链上使其不动的力。“RES19”表示肽链残基数目，因为我们学习使用的oxytocin有9个残基。模拟命令如下：sander–O–imin1.in–omin1.out–poxy.top–coxy.crd–roxy_min1.rst–refoxy.crd&2、整系统能量最优化。配置文件min2.in如下：--------------------------------------------------------------------------------oxytocin:initialminimisationwholesystem&cntrlimin=1,maxcyc=2500,ncyc=1000,ntb=1,ntr=0,cut=10/--------------------------------------------------------------------------------命令如下：sander–O–imin2.in–omin2.out–poxy.top–coxy_min1.rst–roxy_min2.rst&3、有限制的分子动力