利用高斯和Amber生成乙酰辅酶A配体的拓扑文件

最近在做带有乙酰辅酶A为配体的复合物动力学模拟，前期工作摸索了很久，也参考了很多资料，尤其是生成小分子拓扑的过程。这里想记录一下学习步骤。
在利用Amber生成乙酰辅酶A拓扑之前，需要先计算和拟合它的静电势。

打开GaussView软件，导入*.pdb文件，生成计算静电势所需的*.gif文件：

各项关键字的设置可以参考李继存的博文http://blog.sciencenet.cn/blog-548663-942117.html。
熟悉了之后其实可以直接先submit,生成*.gif文件后，打开文本编辑器，在文本中进行关键词的设置和修改：

%nproshared:设置计算的核数，依据自己情况而定。我们服务器用的16核32GB，可能由于乙酰辅酶A体系较大的原因，且含有P和S原子，跑了大概20多个小时
%mem：设置计算的内存
%chk:设置.chk格式输出文件的路径

最后一行的释义依次为：

B3LYP/6-311G**：设置的方法和基组。
em=GD3BJ：表示DFT-D3色散校正
opt：能量最小化，使其收敛
pop=mk iop(6/33=2,6/41=10,6/42=17)：计算静电势的关键词
scrf=(SMD,solvent=Water)：选择隐式溶剂

Lig为设置的标题信息，可自己修改

-4和1分别代表电荷和自旋多重度。普通有机分子、生物分子等，自旋多重度为1。

使用高斯命令行来运行上面生成的.gif文件

g09 < ACE.GJF > ACE.log

运行结束后生成了lig.log文件，默认在运行的工作目录下生成。始终没有生成李继存博客里提到的lig.gesp文件。后来发现利用lig.log也能转成.mol2文件：

antechamber -i lig.log -fi gout -o lig_resp.mol2 -fo mol2 -pf y -c resp

得到mol2文件以后，就可以继续进行生成配体拓扑的步骤了。
首先将mol2文件里的坐标替换成初始坐标（这里解释一下，在未进行替换之前，直接进行模拟发现，配体坐标并不合理。在计算静电势过程中涉及到结构优化，坐标改变导致）

parmchk检查缺失参数：

parmchk -i ACE_resp_alter.mol2 -f mol2 -o ace.frcmod #输入文件为上一步得到的mol2文件

tleap生成amber参数及坐标文件：

tleap -f oldff/leaprc.ff99SB
>source leaprc.gaff
>LIG = loadmol2 ACE_resp_alter.mol2
>list
>check LIG
>loadamberparams ace.frcmod
>check LIG
>saveoff LIG ace.lib
>saveamberparm LIG ace.prmtop ace.inpcrd
>quit

最后生成gro和top文件：

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

参考资料：

1. RESP拟合静电势电荷的原理以及在Multiwfn中的计算：http://sobereva.com/441
2. DFT-D色散校正的使用：http://sobereva.com/210
3. 谈谈量子化学中基组的选择：http://sobereva.com/336
4. 《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》：http://sobereva.com/327
5. Gaussian中有用的IOp一览：http://sobereva.com/93
6. 使用AmberTools+ACPYPE+Gaussian创建小分子GAFF力场的拓扑文件：http://blog.sciencenet.cn/blog-548663-942117.html
7. 如何计算高斯输入文件的自旋多重度：http://blog.sciencenet.cn/blog-485752-1115385.html