使用sobMECP程序結合Gaussian程序搜索極小能量交叉點

使用sobMECP程序結合Gaussian程序搜索極小能量交叉點

文/Sobereva @北京科音

First release: 2015-Mar-18    Last update: 2022-May-27

 

本文主要介紹怎么用sob修改版的Harvey的MECP程序（稱sobMECP）結合Gaussian搜索兩種不同自旋多重度的態之間的極小能量交叉點(MECP)。第一節先簡單介紹一下基本理論和算法方面的知識，第二節介紹一下sobMECP程序的特點和使用流程，第三節給出一些具體實例。

1 原理

1.1 MECP的相關理論知識

我們先看看什么情況下兩個態的勢能面才會交叉。設兩個態波函數分別為|1>和|2>，體系哈密頓算符為H，則這樣的兩態模型構成2*2哈密頓矩陣，其四個矩陣元包括兩個對角元H11=<1|H|1>、H22=<2|H|2>，以及兩個等價的耦合矩陣元H12=<1|H|2>=H21=<2|H|1>。讓這兩個態的能量E1、E2相同需要同時滿足兩個條件：(1) H11=H22 (2) H12=0。

這些哈密頓矩陣元的數值都是體系幾何結構的函數。對于兩個自旋多重度不同的態，在不考慮旋軌耦合情況下，由于波函數自旋部分不同，導致耦合矩陣元H12總為0，因此只需要滿足H11=H22就夠了。對于非線性分子，坐標是3N-6維的，我們可以對坐標進行某種變換，使得H11-H22只依賴于某一維坐標，只要這個坐標恰為某個值時就能滿足H11=H22，而不管其它3N-7個坐標的數值是如何的。這也就是說，不考慮旋軌耦合下兩個自旋多重度不同的態的交叉并不是一個特定的點（對應某個具體的結構），而是在3N-7維的超曲面中交叉的。換句話說，能讓兩個態能量相同的結構有無限個，自然不可能都考察，我們研究的只是其中最有意義的一個點，這一般就是極小能量交叉點(MECP, Minimum energy crossing point)，它是指兩個態能量相同的3N-7維超曲面中能量最低的那個結構。但也不一定MECP就只有一個，正如同幾何優化時往往可以得到多個能量極小點，兩個態之間也可以有多個MECP，能找到哪個取決于MECP搜索時的初猜結構是否離它比較近。

一般搜索MECP的時候是不考慮旋軌耦合的。如果在此之后在計算能量時將旋軌耦合算符引入哈密頓，由于此時兩個態之間耦合矩陣元H12=H21不再為0，兩個態就會發生混合，產生A、B兩個新的態，波函數表示為c1*|1>+c2*|2>形式，且能量相應地出現小幅分裂（即E_A≠E_B）。這種情形類似于自旋多重度相同的兩個態之間的避免交叉極小點（Avoided crossing minimum），因為都是H12=H21≠0.

系間竄越(Intersystem crossing)是重要的光化學過程，是指電子激發后，由于不同自旋多重度的態的勢能面之間存在交叉，導致體系經歷這樣的結構時以非輻射方式改變自旋多重度。因此，尋找MECP的結構對光化學研究十分重要。

1.2 MECP結構的搜索方法

搜索MECP結構涉及兩個問題：（1）理論方法與基組（2）優化算法。

搜索MECP可以用的理論方法有兩類，一類是做態平均CASSCF計算（或MS-CASPT2、MRCI等進一步考慮動態相關的），同時計算多個根，并搜索指定的兩個自旋多重度不同的根之間的MECP。另一類是用普通的理論方法，如DFT、MP2、CCSD，分別對兩個不同自旋多重度的態做單點計算，并由此搜索這兩個態的MECP。一般來說基于DFT搜索MECP是首選，便宜、方便，結果也不錯，可以用在頗大的體系。泛函、基組的選擇和常規DFT計算并無區別。如果對結果要得很精確，可以再考慮CCSD或CASSCF。
（PS：這里說的都是搜索自旋多重度不同的態之間的MECP。搜索自旋多重度相同的態之間的MECP或避免交叉極小點只能用CASSCF，但如果能量分裂較大用TDDFT等單參考態方法也湊合，這些與本文無關就不多說了。）

