談談Gaussian產生downhill路徑的功能
談談Gaussian產生downhill路徑的功能
文/Sobereva@北京科音 2020-Sep-14
1 前言
在Gaussian程序的IRC關鍵詞中，有一個叫downhill（下坡）的選項，很多人都沒有留意過，在本文就介紹一下此選項的實際用處。我們一般跑IRC都是從過渡態開始跑的，詳見《在Gaussian中計算IRC的方法和常見問題》（http://www.shanxitv.org/400），而IRC=downhill任務則可以以任意位置作為起點，Gaussian會按照IRC的算法得到一條downhill軌跡。這條軌跡相當于在質權勢能面上體系不斷沿著受力方向行進而走出來的理想軌跡，也可以認為是原子運動速度無窮慢的情況下原子走的真實運動軌跡。利用這個功能，我們可以考察體系從非平衡狀態（不是當前勢能面極小點的狀態）向極小點結構弛豫的真實過程。下面就舉兩個例子說明其用法和意義。
下文涉及到的文件都可以在http://www.shanxitv.org/attach/571/file.zip下載。
2 吡咯在T1激發態勢能面的結構弛豫
吡咯的基態是單重態，在單重態極小點結構下此體系是C2v點群的平面體系。但是在最低的三重態（T1態）下，這個體系就不是平面的了。如果想研究吡咯在T1態勢能面的S0極小點結構的位置開始逐漸弛豫向T1態極小點結構的過程，就可以跑downhill軌跡。相應輸入文件是本文文件包里的pyrrole_downhill.gjf，關鍵詞是# b3lyp/6-31g(d) IRC(calcall,downhill,maxpoints=100)，結構用的是B3LYP/6-31G*下優化的基態極小點結構。這里利用了maxpoints=100是讓downhill軌跡盡可能跑得長一些，從而能走到盡可能接近極小點的位置。由于Gaussian（至少對于G09和G16）默認用的IRC算法HPC很不穩定，特別容易出現校正步不收斂的報錯，所以這里用了calcall來很大程度上避免，同時也能增加生成的路徑的精度。雖然這會造成耗時增加很多，但由于當前體系很小、計算級別又低，所以就無所謂了（眾所周知，用LQA也可以避免校正步不收斂的報錯，但都會造成downhill任務卡住，至少對于G09 D.01和G16 A.03經測試都有這個問題。若改用另一種IRC算法GS2來跑雖然原理上也可以，但我發現對于這個例子會失敗）。
跑出來的downhill軌跡如下。和觀看IRC一樣，將輸出文件載入到GaussView里，用Results - IRC就可以觀看。
可見這個軌跡很好地描繪了體系在T1勢能面下自發弛豫，從平面到彎曲的過程。如果你對軌跡末端的結構用0 3結合opt進一步做優化和振動分析，會發現上圖最后的結構和進一步優化得到的結構差異甚微。
一定要注意downhill路徑和幾何優化路徑并不相同，沒有可替代關系。一方面，幾何優化過程的步長是動態的，而不像一般的IRC或downhill軌跡那樣是固定的；而且幾何優化目的只是去找精確的極小點結構，為達到這個目的不擇手段，并不要求過程本質是沿著質權勢能面上體系受力方向走的，因此能量的變化也明顯不是像downhill任務那樣總是平滑下降的。此外，downhill任務至多只能走到一個很接近于與初始結構最近的極小點結構的位置，但不可能恰好落在這個極小點上，必須進一步做幾何優化才能精確定位之。
3 18碳環在電場下的結構弛豫
在《一篇文章深入揭示外電場對18碳環的超強調控作用》（http://www.shanxitv.org/570）文章中介紹了筆者的關于電場對18碳環各方面影響的研究工作，其中筆者就計算了18碳環在電場下的downhill路徑。相應的輸入文件是本文文件包里的orggeom_field_downhill.gjf，結構用的是沒有外電場下優化的結構，關鍵詞寫的是
#p wb97xd/6-311+G(2d) nosymm field=x+275 IRC(calcfc,ReCorrect=never,downhill,maxpoints=100)
由于這個體系和基組不是很小、計算耗時并不低，所以為了避免出現校正步不收斂報錯，沒有用calcall，而LQA如前所述又沒法用，所以用了ReCorrect=never，這直接去掉了校正步，不僅完全避免了報錯，還能節約一些時間，但代價就是得到的downhill路徑在有的地方不那么平滑，不過實際發現這對當前任務的不良影響不明顯。這個任務給18碳環沿著X方向（平行于碳環的一個方向）加了0.0275 a.u.的外電場，并且為了讓電場加的方向和GaussView里載入gjf文件并打開笛卡爾坐標軸時看到的情況對應，用了nosymm避免Gaussian自動把體系調整到標準朝向。由于加了較強外電場可能導致電子被極化得偏離體系顯著，因此我的研究中加了彌散函數以確保描述精度。
這個任務跑出來的軌跡如下，是讓GaussView把數據導出到文本文件后用Origin繪制的，圖中藍色箭頭是外加電場的方向。圖中綠色曲線是最長軸距離，相當于碳環最遠的兩個碳之間的距離。想得到這個曲線的最簡單的做法是在GaussView的IRC觀看界面左上角點Plots，選Plot Molecular Property，然后選Bond，輸入最遠的兩個碳的序號，然后就會出現相應的距離變化圖，然后在相應圖上點右鍵后選保存數據即可。
這個是相應的動畫：
從碳環的軸距曲線和動畫可見，剛加外電場時，碳環沒有馬上發生整體的變形，而是一些鍵長、個別鍵角發生了調整，之后逐漸地開始變形，在電場的驅動下被約拉越長。變形的內在原因、驅動力是什么，筆者在前述的《一篇文章深入揭示外電場對18碳環的超強調控作用》中從能量和受力角度做了深入的分析，非常建議一看。