將Gaussian與Grimme的xtb程序聯用搜索過渡態、產生IRC、做振動分析

將Gaussian與Grimme的xtb程序聯用搜索過渡態、產生IRC、做振動分析

First release: 2018-May-29 Last update: 2022-Sep-15

0 前言

在《盤點Grimme迄今對理論化學的貢獻》（http://www.shanxitv.org/388）一文中筆者曾簡單提到GFN-xTB方法，說白了就是類似于DFTB那種半經驗意味的DFT，精度不錯普適性也好，對與之耗時相仿佛的半經驗方法和DFTB帶來了極大的沖擊。Grimme開發了名為xtb的專門做GFN-xTB計算的程序。xtb程序用過的人都說好，速度很快，普適性挺好，推出不久已經開始有很多人使用了，在《2018年度計算化學公社杯最常用的量子化學程序和DFT泛函投票結果統計》（http://www.shanxitv.org/420）里已經有一些得票率了，而且一些第三方程序已經支持xtb了，比如Multiwfn (http://www.shanxitv.org/multiwfn)可以讀取xtb的振動分析的輸出來繪制紅外光譜、基于xtb產生的波函數做波函數分析和繪制分子軌道，molclus程序可以結合xtb來做團簇構型和分子構象搜索（見http://www.keinsci.com/research/molclus.html）。另外計算化學公社上fhh2626還寫了NAMD與xtb結合做QM/MM的界面（http://bbs.keinsci.com/thread-7583-1-1.html）。

xtb程序目前可以做單點、優化、振動分析等任務，但是對于一般計算化學研究者來說，還希望能夠找過渡態、產生IRC，并且希望振動分析的結果（特別是虛頻模式）可以可視化，但這些xtb程序目前還做不到。好在Gaussian從09開始加入了external關鍵詞，在進行極小點/過渡態優化、IRC、振動分析等任務時，可以從外部文件直接讀入能量、受力、Hessian，而外部文件的這些信息可以用任意程序來產生，當然也包括xtb，不過需要自己寫個接口才能實現。因此Gaussian可以被當做一個“optimizer”來使用，這種用法似于ASE（atomic simulation environment）程序。

本文介紹的這個接口筆者稱為gau_xtb，可以從這里下載：http://www.shanxitv.org/soft/gau_xtb。如果大家的實際研究中使用了此接口，請務必這樣引用：Tian Lu, gau_xtb: A Gaussian interface for xtb code, http://www.shanxitv.org/soft/gau_xtb (accessed month day, year) 。此頁面里的腳本始終兼容目前最新版本xtb。當筆者發現xtb的更新導致老的gau_xtb腳本不能用的時候，筆者會在此網頁里對腳本進行更新。因此如果你發現gau_xtb無法正常調用xtb時，首先應當把xtb和gau_xtb都更新為最新版。

筆者還另寫了一篇文章，《Gaussian與ORCA聯用搜索過渡態、產生IRC、做振動分析》（http://www.shanxitv.org/422），將Gaussian和ORCA聯用也帶來很多好處，有興趣可以看看。

下面先介紹xtb的基本用法，然后介紹gau_xtb接口的原理，之后給出具體例子。如果大家已經會用xtb，對接口的原理和細節不感興趣，只想馬上用起來，直接看第4節就行了。

注：經常有人問怎么用本文提供的接口聯用不成功，還懷疑是本文的接口不支持較新版本的xtb。在此明確說明，本文的接口沒任何問題，至少令G16 C.01結合xtb 6.5.1運行都沒問題。遇到跑自己的任務失敗時，起碼先把本文提供的例子重復一遍，死活不成功的話，在讀懂接口腳本的內容基礎上，讓腳本輸出一些中間信息（比如腳本里用ls命令顯示當前目錄下的文件都有哪些、用cp命令將臨時產生的文件復制出來然后人工進行檢查），總能搞明白原因。另外也不要用一些稀奇古怪的運行環境。

1 xtb簡介

借本文的機會簡單介紹一下xtb的相關知識和使用。在https://xtb-docs.readthedocs.io/en/latest/contents.html可以看到在線手冊。

1.1 xtb的安裝

xtb是開源的，可以在https://github.com/grimme-lab/xtb/下載到源代碼，如果想自行編譯的話參看《Grimme的xtb程序的編譯方法》（http://www.shanxitv.org/521）。作為普通用戶，直接用https://github.com/grimme-lab/xtb/releases里提供的預編譯好的版本即可，文件名為比如xtb-191025.tar.xz這樣的形式，代表2019年10月25日發布的版本。這是Linux版，目前沒有預編譯的Windows版。

創建一個目錄，比如/sob/xtb。將比如xtb-191213.tar.xz放進去，在此目錄下執行tar -xJf xtb-191213.tar.xz目錄解壓之，之后應當會看到此目錄下出現了bin、lib64等目錄。

用gedit或vi等工具編輯~/.bashrc文件，加入以下語句
export PATH=$PATH:/sob/xtb/bin
export XTBPATH=/sob/xtb/share/xtb
export OMP_NUM_THREADS=N
export MKL_NUM_THREADS=N
export OMP_STACKSIZE=1000m
ulimit -s unlimited
其中N是并行計算時使用的CPU核心數，不要超過CPU的物理核心數。

然后保存.bashrc文件，重新進入終端，xtb就可以用了。

注：根據我的諸多測試，發現xtb（至少對于GFN2-xTB真空下計算而言）做優化、動力學的速度在大約12核的時候就封頂了，用更多的核來并行對速度的提升微乎其微，某些情況下反倒速度還稍微變得更慢。因此，如果你的機子的物理核心數很多，比如三、四十核，那么建議把N設為12，同時跑多個任務來充分利用計算資源。如果你的CPU比如只有6核，那么N就設6核就好了。但是，如果你要做振動分析，那么用幾十核還是值得的，比如36核速度可以達到12核的兩倍出頭（這是由于振動分析是基于有限差分做的，原理上并行效率可以做得較高）。

1.2 xtb的基本使用

xtb的輸入文件就是一個xyz文件，這是最常用的記錄分子結構的格式之一，很多程序都可以產生。用Multiwfn產生也可以，可以把fch、pdb、mol、wfn、wfx、molden等Multiwfn能認的格式載入Multiwfn，然后進入主功能100的子功能2，選擇輸出xyz文件。由于這個文件格式非常簡單，比如自己把Gaussian的.gjf文件編輯一下產生也可以。

xtb的詳細使用說明也可以通過xtb -h查看。幾個比較常用的選項如下
-c或--chrg：設定體系凈電荷
-u或--uhf：設定alpha電子數減beta電子數（相當于自旋多重度減1）
-g或--gbsa：使用隱式溶劑模型。目前支持的溶劑有toluene、thf、methanol、h2o、ether、chcl3、acetonitrile、acetone、cs2
--molden：計算結束后產生molden.input，這是Molden輸入文件
--gfn：選擇GFN-xTB理論的版本，可以為0、1、2。如--gfn 0就代表GFN0-xTB。GFN2-xTB物理上最嚴格，多數情況精度最佳，但有時候SCF收斂困難；GFN1-xTB不如GFN2-xTB嚴格，平均精度稍遜一點，但SCF收斂容易（因此明顯更適合SCF難收斂的金屬團簇等情況），耗時也比GFN2-xTB低一些。GFN0-xTB精度最爛但速度也最快，非常適合快速簡單粗暴地搞巨大體系，但對于找過渡態的目的就太糙了而不建議用

常用任務類型：
--sp：計算單點（此為默認，可不寫）
--grad：計算梯度
--opt [級別]：幾何優化。級別默認為normal，更佳的是tight、verytight、extreme
--hess：計算數值Hessian并做振動分析
--ohess [級別]：優化后自動計算Hessian并做振動分析
--md：基于當前結構做分子動力學（目前xtb還支持metadynamics，詳見手冊）
--omd：優化后做分子動力學

例如：
對yoshiko.xyz做真空中的單點計算，電荷為1，自旋多重度為2（alpha比beta電子多1個）：xtb yoshiko.xyz --chrg 1 --uhf 1 -sp
對yohane.xyz做甲苯溶劑下優化和振動分析。體系是默認的中性單重態：xtb yohane.xyz --ohess --gbsa toluene

xtb運行時一方面會在屏幕上輸出信息，同時也會在當前目錄下產生一大堆文件。這些文件的含義在自帶的文檔里有說明。

xtb目前有解析梯度，但只支持數值Hessian。--hess或-ohess任務做完會輸出g98.out和g98_canmode.out。前者是模仿高斯freq輸出格式來輸出頻率、紅外強度、正則坐標。后者沒用。

--opt任務產生的xtbopt.xyz是最后結構的xyz坐標文件，其中第二行是對應的能量。xtbopt.log是含有優化過程每一幀的多幀xyz文件，后綴改為.xyz后就可以拖入VMD查看優化軌跡。

Multiwfn可以載入xtb用--molden產生的molden.input文件做十分豐富的波函數分析，相關知識看《Multiwfn波函數分析程序的意義、功能與用途》（http://www.shanxitv.org/184）、《Multiwfn入門tips》（http://www.shanxitv.org/167）。

Multiwfn載入--hess或--ohess的輸出文件后，進主功能11，選擇IR或Raman，進入界面后選0就可以繪制出相應的光譜，超級容易。更多信息見《使用Multiwfn繪制紅外、拉曼、UV-Vis、ECD、VCD和ROA光譜圖》（http://www.shanxitv.org/224）。

1.3 xtb的控制文件

xtb運行時還可以載入控制文件（xcontrol），見在線手冊。通過控制文件，可以對xtb的運行細節做更多的控制、實現更多的功能。有些設置（比如設置體系凈電荷）既可以通過上述選項來指定，也可以在控制文件里指定，前者的優先級更高。

控制文件的文件名隨意，通過-I指定。比如控制文件名字叫inp，那就可以比如這樣執行：
xtb rei_ayanami.xyz -I inp --molden --chrg 1

控制文件里可以設置很多字段，每個字段通過$開頭，到下一個$結束。例如計算當前體系時alpha電子比beta電子多兩個，并且讓3,19,20,21,22原子的位置在優化過程中被凍住，則控制文件的內容應當為：
$spin=2
$fix
atoms:3,19-22
force constant=0
如果想通過諧振勢限制住這些原子而不是完全凍結住，force constant應設為要用的力常數。

