amber中生成小分子模板

amber中生成小分子模板

第一步：生成小分子模板

蛋白质中各氨基酸残基的力参数是预先存在的，但是很多模拟过程会涉及配体分子，这些有机小分子有很高的多样性，他们的力参数和静电信息不可能预存在库文件中，需要根据需要自己计算生成模板。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文件中配体的各个原子名称，和我们上面通过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