使用gaussian和antechamber拟合RESP电荷过程

    本来觉得挺简单一操作，谁知道竟还是浪费了些时间在这（被gaussian09坑了），遂记录一下。

    拟合RESP电荷目前知道的方法有使用gaussian和antechamber拟合RESP电荷，以及Multiwfn拟合RESP电荷（很简单也很方便，参考Sob大神写的http://sobereva.com/441）。由于本人喜欢在服务器上搞这些操作，所以在这里使用第一种方案。

    使用gaussian和antechamber拟合RESP电荷的过程大致分为两步：首先通过gaussian计算得到esp电荷，然后使用antechamber拟合resp电荷.

    1，构建分子的结构文件，并存为mol2文件

    2，使用gaussian优化结构

    （1）关键词如下：#p HF/6-31G* SCF Pop=MK iop(6/33=2) iop(6/42=6) iop(6/50=1) opt          ##其中，Pop=MK是生成Antechamber可用的gesp文件的关键词

    （2）并在坐标后面输入两个文件名bcr_ini.gesp和bcr.gesp，（前者为初始结构的RESP电荷，后者为优化后的RESP电荷）如：

 

       这里需要注意：

iop(6/33=2) iop(6/42=6) iop(6/50=1)这几个关键字是要求Gaussian输出RESP Fitting。其中
iop(6/33=2)是进行RESP Fitting并输出到Gaussian的.log文件。
iop(6/42=6)是指定精度（的关键字之一）。
以上两个关键字可以在Gaussian 03及之前的版本中使用。
iop(6/50)=1是Gaussian 09 C.01之后推荐的独立于高斯输出文件的resp文件格式，在antechamber中称为"gesp"，使用时需要在高斯输入文件末尾指定单独的gesp的文件名称。
Gaussian 09B.01（可能还有G09A，没有该版本不知道)“误删”了RESP Fitting的代码，所以以上关键字没一个管用。
Gaussian 09C.01及后续版本恢复了误删的代码并且加上了gesp的代码，所以以上关键字全部可以使用。                               参考（http://bbs.keinsci.com/thread-2019-1-1.html）

………我就是在这里被坑了很长时间，本来一直用g09计算，但是怎么也没有gesp文件，报错为只有resp拟合中心，没有values。 最后换用g16可行。

 

另外，我有一个体系含有一个Fe原子，报错为：L602，GetVDW:  no radius for atom XX atomic number XX.

原因：使用pop=CHELPG 拟合静电势时没有内置相应元素的半径。
解决：pop里用readradii，输入文件末尾写上元素名和指定的半径（一般用范德华半径，可以查得），例如（我这个没有做opt，所以只需要刚开始的bfe_ini.gesp电荷文件就行）：

#p HF/6-31G* SCF Pop=(mk,readradii) iop(6/33=2) iop(6/42=6) iop(6/50=1)

title

1 1
FE   -24.5090   -57.4950    22.0040
C    -25.1075   -58.6543    19.7452
O    -25.2675   -59.0353    18.4992
O    -24.1915   -59.2153    20.4402
O    -25.9125   -57.7653    20.1602
H    -24.0775   -59.4572    18.0504

FE=2.05

bfe_ini.gesp

2022.1.5更新：可使用下面的__link__的方式先使用小基组优化一下结构，再进行计算。

如：142原子的体系，如果直接用HF/6-31G*去优化，时间为18.5个小时；如果使用下面这种方式，时间为1.5+2=3.5个小时！

%nproc=56
%mem=20GB
%chk=pc1.chk
# HF/STO-3G opt=loose            #实测这里换成半经验方法AM1也行，但是PM7会报错

title

0 1
N  114.3100   115.9140   116.8710
C  118.8280   114.5630   119.7800
C  114.9620   114.2630   118.6340
......

H  115.7625   129.4238   133.1536
H  116.4610   130.8773   132.4809
C  116.1996   130.9391   134.6123
H  115.2241   131.3758   134.5616
H  116.2410   130.2528   135.4321
H  116.9283   131.7096   134.7550

--link1--
%nproc=56
%mem=20GB
%oldchk=pc1.chk

%chk=pc1_2.chk
# B3LYP/6-31G* SCF Pop=MK iop(6/33=2) iop(6/42=6) iop(6/50=1) opt geom=check

title

0 1

pc1_ini.gesp

pc1.gesp
##最后这里空一行

 

    3，使用antechamber拟合resp电荷

antechamber -i ligand.gesp -fi gesp -o ligand_resp.mol2 -fo mol2 -pf y -c resp  
或者使用：antechamber -i ligand.out -fi gout -o ligand_resp.mol2 -fo mol2 -pf y -c resp    
    PS：前面两个命令的计算结果一样，因为第二个命令默认拟合最终优化后结构的resp电荷（优化后和优化前拟合出来的resp电荷会有一定的区别）

计算完成后得到ligand_resp.mol2文件中即包含了所有原子的resp电荷。
到这里，RESP电荷结果就出来了！   
      注意：拟合得到的总电荷不一定正好是0，比如有时候会是0.000006.但是一般都是接近0，对结果几乎没影响的，一般不用管。
           但是我是强迫症，如果多了一点，我会从最后一个原子的电荷上减去这一点，保证在tleap过程中总电荷为0。。。
注意！：用这套方案生成的mol2文件的原子名会重新命名，需要用户修改mol2文件中的原子名，与pdb中的原子名对应。

   后面可以用 parmchk2 -i ligand_resp.mol2 -f mol2 -o ligand.frcmod 来生成小分子的键参数。

   

 

参考：http://jerkwin.github.io/2015/12/08/%E4%BD%BF%E7%94%A8AmberTools+ACPYPE+Gaussian%E5%88%9B%E5%BB%BA%E5%B0%8F%E5%88%86%E5%AD%90GAFF%E5%8A%9B%E5%9C%BA%E7%9A%84%E6%8B%93%E6%89%91%E6%96%87%E4%BB%B6/