Robb提出的搜索MECP的算法其實很簡單，和普通幾何優化算法思路頗為類似，這也是最常用的MECP搜索方法。定義兩個有效梯度，f和g。f是(E1-E2)^2對核坐標求導，具體寫為E1-E2乘以兩個態的梯度矢量之差。沿著f移動坐標會令兩個態能量差減小。g是E1對核坐標求導，并且投影掉了它在f方向上的分量以使f和g正交。沿著g移動會令第一個態的能量降低。若將f和g考慮到一起，按照f+g作為體系有效受力進行優化，那么在優化過程中就會一邊降低兩個態的能量差一邊降低體系總能量，顯然最后得到的就是MECP結構。優化過程用的具體數值算法和一般幾何優化一樣是基于牛頓法或贗牛頓法，只不過Hessian矩陣也相應地和普通幾何優化不同，詳細公式見Theor. Chem. Acc., 99, 95(1998)。

MECP搜索一般同時使用位移、受力、兩個態能量差作為收斂判斷標準，三者都非常小的時候才算收斂。但如果把能量差必須接近0這條收斂標準去掉，那么上面介紹的MECP搜索算法也可以直接用于搜索避免交叉極小點。

搜索MECP用的初猜就取平衡結構或過渡態結構就行，離預期的MECP結構越近越好，但并不需要很精確。上述MECP搜索算法還是比較穩健的。如果搜索不成功可以再考慮調整初猜結構。

 

2 MECP程序的使用

Harvey的MECP程序基于Gaussian，只能用于Linux平臺，搜索可以在HF、DFT、MP2等級別下進行，運行過程中會生成Gaussian輸入文件、調用Gaussian進行計算，并且讀取受力，對結構按照搜索算法進行位移。優化算法基于L-BFGS贗牛頓法。MECP原版可以在這里下載：http://pan.baidu.com/s/1o682NB8。使用時應引用Harvey對苯基陽離子搜索MECP的文章Theor. Chem. Acc., 99, 95(1998)。值得一提的是MECP搜索算法并非是Harvey提出來的，MECP程序里用的，也即上一節介紹的算法是Robb提出的標準的MECP搜索方法。之所以Harvey的MECP程序知名度較高，是因為它支持的是最常用的Gaussian，而與此同時Gaussian自身僅支持用CASSCF搜索MECP（眾所周知Gaussian的CASSCF不給力，而且CASSCF本身又相對復雜、耗時），這才導致了MECP程序在文獻中用得很多。

筆者發現這MECP程序用起來頗不方便，每處理一個新體系要改好多地方，用法上也有不少別扭的地方。筆者遂對此程序進行了很多修改，使之變得方便好用得多，修改版稱為sobMECP，下載地址：http://www.shanxitv.org/soft/sobMECP.zip。下文討論全都是對于sobMECP而言的。如果你使用了sobMECP發表了文章，請在文中以這種格式引用：Tian Lu, sobMECP program, http://www.shanxitv.org/286 (accessed month day, year)

2.1 sobMECP程序的基本使用流程

將sobMECP解壓到任意一處，進入其中，然后按下列步驟使用之。

1 設定好Input_Header_A和Input_Header_B文件，內容分別是兩個態的分子坐標之前的部分。%chk必須得寫，chk文件具體怎么命名無所謂。force和guess(read)這兩個關鍵詞是必須寫的，其它關鍵詞，包括方法和基組、輔助SCF收斂的選項等，根據實際情況去寫。必須用#n。

如果需要用gen或genECP從分子坐標末尾讀取基組/贗勢定義，則需要在Input_Tail文件中以常規方式寫明基組/贗勢定義，比如
C H 0
cc-pVTZ
****
Cu 0
SDD
****
[空行]
Cu 0
SDD
[空行]
[空行]
如果不需要自定義基組/贗勢，則Input_Tail文件留空即可。

2 在geom里填好體系初始結構，注意末尾空一行，原子名必須用元素的序號代替。例如：
   6   -0.17110831    0.45776462    0.00000000
   8   -0.01984499   -0.74914326    0.00000000
   1   -1.10138660    1.00356716    0.00000000
   1    0.76276200    0.95429984    0.54191098
   1    0.76276200    0.95429984   -0.54191098
[空行]

3 運行./prepare.sh。這將會在當前目錄下編譯產生MECP.x可執行文件、運行進度文件ProgFile、主腳本runMECP.sh，并且清空當前目錄下的零碎文件和JOBS目錄。

4 運行./runfirst.sh，這將會在初始結構下對兩個態進行計算，產生Job0_A.gjf、Job0_B.gjf、Job0_A.log、Job0_B.log以及相應的chk文件。

5 運行./runMECP.sh，即開始搜索MECP，直到全部標準收斂為止，收斂信息會不斷輸出到屏幕上以便監控。隨著搜索的進行，geom文件會不斷更新為當前結構，ab_initio文件內容會不斷更新為當前兩個態的能量和梯度，ProgFile里記錄的進度信息也會不斷更新。

運行過程中產生的各種細節信息都會輸出到當前目錄下ReportFile文件中。包括每一步的結構、兩個態的能量、收斂情況、有效梯度（即前文提到的f+g，但這里f額外乘了個刻度因子）、差值梯度（兩個態受力之差）、平行梯度（前文的g）。

運行過程中產生的每一步的Gaussian輸入輸出文件都在JOBS目錄下。

MECP搜索過程中產生的每一步的結構都會記錄到traj.xyz下作為軌跡文件，可以用VMD程序打開來方便地觀看搜索過程中的結構變化，第一幀對應于初始結構。

如果運行中途中斷，應重新依次執行prepare.sh、runfirst.sh、runMECP.sh，這會基于geom文件中儲存的最后一步結構重新進行搜索。

2.2 sobMECP程序的相關細節

默認情況prepare.sh是調用gfortran編譯器來編譯MECP程序。如果機子里沒有gfortran，可以改用ifort編譯，做法是在prepare.sh里在gfortran前頭寫上#將之注釋掉，而把ifort那行開頭的#去掉。

如果要調整收斂判斷閾值，應修改temp/MECP.f的第300行。總共有5個標準。TDE是能量變化，TDXMax和TDXRMS是位移最大值和方均根，TGMax和TGRMS是受力最大值和方均根。默認情況下，受力和位移判斷標準比Gaussian09幾何優化默認判斷標準松將近一倍。

如果要調整搜索步數上限，應修改temp/runMECP.sh第30行的數字。默認是100步。

有sobMECP用戶反映把優化過程的置信半徑減小五倍，往往令收斂更容易、降低震蕩傾向。我未充分考證這一點，讀者可以試試。做法是在MECP.f里搜索STPMX = 0.1d0，將數值改小。

進行如上修改后需再次運行./prepare.sh方可生效。

extract_energy是Linux下awk工具的腳本文件，用于從Gaussian輸出文件中提取能量。提取不同理論方法輸出的能量需要用不同的腳本。temp目錄下自帶了四種：
extract_energy_SCF：提取HF/DFT/半經驗的
extract_energy_MP2：提取MP2的
extract_energy_CIS：提取CIS的
extract_energy_TD：提取TDHF/TDDFT的（對于G16用戶，此文件里的TD-KS必須手動改為TD-DFT）
當前計算在哪個級別下進行，就把上述哪個文件拷到上一級目錄下改名為extract_energy。程序包里直接帶的extract_energy對應的是extract_energy_SCF。
如果第一個態和第二個態用的理論方法不同，比如第1個態用TDDFT計算S1，第二個態用UKS計算T1，那么應同時準備extract_energy和extract_energy2，前者會被用于提取第一個態的能量，后者會被 用于提取第二個態的能量。

如果用的Gaussian不是09版，需把sub_script中的g09替換為相應命令。筆者只測試了Gaussian 09 D.01，可以完全正常運行，其它版本G09以及G03、G16應該也可以正常使用。

 

3 MECP搜索實例

下面給出一些sobMECP程序使用例子，涉及到的文件都在sobMECP的test文件夾里。

3.1 CH2（卡賓）

這一節尋找CH2的單-三重態間的MECP。

在sobMECP文件夾中編輯Input_Header_A文件成為如下內容，對應單重態計算
%mem=6GB
%nproc=4
%chk=singlet.chk
#n B3LYP/6-311G** force guess(read)

First State

0 1

然后，類似地將Input_Header_B文件寫為下面這樣，對應三重態計算
%mem=6GB
%nproc=4
%chk=triplet.chk
#n B3LYP/6-311G** force guess(read)

Second State

0 3

注意上面兩個文件末尾無空行，后同。

然后編輯geom文件，將下面的初始結構寫進去
 6                  0.00000000    0.00000000    0.13397933
 1                  0.00000000   -0.92611695   -0.40193800
 1                 -0.00000000    0.92611695   -0.40193800
[空行]

然后用chmod +x *將sobMECP下的文件都加上可執行權限。確保Gaussian 09可以在命令行下用g09直接正常調用。之后依次運行
./prepare.sh
./runfirst.sh
./runMECP.sh (運行這步之前可能還得執行一次chmod +x *，因為MECP.x是新編譯出來的)

收斂情況隨著計算不斷在屏幕上輸出，僅用了5步就收斂了，屏幕上顯示5個YES。此時geom文件中的結構就是MECP結構了。從ab_initio文件中可見，兩個態的能量分別是-39.1443410840和-39.1443771072 a.u.，確實基本一致。搜索過程的詳細信息在ReportFile里可看到。

將traj.xyz拖到VMD程序的主窗口里可以看到搜索過程結構的變化，可見一開始H-C-H角度是119.9度，而在最后一幀，即MECP結構下，角度減小為100.1度。

3.2 FeO+

這一節尋找FeO+陽離子的四-六重態間的MECP。對O用6-311G*，對Fe用SDD贗勢。對應四重態的Input_Header_A文件應當為
%chk=A.chk
#n B3LYP/genecp force guess(read)

First State

1 4

對應六重態的Input_Header_B文件應為
%chk=B.chk
#n B3LYP/genecp force guess(read)

Second State

1 6

還要寫一個Input_Tail文件，用來設定要從分子坐標后面讀取的內容，這里就是自定義基組和贗勢信息：
O 0
6-311G*
****
Fe 0
SDD
****
[空行]
Fe 0
SDD
[空行]

然后在geom里寫入初始坐標
26 0.000000 0.000000 0.000000
8 0.000000 0.000000 0.670000
[空行]

然后也是依次運行prepare.sh、runfirst.sh、runMECP.sh。

初始結構中Fe-O長度為0.67埃，從輸出結果可見，最后Fe-O長度變為了1.322埃。對于這樣只有一個幾何變量的體系，顯然也可以自行對兩個態做勢能面掃描，然后擬合成曲線來確定交叉點。

注：在2.0埃附近實際上還有個交叉點，但是在附近區域，當前理論方法下兩個態的勢能曲線幾乎完全平行，此時MECP搜索算法完全失效，因此沒法讓初猜Fe-O長度在2.0附近來尋找這個點。這個交叉點只能通過自行考察勢能曲線獲得。

3.3 C6H5+

苯基陽離子C6H5+正是使用MECP需要引的Theor. Chem. Acc., 99, 95(1998)這篇文章中研究的體系，文中使用不同方法考察了它的單-三重態MECP結構。類似前面的例子，也是先寫Input_Header_A、Input_Header_B、geom文件，然后依次運行prepare.sh、runfirst.sh、runMECP.sh。涉及的文件在sobMECP/test/C6H5+中可找到，這里就不累述了。

這里用的初始結構就是把苯去掉一個氫而已，在B3LYP/6-31G**下經過十幾步就找到了想要的MECP。體系依然是C2v構型，但是環發生了一定變形。

sobMECP的test目錄下還有兩個單-三重態MECP搜索例子，其中H3CO+就是原版MECP程序中自帶的一個例子，Pt_coord是個含鉑配合物，算是個較大體系的例子。

3.4 激發態單重態與三重態的交叉

前面考察的單-三重態的MECP，單重態都是直接用DFT算的基態。但是激發態單重態與三重態的交叉也很重要。比如S1-T1的交叉對于TADF（熱活化延時熒光）的研究是關鍵，和反系間竄越效率問題密切相關。

test目錄下有幾個這樣的例子，C6H5+_A2-3B1、C6H5+_B1-3A2、Pyrrole_A2-3B2，橫杠前是單重態的電子態，橫杠后是三重態電子態。這些任務計算第一單重態激發態用的是TDDFT，算三重態用的和之前一樣是UKS。計算的時候要把相應目錄下extract_energy和extract_energy2都拷到sobMECP目錄下，前者提取第一個態的能量（TDDFT給出的），后者提取第二個態的能量（DFT算的）。注意像這些高對稱結構，不同初始結構下T1電子態的不可約表示經常不同。找S1與三重態的MECP時由于不斷讀取上一步的初猜，所以跟蹤的是初始結構下的T1態（這個態到了MECP結構下，未必是能量最低的三重態了）。而TDDFT計算時，S1的電子態的不可約表示則可能發生改變。

雖然原理上也可以用TDDFT算三重態能量，但之所以這里用UKS來算，一方面是更省時間，另一方面是結果往往更合理。

4 關于MECP處的自由能

有人問怎么計算MECP的自由能。一般沒必要算它的自由能，絕大部分涉及到MECP的文章中也只討論電子能量。但如果非要說怎么算的話，應當是在MECP處計算兩個自旋態的能量簡并的3N-7維空間里的3N-7個振動模式的頻率，這是有意義的，因為在這3N-7維空間中MECP相當于極小點。有了振動頻率，就可以用常規方式計算熱力學量了，比如可以把振動頻率等信息提供給Shermo來計算，見《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》（http://www.shanxitv.org/552）。然而，這3N-7個頻率據我所知沒有簡單現成的工具可以計算。

網上有文章稱在MECP處用Gaussian里的freq=projected關鍵詞可以得到MECP的自由能（此關鍵詞本來目的是對IRC上的點計算垂直于IRC方向的頻率），這明顯不具備普適性。例如用上面的C6H5+的例子，在MECP處對單重態和三重態分別做freq=projected時，自由能分別為-231.188990和-231.186629 Hartree，可見存在顯著差異，即結果依賴于自旋多重度的選取，結果不唯一。因此對于非特殊情況，靠freq=projected關鍵詞是沒法得到MECP的自由能的。其問題在于，從MECP位置可以對兩個自旋態用IRC關鍵詞跑downhill路徑，見《談談Gaussian產生downhill路徑的功能》（http://www.shanxitv.org/571），對于MECP來說freq=projected算的是垂直于downhill路徑的3N-7維空間的振動頻率。而MECP處兩個自旋態的downhill路徑不僅方向往往不同，而且也往往都不垂直于唯一決定H11-H22的方向，即垂直于downhill路徑的子空間不等同于兩個態能量簡并的3N-7維空間，因此靠freq=projected一般是沒法正確計算MECP的自由能的。