• 將Gaussian與ORCA聯用搜索過渡態、產生IRC、做振動分析

    將Gaussian與ORCA聯用搜索過渡態、產生IRC、做振動分析

    文/Sobereva @北京科音

    First release: 2018-May-31  Last update: 2022-Jul-13
      

    1 前言

    在《將Gaussian與Grimme的xtb程序聯用搜索過渡態、產生IRC、做振動分析》（http://www.shanxitv.org/421）中筆者介紹了Gaussian的external功能的使用，以及如何利用這個功能，借用xtb程序產生的能量、受力、Hessian來進行過渡態搜索、IRC、振動分析任務。類似地，只要自己寫一個接口，也可以令Gaussian借用知名的、用戶數目僅次于Gaussian的ORCA程序產生的這些信息來做這些任務，本文就介紹具體怎么實現。看本文前一定先把上文看了。ORCA的一些粗淺介紹可以參看《大體系弱相互作用計算的解決之道》（http://www.shanxitv.org/214）。

    這種Gaussian與ORCA聯用的做法可以帶來以下好處：
    (1)可以令Gaussian做上述任務的時候支持ORCA才支持的理論方法，比如PBEh-3c、B97-3c、PWPB95-D3、NEVPT2、MRCI、DLPNO-CCSD(T)等。雖然其中有的理論方法并沒有解析梯度，但原理上，我們可以通過自寫腳本以有限差分方式得到Gaussian所需的梯度信息
    (2)雖然Gaussian支持的DFT泛函很廣，速度也快，但是由于ORCA的RI做得十分出色，借用ORCA來產生能量和導數可以使得上述任務耗時更低，尤其是對于大體系、大基組而言
    (3)雖然ORCA也能直接做幾何優化和找過渡態，但Gaussian在這方面算法上明顯更成熟、更穩健，選項也更豐富、更易用
    (4)ORCA目前沒有IRC功能，和Gaussian聯用使得ORCA中的理論方法也可以用來產生IRC
    (5)使得ORCA產生的Hessian所對應的振動模式能夠通過gview來觀看

    實際上，ORCA有個關鍵詞ExtOpt，即在ORCA優化過程中通過外部文件讀入能量、導數信息，這和Gaussian的external關鍵詞頗為相似，但ORCA的optimizer的算法跟Gaussian比還是有差距，而且也只能用來優化，所以沒太大意義。

    本文使用Gaussian 09 E.01版，ORCA是4.0.1.2版，系統是CentOS 7.2。本文涉及的所有文件都可以在此處下載： http://www.shanxitv.org/attach/422/gau_orca.zip

    如果大家的實際研究中使用了本文的接口，請這樣引用：Tian Lu, gau_orca: A Gaussian interface for ORCA program, http://www.shanxitv.org/422 (accessed month day, year)

     

    2 Gaussian與ORCA聯用的external腳本的寫法

    文件包里的orca.sh就是筆者編寫的接口腳本，Gaussian輸入文件里寫上external='./orca.sh'關鍵詞就代表需要能量和導數信息的時候會調用當前目錄下的orca.sh。下面解釋一下此腳本的內容。

    腳本一開始的部分，需要由用戶自己填寫ORCA運行時候是幾核并行、內存用多大，用什么計算級別，以及數值精度方面的設定，比如SCF收斂限、積分格點精度。為了導數計算準確，腳本默認帶了tightscf。還要設ORCA可執行文件路徑。這些信息都需要用戶使用前根據實際情況修改
    #Set the number of parallel cores and maximum memory utilized by each core (MB)
    nprocs=1
    maxcore=1000
    #Set calculation level
    level="BLYP def2-svp def2/J"
    #Set parameters for numerical aspects
    numset="tightscf"
    #ORCA executable path
    orcapath="/sob/orca/orca"

    之后，用read命令從Gaussian自動產生的InputFile中讀取原子數、需要的導數階數、凈電荷、自旋多重度到相應變量
    read atoms derivs charge spin < $2

    下面的語句判斷當前任務要用的關鍵詞。ORCA里面engrad（即energy+gradient）關鍵詞和Gaussian里的force關鍵詞等價，用來計算當前結構下能量和受力，并輸出到當前目錄下與輸入文件同名的.engrad文件中。如果寫上freq，則會計算Hessian并且輸出到當前目錄下與輸入文件同名的.hess文件中。這倆文件都是文本文件。
    if [ $derivs == "2" ] ; then
    task="engrad freq"
    elif [ $derivs == "1" ] ; then
    task="engrad"
    fi

    之后產生ORCA輸入文件mol.inp。由于從InputFile讀過來的坐標單位是Bohr，所以用了BOHRS關鍵詞。
    #Create ORCA input file
    echo "Generating mol.inp"
    cat >> mol.inp <<EOF
    ! $level $numset $task BOHRS nopop
    %pal nprocs $nprocs end
    %maxcore $maxcore
    * xyz $charge $spin
    $(sed -n 2,$(($atoms+1))p < $2 | cut -c 1-72)
    *
    EOF
    值得一提的是，ORCA在計算上面語句產生的mol.inp的時候，會在當前目錄下產生mol.gbw。ORCA默認是開啟了autostart設定的，即如果計算時候發現當前目錄下已經有了后綴為.gbw的與當前輸入文件同名的文件，就會自動從中讀取波函數作為初猜波函數。整個Gaussian任務在運行期間是把mol.gbw一直保留在當前目錄的，因此如果Gaussian和ORCA聯用做幾何優化，ORCA在計算當前結構的時候會自動讀取上一次運算時候產生的.gbw里的波函數當做初猜，這樣一方面節約了計算時間（畢竟比重新產生的初猜好），另一方面有助于保持所處的電子態的連續性。（稍有Gaussian知識的人都知道，Gaussian做幾何優化的時候，每一步的初猜用的就是上一步收斂的波函數，通過orca.sh將Gaussian與ORCA聯用時也借助gbw文件同樣發揮了這種效果）

    再往后是調用ORCA運行mol.inp
    #Run ORCA
    echo "Running ORCA..."
    $orcapath mol.inp > mol.out
    echo "ORCA running finished!"

    下面通過筆者用Fortran自寫的extorca程序從ORCA產生的.engrad和.hess文件中提取能量以及導數信息，寫入到與$3參數對應的OutputFile文件中
    echo "Extracting data from ORCA outputs via extderi"
    ./extorca $3 $atoms $derivs
    extorca的源程序是extorca.f90，內容并不復雜，就不多說了。
     
    最后，orca.sh用以下命令清空ORCA運行中途產生的各種mol和mol_開頭的文件以保持當前目錄的清潔。但為了保留下來mol.gbw，用rm之前臨時改了個名字。
    mv mol.gbw tmp.gbw -f
    rm mol.* mol_* -f
    mv tmp.gbw mol.gbw -f

    通過Gaussian與ORCA聯用可以留下chk文件，但里面顯然是沒有波函數信息的，因此之后沒法用Multiwfn等程序來看軌道、做波函數分析。不過，由于mol.gbw文件被保留了下來，這里面記錄了ORCA計算最后一個結構時產生的波函數信息，因此可以轉換為.molden格式，然后載入到Multiwfn里看軌道和做波函數分析。怎么把.gbw轉換成.molden看《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。


    3 實例：HCN->CNH異構化的過渡態搜索、振動分析以及IRC生成

    在本文的壓縮包的example目錄里有HCN->CNH異構化的所有輸入輸出文件，用的是RI-BLYP/def2-SVP級別。由于這個體系非常小，用的計算級別也較低，因此實際上完全發揮不出借用ORCA計算能量和梯度的優勢，反倒比純粹用Gaussian還慢。而對于較大體系，ORCA仗著開掛般的RI，就完全不是同樣的光景了。

    首先運行找過渡態的任務，要把TS.gjf、orca.sh、extorca三個文件都放到當前目錄下，并且根據實際情況恰當修改orca.sh開頭的變量，特別是ORCA的可執行文件路徑。TS.gjf中的坐標是初猜的過渡態結構，文件前三行為
    %chk=TS.chk
    %nproc=1
    #P opt(nomicro,calcfc,ts,noeigen) external='./orca.sh'
    用諸如g09 < TS.gjf |tee TS.out運行即可。注意Gaussian用external功能時應當確保以串行方式計算，因此刻意寫了%nproc=1，否則在通過orca.sh調用ORCA期間，Link401會有很高無意義的CPU占用。

    過渡態任務完成后，可以做振動分析檢驗有無虛頻。對應的輸入文件freq.gjf總共就三行，其它信息通過geom=allcheck從找過渡態后產生的TS.chk里讀
    %chk=TS.chk
    %nproc=1
    # freq geom=allcheck external='./orca.sh'

    用gview觀看輸出文件freq.out，可確認有且只有一個虛頻，振動模式正是期望的，虛頻數值為1089.0cm-1。同樣在BLYP/def2-SVP下，完全由Gaussian計算所給出的虛頻為1088.4 cm-1，可見相符極好。

    最后來跑一下IRC，對應的輸入文件IRC.gjf內容也只有以下三行。此任務從TS.chk里讀取過渡態結構，由于用了%oldchk又避免了此任務改寫TS.chk。
    %oldchk=TS.chk
    %nproc=1
    # IRC=calcfc geom=allcheck external='./orca.sh'

    讀者可用gview打開IRC.out繪制IRC曲線圖，可以看到曲線很光滑，而且IRC軌跡正常，證明Gaussian和ORCA聯用非常成功。
    久久精品国产99久久香蕉