更多xtb細節請看自帶的文檔，在北京科音(http://www.keinsci.com)的高級量子化學培訓班里會對GFN-xTB的原理和xtb的用法做深入全面講授。

2 Gaussian的external功能簡介

關于Gaussian的external功能的使用，詳見http://www.shanxitv.org/g09/k_external.htm。簡單來說，Gaussian的輸入文件里寫上比如external='./xtb.sh'，在計算時候就會以這樣的方式調用當前目錄下的xtb.sh腳本：
xtb.sh layer InputFile OutputFile MsgFile FChkFile MatElFile
各個參數的含義可以看手冊，后五個是文件名，這里我們主要關心的是其中第二個參數InputFile和第三個參數OutputFile。如果查看Gaussian輸出文件，會發現這樣的提示
Running external command "./xtb.sh R"
         input file       "/sob/g09/scratch/Gau-28355.EIn"
         output file      "/sob/g09/scratch/Gau-28355.EOu"
         message file     "/sob/g09/scratch/Gau-28355.EMs"
         fchk file        "/sob/g09/scratch/Gau-28355.EFC"
         mat. el file     "/sob/g09/scratch/Gau-28355.EUF"
其中"/sob/g09/scratch/Gau-28355.EIn"和"/sob/g09/scratch/Gau-28355.EOu"正是分別傳遞給xtb.sh的InputFile和OutputFile的文件名。

InputFile文件是Gaussian產生的，記錄了當前步的信息，格式為：
原子數需要的導數電荷自旋多重度
原子1元素序號 X Y Z MM電荷 MM原子類型
原子2元素序號 X Y Z MM電荷 MM原子類型
...
原子N元素序號 X Y Z MM電荷 MM原子類型
如果當前任務只需要能量信息，“需要的導數”為0；如果需要受力，則為1（比如幾何優化任務）；如果還需要Hessian，則為2（比如freq任務，以及優化或IRC時用了calcfc等情況）。

OutputFile是要在執行xtb.sh時候由這個腳本來生成的，里面記錄能量、受力、Hessian等，按照手冊的要求格式為：

其中諸如4D20.12是數據格式，稍微懂點Fortran就能明白。按照這個格式產生好OutputFile文件，則Gaussian就會從中讀取當前任務需要的信息開展計算。比如如果當前任務只需要能量，則填上能量和偶極矩那一行即可，而如果比如需要Hessian，則整個文件所有信息都得填上。此文件中的偶極矩、極化率對于本文涉及的優化、走IRC、振動頻率計算都是不需要的，可以直接填0（但相應地，涉及到這些信息的Gaussian輸出，比如偶極矩、freq任務的紅外強度等也將都為0）。

要想將Gaussian和xtb聯用，關鍵就是要恰當編寫xtb.sh，使得這個腳本可以基于InputFile里的信息去調用xtb計算出當前結構下的能量、受力、Hessian，并按照Gaussian要求的格式轉化出OutputFile文件。下一節就介紹怎么實現。

3 Gaussian與xtb的接口xtb.sh的編寫

筆者寫好的xtb.sh文件在本文一開始提到的壓縮包里有。這是bash shell腳本，本節解釋一下腳本內容，需要懂得一些Linux命令和shell編程知識才能完全理解（值得一提的是，寫這種腳本并非必須用bash shell，也并非必須是Linux環境才能用external功能。比如在Windows下也完全可以寫成.bat腳本，例如Windows版Gaussian調用NBO6就是通過external來調用NBO6的.bat文件實現的）。

xtb.sh文件里$2、$3分別對應于xtb.sh接收到的第2、第3個參數，也即InputFile和OutputFile文件名。腳本首先用
read atoms derivs charge spin < $2
從InputFile中把原子數、需要的導數、電荷、自旋多重度分別讀到atoms、derivs、charge、spin四個變量里，然后用以下命令構建一個mol.tmp文件
cat >> mol.tmp <<EOF
$atoms

$(sed -n 2,$(($atoms+1))p < $2 | cut -c 1-72)
EOF
這個mol.tmp文件是xyz文件的雛形，還需要做兩個處理才能變成xyz格式文件，一方面是把從InputFile讀過來的元素序號轉化成元素名，另一方面是把讀過來的以Bohr為單位的坐標轉化成埃。為此，筆者寫了個Fortran小程序genxyz.f90，編譯好的可執行文件是壓縮包里的genxyz。這個小程序稍微懂點編程的人都能看懂，就不解釋了。xtb.sh以下兩行就是調用genxyz把當前目錄下的mol.tmp轉化為mol.xyz
./genxyz
rm -f mol.tmp

之后，腳本設置xtb并行運行時的線程數，并根據讀入的自旋多重度，將之減1，算出來-uhf后面的參數。然后根據當前任務需要的導數信息，來判斷在調用xtb時是用-grad還是用-hess
export OMP_NUM_THREADS=4
export MKL_NUM_THREADS=4
uhf=`echo "$spin-1" | bc` #nalpha-nbeta
if [ $derivs == "2" ] ; then
echo "Running: xtb mol.xyz --chrg $charge --uhf $uhf --hess --grad > xtbout"
xtb mol.xyz --chrg $charge --uhf $uhf --hess --grad > xtbout
elif [ $derivs == "1" ] ; then
echo "Running: xtb mol.xyz --chrg $charge --uhf $uhf --grad > xtbout"
xtb mol.xyz --chrg $charge --uhf $uhf --grad > xtbout
fi

xtb程序默認用的是GFN1-xTB理論，如果想用原理上更好的GFN2-xTB，可以將上述調用時的參數里加上--gfn 2來改用之。

最后，xtb.sh如下調用筆者自編的extderi程序產生OutputFile文件。對應的源代碼extderi.f90也很簡單，就不解釋了
./extderi $3 $atoms $derivs
extderi會讀取xtb運行后在當前目錄下產生的gradient、hessian文件，分別提取受力、Hessian信息，并且從xtbout文件中讀取能量，然后輸出到OutputFile中。傳遞給extderi的三個參數$3、$atoms、$derivs分別告訴這個程序要產生的OutputFile文件的文件名是什么、總共多少原子、要讀取/寫入哪些導數信息。

xtb.sh中還用了一些rm -f命令，用來刪除xtb產生的各種文件，確保Gaussian運算后當前目錄不殘留多余的文件。

4 Gaussian與xtb聯用搜索過渡態、做振動分析、產生IRC應用示例

注：以下數據用的是本文剛發布的時候的xtb程序算出來的，結果和最新版本xtb可能不同。

我們這里將Gaussian與xtb聯用，搜索一下下圖所示的環丙基卡賓異構化過程的過渡態

對應的Gaussian輸入文件是本文壓縮包里的TS.gjf，開頭兩行如下
%chk=mol.chk
#P opt(ts,calcfc,noeigen,nomicro) external='./xtb.sh'

可見這里用的是常用的opt=TS方式搜索過渡態，因此輸入文件里的結構應當是這個與實際過渡態比較像的過渡態初猜結構。opt里必須寫nomicro，否則Gaussian在優化的時候會試圖調用分子力學的optimizer去搞，達不到我們的目的。這里我們刻意保留了chk文件，因為這樣的話之后做IRC、freq任務就可以直接用geom=allcheck從chk文件中讀取已經優化好的過渡態結構來計算了。注意對于當前情況，不能在這里直接寫opt freq，必須把opt和freq拆成兩步做才行，否則freq任務會出錯。另外，由于當前任務的能量、導數都調用xtb來算了，因此理論方法和基組就不需要寫了。

我們確保機子里已經裝好Gaussian了，xtb也已經配置好了從而可以直接通過xtb命令調用了，然后把TS.gjf、xtb.sh、genxyz和extderi都放到當前目錄下，把xtb.sh里的OMP_NUM_THREADS和MKL_NUM_THREADS都設為CPU的物理核心數。再運行諸如g09 < TS.gjf |tee TS.out，就開始計算了。

這里特別強調一點，如果你的Gaussian的Default.route里已通過-M-設置為默認并行做Gaussian計算，一定要改為串行計算（如果不想動這個文件就在.gjf文件里設%nproc=1），否則在Gaussian通過external方式調用外部腳本期間會造成很高無意義的CPU資源消耗，導致總耗時增加。

此任務收斂很順利，13步就收斂了。找到的過渡態精度如何？下圖上半部分是將Gaussian+xtb找到的過渡態（白線）與B3LYP/TZVP下找到的過渡態（紅線）放在一起進行對比，已經按照《在VMD中計算RMSD衡量兩個結構間的幾何偏差》（http://www.shanxitv.org/290）文中的做法將兩個結構進行了Align使之最大程度匹配。下圖下半部分是過渡態結構的球棍圖，便于讀者看清楚結構。

從上面的對比來看，GFN-xTB方法當初雖然沒有特意考慮過渡態問題，但是優化過渡態的結果和較好精度的B3LYP/TZVP比，誤差基本可以接受，至少定性正確。大家可以用Gaussian+xtb嘗試用不同初猜搜索過渡態，等搜出來一個看著基本合理的過渡態，再用DFT去進一步優化（當然，Gaussian+xtb不是干這個的唯一選擇，筆者也嘗試了用PM7半經驗方法優化這個過渡態，結果也一樣定性正確，至于和Gaussian+xtb給出的結構孰優孰劣，從相對于B3LYP/TZVP結構的RMSD偏差上看半斤八兩）

接下來再做一下振動分析，看看虛頻情況。雖然xtb直接就能做振動分析給出振動頻率，但是不便于觀看振動模式，而通過Gaussian+xtb聯用，結果就可以直接用gview看了。輸入文件是本文文件包里的freq.gjf，內容只有兩行，為
%chk=mol.chk
#P freq geom=allcheck external='./xtb.sh'
運行之，結果是壓縮包里的freq.out。虛頻模式的正則矢量如下

從振動動畫上看，過渡態確實找對了。虛頻大小是732cm-1，而B3LYP/TZVP下是686.7cm-1，PM7下是843.1cm-1，可見Gaussian+xtb的結果合理，而且誤差比PM7明顯更小。

最后，我們再走一下IRC。輸入文件是本文文件包里的IRC.gjf，內容為
%oldchk=mol.chk
#P IRC(maxpoints=20,calcfc) geom=allcheck external='./xtb.sh'
這里用%oldchk是避免IRC任務改寫之前的chk。用Gaussian執行之。gview看到的IRC如下

可見IRC曲線很光滑，而且所有點的結構都正常，證明Gaussian和xtb聯用很成功。

對于上面這個普通有機反應，xtb和PM7都可以給出定性正確的結果。但對于比較復雜、比較難算的牽扯過渡金屬的反應情況又會如何？筆者用Sc金屬引發C-H活化的反應做了個測試，結果如下。具體來說，筆者是先用肯定穩妥的PBE0-D3(BJ)優化出過渡態（下圖左上角。涉及反應的原子被高亮了），然后用這個結構作為初始結構，再用Gaussian+xtb和PM7分別在相應級別下優化過渡態。為了便于對比，三個級別得到的結構都擺成了相同視角。C1和C2分別指的是苯基和甲基上與Sc成鍵的碳。

可見，Gaussian+xtb很好地維持住了PBE0-D3(BJ)的過渡態結構，而且虛頻、鍵長都很合理，筆者也跑了IRC，也很順利且結果正常。然而在PM7下過渡態嚴重變形，跑了一百多步還在嚴重震蕩，最后一幀結構也完全看不出能收斂到合理過渡態的跡象。從最后的結構看，Sc沒有和兩個茂環配位，反倒是與其中一個碳配位了；另外，Sc和甲基、苯基的碳的鍵長居然變得相同了。這說明，用PM7哪怕試圖做這種反應的預優化也根本沒法用。

雖然本文只考察了兩個體系，但至少證明Gaussian+xtb用來粗略研究化學反應是充分可行的，值得在實際研究中廣泛使用。對于普通有機體系，這種做法比起直接用Gaussian自帶的PM6/PM7優勢不顯著。但碰到略詭異，尤其是涉及過渡金屬的體系，預感半經驗方法連定性正確結果也給不出的時候，則十分建議改用Gaussian+xtb。

最后提醒一下，雖然GFN-xTB往往很不錯，但也別以為它的普適性和精度能和一般的DFT計算抗衡。例如優化Li2，M06-2X/def2-TZVP下結果是2.7064埃，實驗值是2.6729埃，但xtb優化出來是2.2991埃，誤差還是不小的，盡管比PM7優化出來的1.8106埃已經強得多了。

如果你要用本文的做法計算很重的元素，比如Pt，有一點主要注意，即雖然Gaussian與xtb聯用時Gaussian自己并不會去做任何量化計算，但是由于你沒寫基組，Gaussian程序內部默認當前設的基組是STO-3G，程序會進而判斷你當前體系里的所有元素是否對于STO-3G有定義，如果沒有就報錯。對很重的元素STO-3G都是沒有定義的。為了避免這個報錯，解決辦法就是關鍵詞里寫上UGBS，這是一個幾乎涵蓋整個周期表的基組，因此就不會因為當前的基組對元素沒有定義而報錯了。