AMBER分子动力学简例.docx

AMBER分子动力学简例AMBER分子动力学简例（一）概述以下是使用AMBER包的简单教程，希望对开始学习分子动力学的同学有用处。申明一下，以下教程原版来自网上，是最最基本的教程，同时也非常实用，有非常好的借鉴意义。AMBER分子动力学程序包是加州圣弗兰西斯科大学（University of California San Francisco,UCSF)的Peter A Kollman和其同事编的，程序很全，现在已经发展到版本9.0。AMBER功能涵盖种类非常多的生物分子，包括蛋白、核算以及药物小分子。软件详细情况请浏览.以下是AMBER软件包中四个主要的大程序：Leap：用于准备分子系统坐标和参数文件，有两个程序：xleap：Xwindows版本的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书籍。建立项目目录，并进入该目录：mkdirmyprojcd myproj下载1NPO.PDB文件到该目录下，使用文本阅读器阅读该PDB文件。从在文件的开头部分可以得到该文件是一个二聚物的PDB文件，删除文件中A、C、和D链对应的文本行，剩下的就是oxytocin的晶体机构了，保存为oxyt.pdb文件。使用Swiss PDB viewer分析以下oxyt.pdb文件，可以得知该pdb文件缺少了一个侧链。如果使用Deep View，它会自动给文件加上侧链。重新保存pdb文件为oxyt.pdb文件。xleap和tleap程序功能是一样的，都是准备分子结构的坐标文件和拓扑文件。xleap启动一个X界面，比较慢。tleap纯文本，要快一些，我们选择tleap：tleap-s-fleaprc.ff03其中leaprc.ff03是AMBER的2003力场文件。力场是一个很重要的文件，定义了分子、原子和残疾等的信息，请查阅有关资料。oxy=loadpdb oxyt.pdb该命令读入oxyt.pdb文件到oxy变量中。bondoxy.1.SG oxy.6.SG连接oxy变量的第一和第六个残疾的SG原子，即是连接二硫键。可以使用check命令检查oxy是否完好，没有特殊情况的话，最后应该是“OK”，表示分子系统状态是好的，可以保存。check oxy接下来就是保存分子系统了：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配置文件中定义。加溶剂之后，要计算系统是否为电中性，使用命令：charge oxy如果现实不是电中性，要使用addions添加反性电荷，如Cl或者Na等。最后保存溶剂环境的系统拓扑结构和坐标文件，并退出Leap程序：saveamberparmoxyoxy.topoxy.crdquit也可以把命令集中到一个，让tleap读取。如将上面的命令集中到以下文件中oxy.leaprc-sourceleaprc.ff03oxy=loadpdb oxyt.pdbbond oxy.1.SGoxy.6.SGcheck oxysaveamberparm oxy oxy_vac.top oxy_vac.crdsolvateoct oxysaveamberparm oxy oxy.top oxy.crdquit将以上两横线之间的命令保存为oxy.leaprc，然后使用以下命令一步搞定：tleap-s-foxy.leaprc到此，分子系统准备已经完成，目录中会多了一个leap.log文件，这是Leap程序的记录文件，它包含了leap程序所做的一切，包括自动的和手动命令。现在已经有了可以进行分子动力学的基本文件了，再编写sander的控制文件，就可以进行模拟了。这些将在后文中介绍。AMBER分子动力学简例（二）分子动力学（1）真空模式真空模式分子动力学模拟将使用NVT系宗分两步进行，即系统能量最优化和分子动力学过程。1、系统能量最优化。我们将使用淬火能量最优化解除系统内原子之间的不正常相互作用，这些原子之间的高能量相互作用如果不消除，可能影响后续的分子动力学过程。因为动力学过程是能量梯度变化的，太高的能量壁垒可能让MD局限在某一个能量局部最小化位置中。Sander程序是分子动力学模拟程序，它的主要功能是能量最优化，动力学模拟和NMR优化计算。我们必须给Sander一个运行的配置文件，使其按照我们的要求进行计算。能量最优化的配置文件如下：min_vac.inoxytocin: initial minimization prior to MD&cntrlimin = 1,maxcyc = 500,ncyc = 250,ntb = 0,igb = 0,cut = 12/-Sander程序的基本命令参数如下：sander O i in o out p prmtop c inpcrd r restrt -ref refc x mdcrd v mdvel e mden inf mdinfo所以我们的命令可以如下：sander O i min_vac.in o min_vac.out p oxy_vac.top c oxy_vac.crd r oxy_vacmin.rst &可以使用more命令或者tail命令查看min_vac.out的输出内容，那是能量最优化的记录。2、分子动力学模拟。Sander程序的配置文件为：md_vac.inoxytocin MD in-vacuo, 12 angstrom cut off, 250 ps&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 i md_vac.in o md_vac.out p oxy_vac.top c oxy_vacmin.rst r oxy_vacmd.rst x oxy_vacmd.mdcrd ref oxy_vacmin.rst inf &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-4 kcal/mol.A。更多参数将在后文中解释。AMBER分子动力学简例（三）分子动力学（2）水环境中的分子动力学模拟溶剂环境中的分子动力学模拟分为以下四步进行：1、溶剂环境能量最优化。这一步保持溶质（蛋白）不变，去除溶剂中能量不正常的范德华相 互作用。2、整系统能量最优化。去除整个系统中能量不正常的相互作用。3、有限制的分子动力学。保持蛋白质不动，溶解溶剂的不同层，同时逐渐将系统温度从0K提升到300K。4、整系统分子动力学模拟。在一个大气压，300K的环境下整个系统分子动力学模拟。可以得到成果的分子动力学模拟。#1、溶剂环境能量最优化。该步骤的配置文件min1.in如下：-oxytocin: initial minimisation solvent + ions&cntrlimin = 1,maxcyc = 1000,ncyc = 500,ntb = 1,ntr = 1,cut = 10/Hold the protein fixed500.0RES 1 9ENDEND该过程保持肽链不动，其中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&cntrlimin = 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&cntrlimin = 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 restraints10.0RES 1 9ENDEND-上述参数解释如下：ntb = 1：表示分子动力学过程保持体积固定。imin = 0：表示模拟过程为分子动力学，不是能量最优化。nstlim = #：表示计算的步数。dt = 0.002：表示步长，单位为ps，0.002表示2fs。temp0 = 300：表示最后系统到达并保持的温度，单位为K。tempi = 100：系统开始时的温度。gamma_ln = 1：表示当ntt3时的碰撞频率，单位为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。在分子动力学中，当ntr1时，作用力只需要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 &进行MD的运行时间一般较长，可以使用程序的并行版本提交集群计算。在主频位2.0GHz的P4单机上，大概需要一个钟头。可以随时查看md1.in文件的程序输出。4、整系统分子动力学模拟。这一步中，我们将进行整个系统的分子动力学模拟，而不对某些特定原子位置进行限制。因为知识一个小例子，我们将只进行250ps的MD计算。配置文件md2.in如下：-oxytocin: 250ps MD&cntrlimin = 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：表示分子动力学过程的压力常数。ntp1：表示系统动力学过程各向同性。taup = 2.0：压力缓解时间，单位为ps。pres1：引用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 &模拟时间较长，在P4单机2.0GHz的上需要7.5个小时。到 此，模拟全部完成，接下来要对得到的数据进行分析。主要数据文件.out文件，包含系统能量、温度，压力等等；.mdcrd文件，是分子动力学轨迹文件， 可以求系统蛋白的RMSD，回转半径等等。数据分析根据不同的研究目的不同而不同，我们将在后文中进行一些简单的分析。使用pdbviewer软件 load structure的时候，假设侧链缺失会 有提示错误信息。直接使用xleap产生该结构的top crd后，再使用ambpdb重新产生structure即可加上侧链mber中生成小分子模板（转）第一步：生成小分子模板蛋白质中各氨基酸残基的力参数是预先存在的，但是很多模拟过程会涉及配体分子，这些有机小分子有很高的多样性，他们的力参数和静电信息不可能预存在库文件中，需要根据需要自己计算生成模板。amber中的antechamber 程序就是生成小分子模板的。生成模板要进行量子化学计算，这一步可以由antechamber中附带的mopac完成，也可以由gaussian完成，这里介绍用gaussian计算的过程。建议在计算前用sybyl软件将小分子预先优化，不要用gaussian优化，大基组从头计算进行几何优化花费的计算时间太长。gaussian计算的输入文件可以用antechamber程序直接生成，生成后去掉其中关于几何优化的参数即可将小分子优化后的结构存储为mol2各式，上传到工作目录，用antechamber程序生成gaussian输入文件，命令如下：antechamber -i 49.mol2 -fi mol2 -o 49.in -fo gzmat这样可以生成49.in文件，下载到windows环境，运行gaussian计算这个文件，如果发现计算时间过长或者内存不足计算中断，可以修改文件选择小一些的基组。获得输出文件49.out之后将它上传到工作目录，再用antechamber生成模板，命令如下：antechamber -i 49.out -fi gout -o 49mod.mol2 -fo mol2 -c resp运行之后就会生成一个新的mol2文件，如果用看图软件打开这个文件会发现，原子的颜色很怪异，这是因为mol2的原子名称不是标准的原子名称，看图软件无法识别。下面一步是检查参数，因为可能会有一些特殊的参数在gaff中不存在需要程序注入，命令如下：parmchk -i 49mod.mol2 -f mol2 -o 49mod这样那些特殊的参数就存在49mod这个文件中了第二步：处理蛋白质文件amber自带的leap程序是处理蛋白质文件的，他可以读入PDB各式的蛋白质文件，根据已有的力场模板为蛋白质赋予键参数和静电参数。PDB格式的文件有时会带有氢原子和孤对电子的信息，但是在这种格式下氢原子和孤对电子的命名不是标准命名，力场模板无法识别这种不标准的命名，因此需要将两者的信息删除ATOM 12 1H ARG A 82 12.412 8.891 34.128 1.00 0.00 H在PDB各式里面，氢原子的信息会在第13或者14列出现H字符，可以应用grep命令删除在13或者14列出现H的行命令如下：grep -v .H 1t4j.pdb xgrep -v .H x 1t4j_noh.pdb除了删除氢和孤对电子，还应该把文件中的结晶水、乙酸等分子删除，这些分子的信息常常集中在文件的尾部，可以直接删除。处理过之后的蛋白质文件，只包括各氨基酸残基和小分子配体的重原子信息，模拟需要的氢原子和水分子将在leap中添加接下来需要进一步整理蛋白质文件，主要关注氨基酸的不同存在形态和小分子的原子名称。半胱氨酸有两种存在形态，有的是两个半胱氨酸通过二硫键相连，有的则是自由的，两种半胱氨酸在力场文件中是用不同的unit来表示的，这相当于是两个完全不同的氨基酸，需要手动更改蛋白质文件中半胱氨酸的名字，桥连的要用CYX，自由的用CYS,可以通过查看晶体的PDB文件来查看以二硫键存在的半胱氨酸残基。组氨酸有若干种质子态，和半胱氨酸一样，也需要查阅文献确定这些质子态，并更改残基名称最后需要修改的是配体分子的原子名，这是工作量最大的事情，仔细观察可以发现，一般PDB文件中配体的各个原子名称，和我们上面通过antechamber 生成的49mod.mol2中原子名称并不一致，这会造成leap无法识别这些原子，无法为这些原子赋予力参数和静电参数，因此需要手动更改蛋白质文件中配体分子的原子名称。进行这一步可以同时用看图软件打开49mod.mol2和蛋白质文件，隐藏蛋白质文件中除配体分子以外的所有结构，旋转两个文件，让他们姿态相近，以方便观察，并且在各自均标识原子名称，然后用文字编辑软件打开蛋白质文件，核对看图软件中两个分子对应的原子名称，手动逐一修改。修改原子名称需要极大的耐心和细心，如果发生错误下一步的操作会无法继续。我现在想到的也只有这个笨办法，如果谁还有别的好法子，欢迎告诉我。现在文件的准备工作都已经完成，该开始正式的模拟了第三步：生成拓扑文件和坐标文件用amber进行分子动力学模拟需要坐标和拓扑文件，坐标文件记录了各个质点所座落的坐标，拓扑文件记录了整个体系各质点之间的链接状况、力参数电荷等信息。这两个文件是由leap 程序生成的amber中有两个leap程序，一个是纯文字界面的tleap，一个是带有图形界面的Xleap。但是amber的图形界面做得很差，用远程登录使用图形界面不仅麻烦而且慢，所以我比较偏爱使用tleap，两个leap的命令是完全一样的，其实用哪一个都无所谓。启动tleap，在shell里输入命令tleap，leap就启动了，shell里会显示-I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/prep to search path.-I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/lib to search path.-I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/parm to search path.-I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/cmd to search path.Welcome to LEaP!(no leaprc in search path)这个是leap的提示符下面要调入库文件。amber是模拟生物分子的好手，主要就是依靠专门为蛋白质多糖核酸量身订做的amber力场，力场的所有参数都存储在库文件里，所以打开leap第一件事便是调入库文件。amber提供了很多种库，这里我们只用到两个库，gaff和02库，输入命令：source leaprc.gaffsource leaprc.ff02之后两个库就调入进来了这时可以用list命令看看库里都有什么：这里面罗列的就是库里面的unit，包括20种氨基酸、糖以及核酸还有一些常见离子的参数下面一步是调入配体分子的模板，首先补全参数，输入命令：loadamberparams 49mod然后读入模板文件,输入命令：MOL = loadmol2 49mod.mol2其中MOL是unit的名字，要保证这个名字和pdb文件中配体的残基名完全一致，否则系统仍然无法识别pdb文件中的小分子现在再输入list命令会发现库里面多了一个unit:那个就是配体分子的模板下面读入pdb文件，输入命令：comp = loadpdb 1t4j_noh.pdb如果输入这个命令之后，屏幕上出现了大量的创建unit或者atom的信息，如下所示，则说明上面一步的pdb文件处理没有做好，还需要重新处理，通常这种情况都发生在配体分子上面，有时则是因为蛋白质中存在特殊残基。Creating new UNIT for residue: FRJ sequence: 1Created a new atom named: O36 within residue: .RCreated a new atom named: S33 within residue: .RCreated a new atom named: O35 within residue: .RCreated a new atom named: N34 within residue: .R如果屏幕仅仅显示添加原子，这说明输入的PDB文件中缺失了部分重原子，leap根据模板自动补全了这些缺失的原子，这种情况不会影响进一步的计算Added missing heavy atom: .R.AAdded missing heavy atom: .R.AAdded missing heavy atom: .R.AAdded missing heavy atom: .R.A根据体系的具体情况，还可能要将成对的CYX残基中的二硫键相连，有时候还会链接其他原子，比如将糖基链接在氨基酸残基上，用bond命令可以完成，命令用法如下：bond comp.35.SG comp.179.SG其中comp是刚才读入的分子名称，35和179是残基序号，SG是CYX残基模板中硫原子的名称，用comp.35.SG这样的语法就可以定位一个原子如果希望进行考虑溶剂效应的动力学模拟，可能还需要为体系加上水，加水有很多种命令，这里只列举一个：solvatebox comp TIP3PBOX 10.0solvatebox命令是说要加上一个方形的周期水箱，comp指要加水的分子，TIP3PBOX是选择的水模板名称，10.0是水箱子的半径有的体系总电荷不为0，为了模拟稳定，需要加入抗衡离子，系统会自动计算体系的静电场分布，在合适的位置加上离子，命令如下：addions comp Na+ 0意思是用钠离子把体系总电荷补平，还可以选择其他库里面有的离子。完成这一步之后就可以输出拓扑文件和坐标文件了，但是为了方便起见，在运行最后一步之前要先把leap里加工好的分子单独存成一个库文件，以后可以直接调用这个库文件，免得重复上面的操作：saveoff comp 1taj.off这样就生成了一个off文件在那里，下面输出拓扑文件和坐标文件saveamberparm comp 1t4j.prmtop 1t4j.inpcrdChecking Unit.Building topology.Building atom parameters.Building bond parameters.Building angle parameters.Building proper torsion parameters.Building improper torsion parameters.total 1 improper torsion appliedBuilding H-Bond parameters.Not Marking per-residue atom chain types.Marking per-residue atom chain types.(Residues lacking connect0/connect1 -these dont have chain types marked:res total affectedCMET 1)(no restraints)quit现在准备好拓扑文件和坐标文件，接下来就可以开始能量优化和动力学模拟了。如果愿意的话还可以用ambpdb这个命令生成一个pdb文件，直观地看一看生成了什么东西，命令如下：snowyowlslocalhost actualamber$ ambpdb -p 1t4j.prmtop kankan.pdb| New format PARM file being parsed.| Version = 1.000 Date = 09/08/06 Time = 16:36:09snowyowlslocalhost actualamber$可以下载之后用看图软件欣赏，如果加了溶剂盒子的话，看的时候可要小心，会恨吓人的样子。第四步：能量优化用amber进行分子动力学模拟需要坐标和拓扑文件，这在上一步已经完成了，分别生成了1t4j.prmtop 和1t4j.inpcrd两个文件，下面一步就是动力学模拟之前的能量优化了。由于我们进行的起始结构来自晶体结构或者同源模建，所以在分子内部存在着一定的张力，能量优化就是在真正的动力学之前，释放这些张力，如果没有这个步骤，在动力学模拟开始之后，整个分子可能会因此散架。能量优化由sander模块完成，运行sander至少需要三个输入文件，分别是分子的拓扑文件，坐标文件以及sander的控制文件。现在分子的拓扑文件和坐标文件已经完成，需要建立输入文件，min_1.inInitial minimisation of our structures&cntrl imin=1, maxcyc=4000, ncyc=2000, cut=10, ntb=1,ntr=1, restraint_wt=0.5 restraintmask=:1-283/文件首行是说明，说明这项任务的基本情况； &cntrl与/之间的部分是模拟的参数其中imin=1表示任务是能量优化，maxcyc=4000表示能量优化共进行4000步，ncyc=2000表示在整个能量优化的4000步中，前 2000步采用最陡下降法，在2000步之后转换为共轭梯度法，如果模拟的时候不希望进行方法的转换，可以再加入另一个关键词NTMIN，如果NTMIN =0则全程使用共轭梯度法，NTMIN=2则全程使用最陡下降法，此外还有=3和=4的选项，分别是xmin法和lmod法，具体情况可以看手册。第二行的cut=10表示非键相互作用的截断值，单位是埃， ntb=1表示使用周期边界条件，这个选项要和前面生成的拓扑文件坐标文件相匹配，如果前面加溶剂时候用的是盒子水，就设置ntb=1，如果加的是层水，那就应该选择ntb=0；ntr=1表示在能量优化的过程中要约束一些原子，约束的是哪些原子呢？后面有。第三行和第四行都是关于约束原子的信息，restraint_wt=0.5限定了约束的力常数，在这里约束原子就是把原子用一根弹簧拉在固定的位置上，一旦原子偏离固定的位置，系统就会给他施加一个回复力，偏离的越远，回复力越大，回复力就是由这个力常数决定的，单位是Kcal/(mol*A)。 restraintmask=:1-283表示约束的是从1到283号残基，在这个分子中，1-283号残基是蛋白质上的氨基酸残基，从284号开始，就都是水了，所以这个命令的意思就是，约束整个蛋白质，自由地优化溶剂分子。因为溶剂分子是前面的tleap自动给加上的，所以一定要额外多关注一些。下面运行sander来执行这个优化：snowyowlsl