將Gaussian等程序收斂的波函數作為ORCA的初猜波函數的方法
注:讀者務必使用目前Multiwfn官網上的最新版本,早期版本Multiwfn有某些bug
將Gaussian等程序收斂的波函數作為ORCA的初猜波函數的方法
文/Sobereva@北京科音
First release: 2019-Oct-11 Last update: 2021-Oct-1
1 說明
ORCA是非常強大的量子化學程序,筆者之前也寫過不少相關文章,開發了不少輔助工具,見www.shanxitv.org右側的ORCA分類。ORCA相對于最常用的量子化學程序Gaussian來說,有一個缺點是SCF收斂做得不夠好,很多Gaussian能收斂的情況ORCA難以收斂,而且Gaussian的SCF不收斂的解決方案比較成熟,見《解決SCF不收斂問題的方法》(http://www.shanxitv.org/61)。另外,ORCA在產生初猜波函數方面也沒Gaussian靈活,比如Gaussian用戶可以用GaussView構建片段組合波函數作為初猜,而且利用stable=opt關鍵詞還可以自動優化出穩定的波函數,這對于雙自由基、反鐵磁性耦合體系的研究十分重要,參看《談談片段組合波函數與自旋極化單重態》(http://www.shanxitv.org/82)。雖說ORCA也不是沒法算這類自旋極化單重態體系,利用%SCF FlipSpin也可以實現,但往往明顯不如Gaussian方便。
顯然,如果能讓其它量子化學程序,特別是波函數初猜以及SCF迭代方面做得很好的Gaussian產生的波函數作為ORCA的初猜波函數,就能令ORCA取長補短,從很多惱人的問題中解脫。
Multiwfn(主頁&下載地址:http://www.shanxitv.org/multiwfn)程序支持載入.fch、.molden、GAMESS-US/Firefly輸出文件這些含有基組定義以及軌道信息的格式(在Multiwfn里稱為“基函數信息”),并且可以導出各種含有波函數信息的格式,支持的文件格式詳見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》(http://www.shanxitv.org/379)。Multiwfn支持絕大多數主流量化程序產生的波函數,并可以作為波函數文件格式的轉換器來用。從2019-Oct-10更新的Multiwfn開始,Multiwfn的導出文件功能更是新增了產生.mkl文件的功能。.mkl文件是老版本Molekel的輸入文件,可以被ORCA自帶的orca_2mkl工具轉換成.gbw文件,ORCA計算時可以從.gbw文件中讀取軌道作為初猜波函數。顯然,利用Multiwfn可以直接實現將其它量化程序產生的波函數作為ORCA的初猜波函數的目的。甚至可以說,只要存在一種量化程序有辦法得到當前體系當前級別下收斂的波函數,用ORCA計算也一定能收斂。
有一點要提及的是ORCA是基于球諧型高斯函數做的計算,因此用Multiwfn實現這個目的,其它程序計算時也應當用球諧型高斯函數。不了解這方面的話看《談談5d、6d型d殼層基函數與它們在Gaussian中的標識》(http://www.shanxitv.org/51)。
2 例子:丁烷雙自由基
這里舉個例子,計算丁烷雙自由基C4H8,用以說明如何將Gaussian收斂的波函數作為ORCA的初猜波函數,請大家舉一反三。本例涉及的文件都在http://www.shanxitv.org/attach/517/file.rar里。
本例使用2020-Jul-1更新的Multiwfn,Gaussian使用G16 A.03版,ORCA使用4.2版。
首先用Gaussian運行C4H8.gjf,內容如下。任務是對C4H8在B3LYP/def2-SVP級別下產生最穩定波函數,對此體系也是對稱破缺波函數。使用def2-SVP的時候Gaussian默認用的是球諧型基函數。
%chk=C:\C4H8.chk
# B3LYP/def2SVP stable=opt
[空行]
ub3lyp/6-31g(d) opted
[空行]
0 1
C -0.74400100 1.78566400 0.00000000
H -0.60282700 2.33865300 0.92499500
H -0.60282700 2.33865300 -0.92499500
C -0.74400100 0.30988100 0.00000000
H -1.25452600 -0.08746700 0.88463900
H -1.25452600 -0.08746700 -0.88463900
C 0.74400100 -0.30988100 0.00000000
H 1.25452600 0.08746700 -0.88463900
H 1.25452600 0.08746700 0.88463900
C 0.74400100 -1.78566400 0.00000000
H 0.60282700 -2.33865300 -0.92499500
H 0.60282700 -2.33865300 0.92499500
然后用formchk將chk轉換為fch,載入Multiwfn后依次輸入
100 //其它功能 Part 1
2 //導出文件
9 //導出mkl文件
C4H8.mkl
y //表明產生的mkl文件是給ORCA當初猜用,程序會做特殊處理
現在當前目錄下就有了C4H8.mkl。
在當前目錄下運行orca_2mkl C4H8 -gbw,就把C4H8.mkl轉換為了C4H8.gbw。下面,我們將C4H8.gbw里的軌道作為初猜波函數,用ORCA也在B3LYP/def2-SVP下跑一下這個雙自由基的單點。輸入文件如下,將之命名為C4H8.inp,把它和C4H8.gbw都放在當前目錄下,然后進行計算。運行時ORCA會自動從C4H8.gbw中讀取軌道作為初猜波函數。由于C4H8.gbw里的波函數是非限制性波函數,因此寫了UKS關鍵詞。由于ORCA的B3LYP和Gaussian的定義不同,所以加了/G來和Gaussian一致。此文件里的結構和上面Gaussian輸入文件里的結構精確一致。
! UKS B3LYP/G def2-SVP miniprint nopop
* xyz 0 1
C -0.74400100 1.78566400 0.00000000
H -0.60282700 2.33865300 0.92499500
H -0.60282700 2.33865300 -0.92499500
C -0.74400100 0.30988100 0.00000000
H -1.25452600 -0.08746700 0.88463900
H -1.25452600 -0.08746700 -0.88463900
C 0.74400100 -0.30988100 0.00000000
H 1.25452600 0.08746700 -0.88463900
H 1.25452600 0.08746700 0.88463900
C 0.74400100 -1.78566400 0.00000000
H 0.60282700 -2.33865300 -0.92499500
H 0.60282700 -2.33865300 0.92499500
*
迭代過程信息如下
ITER Energy Delta-E Max-DP RMS-DP [F,P] Damp
*** Starting incremental Fock matrix formation ***
0 -157.0020361024 0.000000000000 0.00099750 0.00002892 0.0004319 0.7000
1 -157.0020389471 -0.000002844760 0.00101224 0.00002810 0.0003439 0.7000
***Turning on DIIS***
2 -157.0020411727 -0.000002225540 0.00299633 0.00007975 0.0002632 0.0000
3 -157.0020467355 -0.000005562801 0.00029108 0.00000721 0.0000572 0.0000
4 -157.0020467759 -0.000000040442 0.00006147 0.00000180 0.0000542 0.0000
**** Energy Check signals convergence ****
*****************************************************
* SUCCESS *
* SCF CONVERGED AFTER 5 CYCLES *
*****************************************************
可見由于用了Gaussian收斂的波函數當做初猜,收斂非常快和順利,第一輪迭代時的能量和密度矩陣變化就已經極小了,之后很快就收斂了。之所以不是一輪就收斂,是因為Gaussian和ORCA用的DFT積分格點有異。如果二者都用的是HF方法的話,ORCA里僅僅一輪就能收斂。
如果gbw和輸入文件不同名,為了能從gbw中讀取初猜波函數,應當寫上moread關鍵詞,然后加上一行諸如%moinp "/sob/Azusa_Nakano.gbw"來指明讀取的gbw文件位置。
為了讓Gaussian收斂的波函數放到ORCA里能夠收斂的概率盡可能高,有以下一些注意事項和建議
(1)確保Gaussian里用的基組和ORCA里精確一致。比如Gaussian里用6-31G系列基組時,默認是用笛卡爾型d基函數,而ORCA總是用球諧型基函數,因此Gaussian計算時要寫5d關鍵詞(不過ORCA用戶極少會用Pople系列基組,所以這無所謂)
(2)讓Gaussian計算時實際坐標與ORCA一致。Gaussian在計算時會自動將輸入朝向擺到標準朝向,因此chk文件最后轉出來的gbw里的笛卡爾坐標(標準朝向的)可能和ORCA計算時的不一樣,因此要么ORCA輸入文件里的坐標就用標準朝向的,要么Gaussian計算時候加上nosymm關鍵詞避免擺到標準朝向下,詳見《談談Gaussian中的對稱性與nosymm關鍵詞的使用》(http://www.shanxitv.org/297)。
(3)Gaussian為了節約電子積分計算的耗時,用Dunning相關一致性基組等(部分)廣義收縮的基組在計算的時候會自動把某些指數重復的primitive GTF給去掉。程序還會把收縮系數很小的primitive GTF給去掉。由于ORCA不自動做這種處理,可能導致個別情況下Gaussian里收斂的波函數拿到ORCA里還是沒法收斂。為避免這個問題,可以在Gaussian計算時加上int=NoBasisTransform關鍵詞避免去掉任何primitive GTF,即基組原始怎么定義的就怎么用。如果你沒有經驗的話,我強烈建議總是帶上這個關鍵詞!
(4)建議帶上IOp(3/32=2)關鍵詞(尤其是帶彌散函數時)避免Gaussian自動去除線性依賴的基函數,以確保和ORCA計算使用的基函數的一致性。
(5)對于DFT計算,如果上面問題都已經考慮了,但放到ORCA里還是不能收斂,可以讓Gaussian和ORCA在計算時都用更好的DFT積分格點精度,并且在ORCA里不開RI。
順帶一提,利用Multiwfn還能將其它量化程序產生的波函數作為GAMESS-US的初猜波函數,因為如果載入的波函數文件含有基函數信息,Multiwfn產生的GAMESS-US輸入文件里可以直接帶初猜信息。Multiwfn也可以將其它量化程序產生的波函數當Gaussian的初猜,因為Multiwfn可以導出fch格式,用Gaussian的unfchk轉成chk后,就可以用guess=read來從中讀取波函數當初猜了。若有不清楚的,參看前述的《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》。
3 chk到gbw的批量轉換腳本
為了讓大家可以更方便地把chk轉換成gbw,筆者寫了一個批量自動轉換的bash小腳本chk2gbw.sh,可以在這里下載:http://www.shanxitv.org/attach/517/chk2gbw.sh。運行前需要先把Gaussian、Multiwfn、ORCA都安裝好,運行過程中會自動調用formchk、Multiwfn和orca_2mkl工具。
運行此腳本后,腳本會對當前目錄下每個chk文件進行轉換。比如有個文件叫popipa.chk,運行此腳本后會得到popipa_Gau.gbw,還會得到popipa.inp,這是ORCA輸入文件,打開一看就會看到里面已經寫了moread關鍵詞和%moinp "popipa_Gau.gbw",坐標和chk里的直接對應,因此用戶現在只需要把定義計算級別的關鍵詞改成實際情況即可開始ORCA計算。可見此腳本方便至極!
運行時輸出信息示例:
Converting B2H6.chk to B2H6.gbw ... (1 of 3)
Converting C2H5F.chk to C2H5F.gbw ... (2 of 3)
Converting CCl4.chk to CCl4.gbw ... (3 of 3)
4 gbw到chk的批量轉換腳本
從2020年8月21日更新的Multiwfn開始,在其examples\scripts目錄下提供了gbw2chk.sh文件,只要運行,就會把當前目錄下的所有gbw文件轉化為同名的chk文件,以便通過Gaussian使用guess=read關鍵詞讀取其中的波函數當初猜。運行前需要先把Gaussian、Multiwfn、ORCA都安裝好。注:如果波函數是UHF/UKS計算得到的,必須用2021年10月1日及以后更新的Multiwfn。