使用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

    (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

    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。。。

   后面可以用 parmchk -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/


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM