AMBER分子动力学简例
AMBER分子动力学简例(一)
概述
以下是使用AMBER包的简单教程,希望对开始学习分子动力学的同学有用处。申明一下,以下教程原版来自网上,是最最基本的教程,同时也非常实用,有非常好的借鉴意义。
AMBER分子动力学程序包是加州圣弗兰西斯科大学(University of California San Francisco,UCSF)的Peter A Kollman和其同事编的,程序很全,现在已经发展到版本9.0。AMBER功能涵盖种类非常多的生物分子,包括蛋白、核算以及药物小分子。软件详细情况请浏览http://amber.scripps.edu.
以下是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书籍。
建立项目目录,并进入该目录: mkdir
myproj
cd 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 -f
leaprc.ff03
其中leaprc.ff03是AMBER的2003力场文件。力场是一个很重要的文件,定义了分子、原子和残疾等的信息,请查阅有关资料。
oxy=loadpdb oxyt.pdb
该命令读入oxyt.pdb文件到oxy变量中。 bond
oxy.1.SG oxy.6.SG
连接oxy变量的第一和第六个残疾的SG原子,即是连接二硫键。可以使用check命令检查oxy是否完好,没有特殊情况的话,最后应该是“OK”,表示分子系统状态是好的,可以保存。
check oxy
接下来就是保存分子系统了: saveamberparm oxy oxy_vac.top
oxy_vac.crd
其中oxy_vac.top和oxy_vac.crd分别是分子系统真空状态下的拓扑和坐标文件。
接下来要给分子系统添加水环境, solvateoct oxy
TIP3PBOX9.0
该命令让Leap程序使用TIP3PBOX水模型,水环境距离分子为9纳米。也可以用solvatebox命令,但是solvatebox添加的水环境是一个立方体,不是八面体。一般要求水表面到蛋白距离要达到8.5纳米,避免MD过程中蛋白冲出水环境,但是水环境越大计算时间将越长。如果MD过程内含PME,那么水箱长度要大于2倍截矩(cutoff),截矩在MD配置文件中定义。加溶剂之后,要计算系统是否为电中性,使用命令:
charge oxy
如果现实不是电中性,要使用addions添加反性电荷,如Cl-或者Na+等。
最后保存溶剂环境的系统拓扑结构和坐标文件,并退出Leap程序: saveamberparm oxy oxy.top
oxy.crd
quit
也可以把命令集中到一个,让tleap读取。如将上面的命令集中到以下文件中
--------------------------------------------------------------------------------
source
leaprc.ff03
oxy=loadpdb oxyt.pdb
bond oxy.1.SG
oxy.6.SG
check oxy
saveamberparm oxy oxy_vac.top oxy_vac.crd
solvateoct oxy
saveamberparm oxy oxy.top oxy.crd
quit将以上两横线之间的命令保存为
oxy.leaprc
到此,分子系统准备已经完成,目录中会多了一个leap.log文件,这是Leap程序的记录文件,它包含了leap程序所做的一切,包括自动的和手动命令。
现在已经有了可以进行分子动力学的基本文件了,再编写sander的控制文件,就可以进行模拟了。这些将在后文中介绍。
AMBER分子动力学简例(二)
分子动力学(1)
真空模式
真空模式分子动力学模拟将使用NVT系宗分两步进行,即系统能量最优化和分子动力学过程。
1、系统能量最优化。我们将使用淬火能量最优化解除系统内原子之间的不正常相互作用,这些原子之间的高能量相互作用如果不消除,可能影响后续的分子动力学过程。因为动力学过程是能量梯度变化的,太高的能量壁垒可能让MD局限在某一个能量局部最小化位置中。
Sander程序是分子动力学模拟程序,它的主要功能是能量最优化,动力学模拟和NMR优化计算。我们必须给Sander一个运行的配置文件,使其按照我们的要求进行计算。能量最优化的配置文件如下:
oxytocin: initial minimization prior to MD
&cntrl
imin = 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 &cntrl
imin = 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 mdvac.info &
MD的计算过程一般比较久,真空相对与溶剂中要快以下,250ps的模拟大概在一个主频为2。0GHz的Linux单机上运行5分钟。
以上配置文件中用到很多参数,这些参数这是sander程序参数的一小部分,以下将相应解释。如果要了解sander的其他参数,请阅读AMBER用户指南。
&ctrl和
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
&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(我们需要引用一个坐标文件做分子动力学过程的比较,我们需要使用
运行命令如下:
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,回转半径等等。数据分析根据不同的研究目的不同而不同,我们将在后文中进行一些简单的分析。
使用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 > x
grep -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.gaff
>source 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: 1
Created a new atom named: O36 within residue: .R
Created a new atom named: S33 within residue: .R
Created a new atom named: O35 within residue: .R
Created a new atom named: N34 within residue: .R
如果屏幕仅仅显示添加原子,这说明输入的PDB文件中缺失了部分重原子,leap根据模板自动补全了这些缺失的原子,这种情况
不会影响进一步的计算
Added missing heavy atom: .R.A
Added missing heavy atom: .R.A
Added missing heavy atom: .R.A
Added 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.0
solvatebox命令是说要加上一个方形的周期水箱,comp指要加水的分子,TIP3PBOX是选择的水模板名称,10.0是水箱子的半径
有的体系总电荷不为0,为了模拟稳定,需要加入抗衡离子,系统会自动计算体系的静电场分布,在合适的位置加上离子,命令如下:
>addions comp Na+ 0
意思是用钠离子把体系总电荷补平,还可以选择其他库里面有的离子。 完成这一步之后就可以输出拓扑文件和坐标文件了,但是为了方便起见,在运行最后一步之前要先把leap里加工好的分子单独存
成一个库文件,以后可以直接调用这个库文件,免得重复上面的操作: >saveoff comp 1taj.off
这样就生成了一个off文件在那里,下面输出拓扑文件和坐标文件 >saveamberparm comp 1t4j.prmtop 1t4j.inpcrd
Checking Unit.
Building topology.
Building atom parameters.
Building bond parameters.
Building angle parameters.
Building proper torsion parameters.
Building improper torsion parameters.
total 1 improper torsion applied
Building H-Bond parameters.
Not Marking per-residue atom chain types.
Marking per-residue atom chain types.
(Residues lacking connect0/connect1 -
these don't have chain types marked:
res total affected
CMET 1
)
(no restraints)
>quit
现在准备好拓扑文件和坐标文件,接下来就可以开始能量优化和动力学模拟了。如果愿意的话还可以用ambpdb这个命令生成一
个pdb文件,直观地看一看生成了什么东西,命令如下:
[snowyowls@localhost actualamber]$ ambpdb -p 1t4j.prmtop
kankan.pdb
| New format PARM file being parsed.
| Version = 1.000 Date = 09/08/06 Time = 16:36:09
[snowyowls@localhost actualamber]$
可以下载之后用看图软件欣赏,如果加了溶剂盒子的话,看的时候可要小心,会恨吓人的样子。
第四步:能量优化
用amber进行分子动力学模拟需要坐标和拓扑文件,这在上一步已经完成了,分别生成了1t4j.prmtop 和1t4j.inpcrd两个文件
,下面一步就是动力学模拟之前的能量优化了。
由于我们进行的起始结构来自晶体结构或者同源模建,所以在分子内部存在着一定的张力,能量优化就是在真正的动力学之前,
释放这些张力,如果没有这个步骤,在动力学模拟开始之后,整个分子可能会因此散架。
能量优化由sander模块完成,运行sander至少需要三个输入文件,分别是分子的拓扑文件,坐标文件以及sander的控制文件。
现在分子的拓扑文件和坐标文件已经完成,需要建立输入文件,min_1.in
Initial 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来执行这个优化:
[snowyowls@localhost actualamber]$ sander -O -i min_1.in -p 1t4j.prmtop -c 1t4j.i
npcrd -ref 1t4j.inpcrd -r 1t4j_min1.rst -o 1t4j_min1.out
命令中,-O表示覆盖所有同名文件,-i min_1.in表示sander的控制文件是min_1.in,-p 1t4j.prmtop表示分子的拓扑文件
,-c 1t4j.inpcrd表示坐标文件,-ref 1t4j.inpcrd是参考坐标文件,只有在控制文件中出现关键词ntr=1的时候才需要给
sander指定-ref文件,这是约束原子的参考坐标,- ref 1t4j.inpcrd就是说以1t4j.inpcrd中的坐标为准进行约束原子的优化。
以上这四个是输入文件。-r 1t4j_min1.rst表示经过模拟之后新的原子坐标会输出到1t4j_min1.rst文件中,-o 1t4j_min1.out
则表示优化过程中的相关信息都会写入到1t4j_min1.out文件中。 运行起这个命令之后,等拿到 1t4j_min1.rst文件就意味着对溶剂的优化已经差不多了,显然下面还需要对蛋白质本身进行优化,
这个优化还要分两步进行,控制文件分别是min_2.in 和min_3.in min_2.in
Initial minimisation of our structures
&cntrl
imin=1, maxcyc=5000, ncyc=2500,
cut=10, ntb=1,ntr=1,
restraint_wt=0.5
restraintmask=':1-283@CA,N,C'
/
在这里发生变化的是约束原子的范围, ':1-283@CA,N,C'表示1-283号残基中名叫CA,N和C的原子,这些原子实际上是蛋白质主链上
的原子,也就是说这一次的优化是约束了蛋白质主链上的原子之后,对溶剂和侧链原子进行自由优化。
min_3.in
Initial minimisation of our structures
&cntrl
imin=1, maxcyc=10000, ncyc=5000,
cut=10, ntb=1,
/
在这里不再进行约束原子的优化了,对整个分子进行全原子优化。 三步优化的命令如下:
[snowyowls@localhost actualamber]$ sander -O -i min_1.in -p 1t4j.prmtop -c 1t4j.inpcrd -ref 1t4j.inpcrd -r 1t4j_min1.rst -o 1t4j_min1.out
[snowyowls@localhost actualamber]$ sander -O -i min_2.in -p 1t4j.prmtop -c 1t4j_min1.rst -ref 1t4j_min1.rst -r 1t4j_min2.rst -o 1t4j_min2.out
[snowyowls@localhost actualamber]$ sander -O -i min_3.in -p 1t4j.prmtop -c 1t4j_min2.rst -r 1t4j_heat0.rst -o 1t4j_min3.out 接下来就是真正的分子动力学模拟了……
转自: