使用AmberTools+ACPYPE+Gaussian创建小分子GAFF力场的拓扑文件

类别:    标签: gmx   阅读次数:   版权: (CC) BY-NC-SA

使用AMBER的GAFF力场处理有机小分子有很大的优势, 可惜的是GROMACS没有自带GAFF力场, 所以需要组合使用各种软件组合来实现. 这里我们使用的是AmberTools和ACPYPE.

要想使用AmberTools和ACPYPE创建小分子的GAFF力场, 需要先安装这两个工具. 但Windows下AmberTools的安装并不容易. 在这里我提供一个已经编译好的AmberTools+ACPYPE, 它来自Chimera中所带的amber14, 并增加了我自己编译好的RESP程序与ACPYPE程序, 因为Chimera中并没有包含这两个程序. 具体的整合过程请参看Windows下的AmberTools+RESP+ACPYPE.

点击这里下载AmberTools+ACPYPE, 下载后解压到某一目录(路径中不要包含中文字符), 然后新建环境变量AMBERHOME并将其设置为amber14的路径即可.

具体操作如下:

右键我的电脑->属性->高级->环境变量

在用户变量中新建AMBERHOME, 并其值设为amber14的路径, 如C:/amber14. 注意, 要使用/作为路径分格符.

如果你想在任意目录下都可以使用这些工具, 可以在用户变量中新建Path变量(如果存在就选中编辑它), 将其值设为%Path%; C:\amber14\bin(如果已存在Path变量, 只要增加C:\amber14\bin即可).

如果你要使用Python版本的ACPYPE, 而不是(或不能)使用我编译好的二进制版本, 到Python官方网站下载Python 2.x安装即可.

对电荷的处理, 目前有两种比较合理的方法, AM1-bcc电荷与RESP电荷. AM1-bcc电荷可以使用AmberTools中的antechamber程序直接得到(或使用Chimera软件, 它集成了antechamber). RESP电荷则需借助第三方量子化学程序来得到, 如Gaussian, GAMESS等. 因此, 如果你想使用RESP电荷, 那你还需要安装一种量子化学程序. 由于Gaussian使用方便, 所以使用更广泛些. Windows版本的Gaussian及其自带的GaussView界面安装很容易, 网上资料也很多, 这里不做介绍.

使用AmberTools+ACPYPE+Gaussian创建小分子GAFF力场拓扑文件的整个流程如下:

创建过程的大致步骤是先利用Gaussian得到RESP电荷, 然后利用AmberTools得到AMBER的参数文件. 由于GROMACS和AMBER参数文件的格式很不一样,所以最后需要使用ACPYPE将AMBER参数文件转换为GROMACS可识别的.gro和.top文件. 具体步骤说明如下:

1. 使用GaussView创建分子构型并做初步优化

使用熟悉的分子编辑软件创建分子构型. 可用的软件非常多, 一般只要支持输出为.pdb或.mol2格式即可. GaussView, Chimera, Chem3D, VMD等等都可以. 由于我们需要使用Gaussian计算静电势以用于拟合RESP电荷, 所以使用与Gaussian配套的GaussView来构建分子并准备输入文件更方便些. 因此, 如果你对分子编辑软件还没有形成什么偏好的话, 那我推荐你试一试GaussView.

我这里使用的Gaussian 09 C.01版本, 建议你至少使用这个版本, 或使用更新的D.01版, 不要使用更低的版本, 因为下面的有些操作更低的版本不支持. 详细信息见参考资料中的博文.

构建好分子之后, 一般要做下粗略的优化, 使搭建出来的分子构型看起来更合理一些. 这在GaussView中很容易完成, 只要点击Edit -> Clean就可以了.

注意, Gaussian以及GaussView的安装路径中不能含有中文, 输入输出文件的保存路径中也不能使用中文

下面就是我们构建出来的分子结构


视图: 投影 正交    速度:
模型: 球棍 范德华球 棍状 线框 线型   名称
左键: 转动   滚轮: 缩放   双击: 自动旋转开关   Alt+左键: 移动

Fig.1

2. 使用Gaussian进一步优化分子构型并计算静电势

粗略优化之后, 我们要使用Gaussian进行进一步的优化, 同时计算其静电势, 用于后面拟合RESP电荷. 如果你熟悉Gaussian的计算流程, 可以先保存文件, 然后修改输入文件后在命令行中运行. 如果不熟悉的话, 你可以在GaussView的菜单中操作.

点击Calculate->Gaussian Calculation Setup, 打开输入文件编译界面

先将计算类型修改为优化,

设定优化使用的方法, 基组, 以及体系的电荷, 自旋多重度

标题段可改可不改

Link 0部分可设定计算时所用的内存和核数. Windows下Gaussian最多可使用1 GB内存, 我的电脑是四核的, 这里我设置使用2个核

