在Gaussian中計算IRC的方法和常見問題
在Gaussian中計算IRC的方法和常見問題
文/Sobereva @北京科音
First release: 2018-Jan-4 Last update: 2020-Sep-14
0 前言
IRC是量子化學研究化學反應的重要概念，它是質權坐標下連接勢能面相鄰兩個極小點的能量最低路徑，描述了化學過程在不考慮熱運動因素下最理想的結構變化軌跡，對于討論微觀化學過程至關重要，而且也是驗證過渡態找沒找對的最決定性的方法。在無數年前筆者寫的《過渡態、反應路徑的計算方法及相關問題》（http://www.shanxitv.org/44）中有很多介紹和討論，建議看看。
在Gaussian中計算IRC的問題是筆者在網上被問及最多的問題類型之一，在計算化學公社論壇(http://bbs.keinsci.com)和思想家公社QQ群里屬于半周經問題。為了免得日后反復重復回答大量相似問題，決定寫個文章專門說一下。下文先說一下Gaussian中產生IRC的方法，然后再集中回答一下我常被問及的各種與IRC有關的問題。還有很多其它和做IRC有關的知識，限于時間精力以及表達形式的限制就不寫在這里了，筆者會在北京科音的初級量子化學培訓班（http://www.keinsci.com/workshop/KEQC_content.html）和基礎量子化學培訓班（http://www.keinsci.com/workshop/KBQC_content.html）里詳細講，并且給出極其豐富的研究例子。
筆者還另外寫了一篇和IRC有密切關系的文章《談談Gaussian產生downhill路徑的功能》（http://www.shanxitv.org/571），是Gaussian的IRC功能的一個特殊用法，推薦在看完本文后閱讀。
本文內容適用于Gaussian09和16。
1 在Gaussian中產生IRC的方法
在產生IRC之前，必須先進行過渡態(TS)優化。必須以過渡態的結構為初始結構，才能進行IRC任務。既可以將優化出的TS結構直接寫在IRC輸入文件里，也可以用geom=check關鍵詞從優化TS任務的.chk文件里讀取最后一幀結構。下面是一個最簡單、典型的找TS+走IRC的過程：
(a)在gview里憑借化學直覺將體系擺成盡可能接近于過渡態的結構，保存輸入文件為TS.gjf。
(b)將TS.gjf里的關鍵詞設為# B3LYP/6-31G* opt(calcfc,noeigen,TS)，用Gaussian執行之
(c)打開上一步任務的輸出文件，保存成輸入文件IRC.gjf
(d)將IRC.gjf里的關鍵詞設為# B3LYP/6-31G* IRC=calcfc，用Gaussian運行即可得到IRC。此任務出現報錯幾率很大，仔細耐心閱讀完本文就知道怎么辦了。
IRC任務特別需要注意的一點是，產生IRC和找過渡態用的計算級別必須嚴格一樣！確切來說，任何影響勢能面的設定必須嚴格相同（如果你不知道哪些會產生影響，比較修改前后單點能是否有差異便知）。比如找過渡態用了int=ultrafine scrf(SMD,solvent=ethanol)，那么產生IRC的時候也必須寫上int=ultrafine scrf(SMD,solvent=ethanol)，否則相當于兩個任務的勢能面并不相同，走IRC也就不是從勢能面的準確的一階鞍點來走的，此時IRC任務要么出錯，要么結果無意義。
Gaussian09/16支持多種產生IRC的算法，比較主要的三個在這里說一下：
LQA(Local quadratic approximation)：這是非常傳統的一種IRC產生算法，1988年提出了，每一步都需要利用Hessian矩陣。
HPC(Hessian-based Predictor-Corrector)：2004年由Gaussian作者之一Schlegel提出，這是G09/16默認的算法。這個方法產生IRC每一個點的時候，是將LQA方法作為預測步(predictor)，再用修改的Bulrisch-Stoer方法做為校正步(corrector)。這樣“預測+校正”結合起來，精度比只用LQA更好，而耗時增加不多。此方法校正步的計算過程是一個迭代過程，眾所周知，只要有迭代，就會伴隨著不收斂的可能。HPC方法的一個最大問題，正是做校正步的時候經常不收斂，導致網上問Gaussian的IRC問題的多半都是問這個。
GS2：1989年提出，GS是Gonzalez-Schlegel的縮寫。這是G03默認的算法，在《過渡態、反應路徑的計算方法及相關問題》里有詳細介紹。此方法之所以不再是如今版本默認的，是因為此方法耗時比HPC高很多，每一步都相當于做一次限制性優化。
Gaussian里還支持其它方法，諸如DVV、EulerPC等，由于平時用不到，所以就不提了。
走IRC時，對初始結構需通過calcfc關鍵詞來產生精確的Hessian矩陣。LQA和HPC方法走每一步的時候也都需要Hessian矩陣，默認情況下是通過Bofill方法基于梯度和前一步的Hessian近似產生的。如果你在IRC里用calcall關鍵詞代替calcfc，則每一步的Hessian矩陣都會精確生成，顯然IRC走得也會更準確，但是耗時會暴增一個數量級甚至更多。也可以在寫calcfc的同時用recalc=n關鍵詞，此時每n步才會精確重算一次Hessian矩陣。顯然，n越小越準確，耗時也越高，n=1時和calcall是等價的。
走IRC默認是先往反應路徑一邊走，之后再往另一邊走，兩邊的路徑和初始的TS結構合在一起就是完整的IRC。IRC每一幀結構對應一個反應坐標和一個能量，通常將IRC以“能量vs反應坐標”繪制成曲線圖表示。在Gaussian中IRC任務的初始結構（過渡態結構）被作為反應坐標的零點，往反應物一側反應坐標為負，往產物一側反應坐標為正（但有時也可能是和實際情況反過來的，看后文）。Gaussian給出的反應坐標單位是amu^1/2 Bohr，這里amu是指atomic mass unit，有amu^1/2這項是因為IRC是在質權坐標下定義的。
由于IRC是在質權坐標下定義的，因此結果受同位素質量設定的影響，這里有具體例子展現同位素產生的影響《談談溫度、壓力、同位素設定對量子化學計算結果產生的影響》（http://www.shanxitv.org/423）。如果在IRC里寫上Cartesian關鍵詞，代表在笛卡爾坐標下執行IRC任務，此時結果就不叫IRC了，確切來說叫做MEP（minimal energy path，極小能量路徑）。筆者時不時看到有的菜鳥莫名其妙地在IRC里寫了個cartesian關鍵詞，明明他都不知道這是什么含義就瞎寫，明顯是被以訛傳訛了，而且還是不求甚解就知道盲目模仿。
在IRC里可以用Reverse或Forward關鍵詞讓程序只往逆方向或者正方向一側走IRC。經常有人問諸如怎么我寫的是Reverse，但是Gaussian反倒往產物方向走了？或者，我讓程序往兩邊走，走出來的IRC怎么左邊是產物右邊是反應物？那是因為哪邊是反應物哪邊是產物只有研究者自己知道，Gaussian知道的只是那是兩個不同方向而已。如果你想明確讓Gaussian知道指定哪邊是正方向哪邊是逆方向，需要在IRC里利用Phase關鍵詞定義。
走完IRC后，可以用gview打開輸出文件，在窗口左上角切換幀號可觀看IRC每一幀結構，或者點綠色圓點播放IRC軌跡。IRC軌跡動畫可以通過File - Save Movie - Save Movie File保存出來。IRC的“能量vs反應坐標”圖可以用Results - IRC/Path來觀看，點擊圖的左上角的Plots - Plot Molecular Property還可以把幾何參數等信息隨反應坐標的圖繪制出來，在圖上點右鍵還可以將數據點導出，用來在Origin之類程序里繪制和調節達到更好效果。
走IRC有兩個參數非常關鍵：
(1)Maxpoints：設定的是每個方向走的步數，默認為10。比如設Maxpoints=20且讓IRC往兩邊走，則IRC總共最多走1+2*20=41個點，其中1是過渡態那個點。注意maxpoints設的是步數上限，因為走IRC時，如果已經走到了被程序自動判斷是離極小點很近的位置了，則IRC就認為這個方向的IRC已經產生完畢了，會提示PES minimum detected，然后就開始走另一頭了。如果兩頭都走完了，就輸出匯總信息然后正常結束了。maxpoints設的越大，顯然IRC可能走的點數就越多，總耗時也會越高。
(2)Stepsize：設的是IRC步長，默認為10。如果設的是正值，則單位是0.01 Bohr，如果設的是負值，則單位是0.01 amu^1/2 Bohr。比如你設stepsize=20，則步長就是0.2 Bohr，由于默認情況下IRC任務是在質權坐標下做的，因此程序會自動轉換為以amu^1/2 Bohr為單位的質權坐標下的步長，仔細看輸出文件就能找到相應信息。stepsize設得越小，IRC越準確、曲線越光滑，IRC軌跡震蕩情況越不容易出現；stepsize設得越大，則IRC越不準確，曲線越可能出現褶皺，觀看IRC軌跡時也越可能看到結構震蕩現象，在用HPC方法的時候還越容易出現校正步不收斂而報錯的問題。
為了更直觀看到stepsize對IRC產生的影響，以及LQA和HPC的差異，這里給出HPC原文里的圖，這是對一個模型勢能面Muller-Brown surface走IRC的情況：
由圖可見，相同步長情況下，右圖的HPC方法跑出的IRC的點比左圖的LQA方法更接近于精確的IRC曲線（實線），因為引入了校正步。而帶來的代價就是在實際研究中很容易出現校正步不收斂問題。從上圖也看出，IRC步長設得越小，走的IRC點也越接近精確路徑，反之誤差越大。在LQA方法結合0.2 Bohr步長（默認的stepsize的兩倍）時走出的IRC已經誤差挺明顯了。
maxpoints和stepsize共同決定IRC最多能走多長，即每一側長度上限是maxpoints*stepsize。默認設定下，IRC只能跑出與TS比較接近的一段，如果你想讓IRC跑得更長，顯然要么增加maxpoints要么增大stepsize。增大maxpoints會增加耗時，而增大stepsize，則會導致IRC精度下降、HPC方法容易因校正步不收斂而出錯。
不同反應的反應路徑的總長度是不同的，可能某個設定下對A反應已經把IRC跑得很完整了，但是對于B反應，IRC的兩個端點距離反應物和產物極小點還尚有不小距離。要判斷IRC是否跑得比較完整，可以看“能量vs反應坐標”兩端曲線是否已經接近水平了，也可以看gview在“能量vs反應坐標”圖下方給出的“受力vs反應坐標”圖，如果兩端的點的受力已經比較接近0了，也說明接近極小點了，跑得較完整了。
到底maxpoints和stepsize應該怎么設合適？這要看你跑IRC到底想干嘛，以下是筆者的建議：
(a)驗證過渡態找沒找對：對于這個目的，IRC只要跑一小段就足夠看出趨勢了，也不需要IRC跑得質量多高。maxpoints用默認即可，stepsize可以設大到15或20，從而在不增加計算量的前提下能比默認時稍微跑長一些。建議加上LQA關鍵詞避免HPC方法因校正步不收斂而中斷。
(b)獲得近似的反應物和產物結構：一切同上，但是把maxpoints設50甚至更大，從而使IRC盡可能跑完整。
(c)獲得高質量、完整的IRC曲線：這主要是用于發文章目的，粗糙、不光滑的IRC曲線圖放到文章中肯定會遭人嫌，我們也希望IRC盡量跑完整從而能夠充分描述整個反應歷程。為了讓IRC質量高，可以把stepsize設小到5。由于設小了步長，又想要IRC完整，maxpoints必須設很大，比如200。如果此時跑出來的IRC還是不平滑，或者因HPC方法的校正步不收斂而中途報錯中斷，應再加上calcall重新跑。如果計算能力不足用不起calcall，也可以用比如recalc=3或5。
2 一些與IRC相關的常見問題
看這一節之前應確保已經仔細閱讀、充分理解了上一節的文字。Q：IRC任務報錯啦！末尾提示Maximum number of corrector steps exceded咋辦？
A：這就是前面反復說的HPC方法校正步不收斂。有以下辦法：
(1)用LQA，比如IRC(calcfc,LQA)，由于此時不涉及校正步了，因此100%解決問題。但這容易造成IRC不夠準確、不光滑。
(2)如果你不希望用LQA在避免報錯的同時犧牲IRC精度，則嘗試減小步長（越小避免報錯的幾率越高），比如用IRC(calcfc,stepsize=5)。或者加上calcall，若嫌太昂貴就改用recalc=x（x越小避免報錯的幾率越高，但也越昂貴）。
(3)改用GS2算法，即IRC(calcfc,GS2)，可完全避免以上報錯。耗時比LQA高得多，但精度也比LQA好。
(4)IRC里加上ReCorrect=never，這使得HPC方法不做校正步，故也完全避免以上報錯。此時耗時和LQA相同，但所得IRC精度不如LQA，因此強烈不建議用。
(5)IRC里加上maxcyc=N（N應大于默認的20）來加大HPC校正步迭代次數上限。筆者時常在菜鳥的IRC輸入文件里看到這關鍵詞，筆者強烈不建議用。雖然此方法不是說解決問題的可能性精確為0%，但可能性實在甚微，沒有試的必要。“不收斂就直接加大循環次數上限”是菜鳥最常見的思維方式。
順帶一提，有時筆者看到有人的IRC輸入文件里在IRC里用了tight，這是用來把校正步收斂限設得更嚴的，用這個完全是莫名其妙，也不知道從哪里學來的。本來默認的收斂限下就容易不收斂，居然還給設得更嚴，明顯會導致出現上述報錯的幾率大增。
Q：怎么IRC剛走了幾步就正常結束了？怎么IRC走出來的兩側的曲線是相同的？
A：此問題是繼上一個問題在網上被問得最多的與IRC有關的問題。出現這種問題都是因為優化過渡態時定位準確度不夠。看下圖，當優化出的過渡態位置不準確時，結構就不是在IRC的極大點了，而是稍微偏離一些的紅球的地方
出現這種情況時，往右邊產生IRC能正常產生，但是從紅球位置往左邊產生IRC時，還沒怎么走，程序就發現能量升高了，誤以為IRC已經走到了離極小點很近的位置，于是就不再繼續走了，就正常結束了。還有一種情況，是剛往左邊走IRC，由于體系受力是沖著右邊的，導致馬上轉了個彎就往右邊走了，就呈現了IRC左右兩邊曲線都一樣的結果。
對這個問題，應按照以下方式排查和嘗試解決
(1)先確保初始結構是之前優化TS得到的結構，而且過渡態優化和走IRC都是在嚴格相同級別下進行的。
(2)提高過渡態定位精度。在找過渡態時候用tight，對于DFT再同時結合int=ultrafine（此時產生IRC也必須用int=ultrafine）。如果還不行，優化過渡態時用calcall（或者用諸如recalc=3）。
(3)如果反復嘗試了(2)的方法還是不行，或者你不想嘗試(2)，畢竟會增加很多耗時，那也可以嘗試增大IRC步長，比如20乃至30。由于步長大了，從上圖紅球的位置往左走的時候可能一下子就越過了TS，之后就能正常繼續往左產生IRC了。不過步長大了容易導致HPC校正步不收斂、IRC不準確不平滑等問題，怎么考慮和處理前面已經說了。
另外，出現這種問題還有一種可能是在IRC任務中，基于自動初猜的波函數做SCF后收斂到的波函數與找過渡態任務最終得到的波函數不同，此時相當于IRC任務所在的勢能面和過渡態搜索任務所在的勢能面不同，這也會導致IRC異常，因為類似于違背了前述的走IRC的“任何影響勢能面的設定必須嚴格相同”的這個前提。出現這種情況時，你會發現IRC任務第一次輸出的SCF Done能量和找過渡態最后一步的SCF Done能量明顯不同。為解決此問題，走IRC的時候可以用guess=read關鍵詞，從優化過渡態的chk文件中讀取最后的波函數（并且最好用forward和reverse關鍵詞通過兩個任務分別跑正向和逆向IRC），這樣通常可以確保IRC任務所在的勢能面和優化過渡態時相同。
用SMD溶劑模型時，也可能個別時候由于數值噪音問題出現IRC走幾步就停了的現象。可將優化和走IRC用的溶劑模型都改為IEFPCM再試，說不定能解決。
Q：我的IRC怎么成這樣了？怎么有個點跑到下頭去了？
A：即便IRC任務還在算著，也可以用gview打開輸出文件看看走到第多少步了，檢查是否正常。當IRC還沒有走完的時候就打開輸出文件，或者看的是IRC失敗的任務的輸出文件，往往就會看到中間有一個點突然掉下去的現象。任務正常結束后再看就不會有這個現象。
Q：走出來的IRC最兩端的點是反應物和產物結構么？
A：否。IRC和幾何優化不一樣，幾何優化的步長是隨機應變的，而IRC的步長是固定的。從TS開始走，即便把maxpoints設得非常大，stepsize設得很小，來讓IRC跑得又完整質量又高，也注定不可能恰好有一個IRC點正好落在與反應物或產物對應的極小點的位置。因此，必須對IRC兩端的點進一步做幾何優化才能得到準確的反應物或產物結構。
Q：為什么我看IRC趨勢是對的，我對IRC兩邊的點優化，優化得到的結構卻不是我想要的極小點結構？
A：有兩種可能：
(1)與這個TS直接連通的極小點本來就是那樣。由于你的IRC跑得還太短，在加上你的化學直覺不準，導致并沒有估計對真正極小點結構。可以讓IRC盡量跑完整一些來檢驗。
(2)優化的時候由于走得步子太大，或者Hessian矩陣不準確，導致優化到了其它鄰近的極小點去。此時可以在優化時用較小的步長上限并結合notrust（相關信息看《量子化學計算中幫助幾何優化收斂的常用方法》http://www.shanxitv.org/164），或者結合calcall/recalc=x提供精確Hessian矩陣，或者先把IRC跑得盡量長一些，使兩端結構盡可能接近極小點時再做幾何優化。
Q：柔性掃描和IRC有什么區別？柔性掃描能代替IRC么？
A：二者概念完全不同，不能代替。IRC是從TS開始順著虛頻方向，沿著梯度負方向移動結構走出來的軌跡，所有坐標是同步變化的。而柔性掃描，是令某個（或某些）坐標以特定方式不斷改變，每次都優化其余所有變量所得到的（限制性優化），故所有變量間的變化不是協同的，得到的能量變化曲線與IRC也是大相徑庭的。柔性掃描曲線及其軌跡經常會出現突躍，而不像IRC那樣總是平滑變化。不過，當某個化學過程主要只涉及一個變量發生變化的話，柔性掃描這個變量和走IRC產生的曲線倒是比較相近，此時如果你用的計算級別沒有解析Hessian而不容易走IRC的話，用柔性掃描也未嘗不可（比如當年沒有TDDFT解析Hessian的時代，由于不方便找TS走IRC，對激發態質子轉移過程的研究經常是靠柔性掃描來代替）。
關于柔性掃描，建議閱讀《詳談使用Gaussian做勢能面掃描》（http://www.shanxitv.org/474）。
Q：我的IRC曲線在過渡態位置怎么特別尖/不平滑，貌似不對勁？
A：應當確認初始結構是優化過的過渡態結構，并且找過渡態和走IRC用的計算級別（包括任何能夠影響單點能的因素）嚴格相同。如果你就是這樣做的，那么嘗試增加過渡態的定位精度，比如優化過渡態時用tight收斂限、calcall/recalc=x、增加積分格點精度等。如果這些也都試了，走IRC時用GS2算法，或加上calcall試試。另外也要確保IRC任務第一步收斂到的波函數與過渡態搜索最后一步相同，這在前面已經說過了。
Q：IRC任務如何續算？
在IRC里添加restart關鍵詞即可，例如IRC(restart,maxpoint=15) b3lyp/6-31g(d)。此時%chk應當和之前中斷的IRC任務相同。
Q：我發現之前IRC跑得不夠長，如何在不完全重算的前提下在兩側再增加一些點？
A：比如maxpoints=15的設定下IRC任務已經正常跑完，現在還想兩個方向各增加5個，執行IRC(restart,maxpoint=20)，%chk應當和之前的IRC任務相同。
注意gview觀看restart得到的IRC軌跡要打開chk/fch文件，而不要打開out/log文件，否則看到的只是續跑出來的點。
Q：IRC任務出現L502錯誤怎么辦？
A：這是走IRC過程中出現了SCF不收斂所導致。按照常規解決SCF不收斂的方法嘗試添加適當的關鍵詞解決，見《解決SCF不收斂問題的方法》（http://www.shanxitv.org/61）。除此以外，如果是IRC任務第一次做SCF就沒收斂，而且之前opt任務的chk文件還在，可以做IRC任務時加上guess=read讀取其中收斂的波函數當初猜。如果是IRC中途某個點SCF不收斂，也可以減小IRC步長，因為IRC的每個點會自動沿用上個點收斂的波函數當初猜，顯然二者結構相差越小，那么上個點收斂的波函數就是這個點越好的初猜波函數。
Q：我的IRC看起來比較怪，合理么？怎么邊上還凸起來一塊？（下圖為筆者隨便找的例子）
Q：我的理論方法支持解析梯度但不支持解析Hessian（比如CCSD），怎么跑找TS并且跑IRC？
A：先說此時怎么找TS。有兩類做法：
(1)純依賴于梯度的方法
用CCSD/cc-pVDZ opt(TS,modRedundant,noeigen)或者opt(TS,gediis,noeigen)或QST2/3，這些過渡態搜索關鍵詞都不需要提供初始Hessian
(2)基于半數值Hessian的方法
把初猜TS結構的Hessian矩陣存到chk文件里：# CCSD/cc-pVDZ freq（這步耗時很高）
從chk里讀Hessian矩陣并優化TS：# CCSD/cc-pVDZ opt(TS,noeigen,readfc)
找到TS結構后，使用這個關鍵詞即可生成IRC：CCSD/cc-pVDZ IRC(gradientonly,euler)。這里gradientonly代表使用只依賴于梯度的IRC算法，默認是DVV方法，由于結果爛到爆，所以改成相對好點的Euler方法，但質量也不怎么樣，經常震蕩得很厲害。
Q：我怎么對IRC過程做波函數分析，考察電子結構隨反應過程的變化？
A：參看《產生Gaussian的IRC和SCAN任務每個點的波函數文件的工具》（http://www.shanxitv.org/199）、《通過鍵級曲線和ELF/LOL/RDG等值面動畫研究化學反應過程》（http://www.shanxitv.org/200）。
Q：我之前對過渡態已經做過振動分析了，chk里已經有了過渡態結構下的Hessian矩陣，能否在IRC計算時直接讀取之，避免用calcfc計算Hessian以節約時間？
A：完全可以，對較大體系也推薦這么做。在IRC里加上rcfc關鍵詞即可，此時就不需要寫calcfc了，程序會從%chk文件中直接讀取Hessian矩陣來用。