博文
使用AmberTools+ACPYPE+Gaussian创建小分子GAFF力场的拓扑文件
|
使用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:amber14bin(如果已存在Path变量, 只要增加C:amber14bin即可).
如果你要使用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的安装路径中不能含有中文, 输入输出文件的保存路径中也不能使用中文
下面就是我们构建出来的分子结构
Fig.12. 使用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.geps和Lig.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.frcmodparmchk2是原先parmchk的增强版, 可以检查输入分子构型中GAFF的缺失参数, 并生成相应的补充参数文件Lig.frcmod.
5. 使用sleap生成AMBER参数文件及坐标文件编写一个leap.in文本文件, 内容如下:
source leaprc.ff14SBsource leaprc.gaffloadamberparams Lig.frcmodlig=loadmol2 Lig.mol2check ligsaveamberparm Lig.prmtop Lig.inpcrdquit然后运行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一步执行成功.
参考资料
◆本文地址: http://jerkwin.github.io/2015/12/08/使用AmberTools+ACPYPE+Gaussian创建小分子GAFF力场的拓扑文件/, 转载请注明◆
https://blog.sciencenet.cn/blog-548663-942117.html
上一篇：Windows下GROMACS程序的编译
下一篇：力场与拓扑之一：力场的概念分类及发展方法