General部分设置使用二次收敛的SCF方法以确保收敛, 选择忽略对称性, 不将连接信息写到输入文件中. 其实这些设置影响不大, 不改一般也没事.

最后, Add. Inp.部分添加静电势输出文件的名称, 并在``中添加计算静电势的关键词

点击Submit..保存为.gjf文件, 并输出直角坐标和附加输入

打开产生的输入文件, 内容如下:

有关输入文件中关键词的解释, 见参考资料中的博文.

注意, Gaussian输入文件中的空行非常重要, 要严格按图上的格式来. 比如, Lig_ini.gesp前面只能有一个空行, 而Lig.gesp后面至少要有两个空行.

得到了输入文件后, 我们就可以使用Gaussian来运行它了. 你可以使用Gaussian自带的用户界面来运行, 也可以使用命令行.

由于这个分子较大, 在我的机器上, 运行了3个多小时才正常结束.

运行正常结束后, 会产生Lig.chk, Lig.out, Lig_ini.gepsLig.gesp等文件, 我们只需要其中的Lig.gesp文件, 它是优化后构型的静电势文件, 我们就使用它来拟合RESP电荷.

另外, 我们也可以利用antechamber命令来产生运行Gaussian需要的输入文件, 像下面这样

antechamber -i Lig.mol2 -fi mol2 -o Lig.gjf -fo gcrt -pf y -gm "%mem=1GB" -gn "%nproc=2" -nc 1 -gk "#HF/6-31G* SCF=tight Pop=MK iop(6/33=2,6/42=6,6/50=1)"

3. 使用antechamebr拟合RESP电荷

使用上一步得到的Lig.gesp文件, 运行antechamebr拟合RESP电荷

antechamber -i Lig.gesp -fi gesp -o Lig.mol2 -fo mol2 -pf y -c resp

我们只需要Lig.mol2输出文件, 它包含了构型以及RESP电荷, 其他文件ANTECHAMBER*, ATOMTYPE.INF BCCTYPE.INF NEWPDB.PDB PREP.INF, esout, qout, punch等都可删除.

4. 使用parmchk2检查GAFF参数并生成缺失参数文件

使用上一步得到的Lig.mol2文件, 运行parmchk2命令

parmchk2 -i Lig.mol2 -f mol2 -o Lig.frcmod

parmchk2是原先parmchk的增强版, 可以检查输入分子构型中GAFF的缺失参数, 并生成相应的补充参数文件Lig.frcmod.

5. 使用sleap生成AMBER参数文件及坐标文件

编写一个leap.in文本文件, 内容如下:

bash
1
2
3
4
5
6
7
source leaprc.ff14SB
source leaprc.gaff
loadamberparams Lig.frcmod
lig=loadmol2 Lig.mol2
check lig
saveamberparm lig Lig.prmtop Lig.inpcrd
quit

然后运行sleap命令

sleap -f leap.in

这样就拿到了分子的AMBER参数文件Lig.prmtop, 结构文件Lig.inpcrd.

也可直接打开sleap依次执行leap.in文件中的每一行, 只不过麻烦一些.

6. 使用ACPYPE将AMBER文件转换为GROMACS文件

运行ACPYPE命令

acpype -p Lig.prmtop -x Lig.inpcrd -d

这样就得到了GROMACS支持的Lig_GMX.gro, Lig_GMX.top, em.mdp, md.mdp等文件. 一般我们只需要前面两个文件.

如果想将.top文件进行处理生成.itp文件,以便在蛋白质的拓扑文件中包含, 可以除去表头, 改动原子类型, 再除去后面的附加信息.

实际上, 上面的3, 4, 5, 6这几个步骤可以使用ACPYPE一步完成, 但在Windows下由于路径的原因很容易出问题. 像上面这样分开做的话, 一步一步完成的话, 出错了容易定位具体的出错步骤. 如果你对这个过程很熟悉了, 可以将这些流程写在一个脚本中自动执行, 或是研究一下如何使用acpype一步执行成功.

【2016-03-28 补注】

Google上搜索acpype code, 主要有两个版本供下载: GitHub版SF版. 前者是个旧版本, 与后者的主要区别在于生成的拓扑文件中二面角的默认函数类型. 旧版本生成的拓扑文件中, 二面角的默认函数类型为31, 但可通过-r选项改变. 而新版本中, 二面角的默认函数类型为94, 可通过-z选项改变. 由于GAFF力场中二面角的函数类型为94, (文档说明), 所以建议使用新版本的ACPYPE. 我提供的压缩包中包含的是新版本. 有关新版本的说明可参考这里.

参考资料

◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: Windows下的AmberTools+RESP+ACPYPE
后一篇: 力场与拓扑之一:力场的概念分类及发展方法

访问人次(2015年7月 9日起): | 最后更新: 2024-01-20 10:40:28 UTC | 版权所有 © 2008 - 2024 Jerkwin