談談諧振頻率校正因子
1 頻率校正因子的基本知識
通常理論計算的頻率是諧振頻率,由于忽視了非諧振效應,而且理論方法本身也有誤差,所以直接基于它計算基頻、零點能、焓、熵都可能有較大誤差。而將頻率校正因子(后文簡稱校正因子)乘在理論計算出的諧振頻率上,基頻以及熱力學校正量就變得較為準確了,往往能減小2、3倍誤差。例如一般雜化泛函算出的基頻誤差能到100多cm-1,而校正后一般只有幾十cm-1的誤差。由于忽略了非諧振效應會高估頻率,因此基頻的校正因子一般都小于1.0。
校正因子分為多種,通常有
(1)用于獲得準確基頻(fundamental)的校正因子。通常對于所有頻率都是相同的校正值,但有些人為了校正得更準確,把頻率劃分成不同范圍(比如高頻和低頻范圍),分別給出不同的校正因子
(2)用于獲得準確諧振頻率的校正因子。注意諧振頻率不是實驗可觀測的而只是個概念,但可以通過一定手段將實驗數據轉化得到
(3)用于獲得準確的ZPE的校正因子
(4)用于獲得準確的焓變ΔH=H(T)-H(0)的校正因子。注:H(0)=U(0)=電子能量+ZPE
(5)用于獲得準確的熵的校正因子
這5種校正因子各不相同,不要隨意混用。不過(1)/(3)以及(2)/(3)的數值對各種情況倒是基本總是固定的,分別為0.974和1.014。對HF或DFT,(4)和(5)的數值比較接近,混用也無妨。ΔH和熵的數值是依賴于計算時設定的溫度的,其校正因子也因此受溫度影響。
由于諧振子模型下ZPE是頻率的線性函數,即ZPE=∑hν/2,所以第(3)類校正因子乘到原始頻率上還是乘到基于原始頻率算出的ZPE上效果都一樣。但是對于(4)和(5)類校正因子,要分別乘到原始頻率上再算ΔH和熵,而切不要用原始頻率算完了之后再對所得的ΔH和S乘上校正因子。
還有一點要說明,由于ZPE=∑hν/2,表面上看仿佛只要對振動頻率ν乘上校正因子成為準確基頻后ZPE就應該準確了,因此(1)和(2)類校正因子仿佛應當是相同的,這是誤解!因為這種ZPE的計算公式是對于諧振模型而言的,既然校正成了實際的非諧振頻率,那么當然ZPE也得按照非諧振子來算才行。嚴格的非諧振的ZPE必須直接做非諧振計算才能得到。ZPE的校正因子則實際上是先乘到諧振頻率上將之轉化成一種特殊頻率,然后再按照ZPE=∑hν/2來算的時候,得到的結果就能正好和實際的ZPE很接近。
2 頻率校正因子的獲取
校正因子一般都是用先人擬合好的,通過向相應類型實驗數據擬合來得到。用于擬合的分子少則幾種,多則上百種。常用的校正因子可以在這些地方得到:
1 http://comp.chem.umn.edu/freqscale/version3b2.htm。這是Truhlar組擬合的一大堆校正因子,相當全面,包括前述的(1)、(2)、(3)類校正因子。不過很遺憾的是沒有給出ΔH和熵的校正因子。
2 http://cccbdb.nist.gov/vibscalejust.asp。這是Computational Chemistry Comparison and Benchmark DataBase(CCCBDB)數據庫中的一頁,給出的是基頻的校正因子。
3 J. Phys. Chem., 100, 16502-16513 (1996)中的表10,給出的是前述(1)、(3)、(4)、(5)類校正因子。注意表10的ΔH和S的校正因子是298.15K下的,其它溫度下的見文中圖2。此文作者后來在2007年又發了一篇文章J. Phys. Chem. A, 111, 11683-11700 (2007),給出了更多方法和基組(主要是pople基組)下的校正因子。
4 J. Phys. Chem. A, 108, 9213-9217 (2004)。給出的是HF、MP2和B3LYP結合Dunning相關一致性基組的各類校正因子。
5 J. Comp. Chem., 33, 2380 (2012)。給出了一些泛函結合pc系列基組的高頻、低頻、焓、熵、ZPE的校正因子。
6 Theor. Chem. Acc., 105, 413 (2001)。給出了HF、幾種泛函和MP2結合Sadlej pVTZ基組時的高頻和低頻校正因子
網友Liyuanhe整理了一個很好的表格,匯總了各種來源(包括上述的)的校正因子,查閱方便,值得下載留存,見http://bbs.keinsci.com/forum.php?mod=viewthread&tid=3805。
雖然這些來源擬合校正因子過程所用的測試集和方法有異,但對于同種級別,結果是接近的,因此它們是有可比性的。還有一些校正因子零零碎碎地在一些文獻中報道,用得較少,可靠性難說。
對于HF和各種GGA泛函,基頻的校正因子分別接近0.9和1.0。這即是說GGA泛函由于理論方法誤差和非諧振效應抵消,不做校正算出來的基頻也挺準確。而HF由于顯著高估力常數,所以算出來的頻率太高,算出來的3000cm-1實際上只有約2700cm-1而已。雜化泛函的基頻校正因子則介于0.9~1.0之間,而且數值和HF成分是近似呈線性反比關系的。
3 關于頻率計算級別的選擇
注意振動頻率和熱力學校正量的計算絕對不是用的理論級別越高越好。CCSD(T)/cc-pVTZ可能還不如B3LYP/6-31G*算出來的頻率準確,因為誤差既來自于理論方法給出的力常數的誤差,也來自于忽略了非諧振效應的誤差。在某些低水平下,可能恰好這兩部分抵消得好,導致結果不錯,GGA泛函算的基頻幾乎不需要校正就是這個原因。而如果盲目地用很高的計算級別,雖然基本消除了理論方法本身的誤差,但非諧振效應缺乏較好的抵消,那么所得的基頻以及熱力學校正量很可能比很爛的級別算出來的還糟。另外還有關鍵的一點,就是在使用了校正因子后,不同方法間的差距會很大程度上拉平,原本差的方法在校正后結果可能很好。這正是為什么熱力學組合方法G3,雖然目標是QCISD(T)的精度,卻采用HF/6-31G*(經過校正因子校正)這樣表面看上去如此爛的級別來算ZPE。所以說,用很高級別的方法去做振動分析是⑨的行為,何況計算二階導數還那么費時,純屬白耽誤功夫。半經驗方法也有對應的校正因子,但是即便是經過校正后,誤差還是較大,大于從頭算方法很多,所以除非迫不得已不宜用半經驗方法算頻率。基于分子力場也可以算頻率和熱力學校正量,如果分子力場選得合適,那么結果可能很不錯,且不需要做校正,例如MMFF94和MM3對于典型有機分子基頻平均誤差在60cm-1左右。
通常來說,使用B3LYP/6-31G*結合頻率校正因子就能得到很好的計算結果,計算量小,而且比起其它多數方法(哪怕是更昂貴的)經過校正后的結果更好,所以完全不需要考慮其它的。雖然Truhlar的文章表明VSXC經過校正后更準確一點,但畢竟是相對小眾向泛函,寡人不推薦。BMC-CCSD經過校正后特別準確,但畢竟比較貴,除非是要做高精度計算否則就不推薦了。
雖然頻率校正因子校正后的B3LYP/6-31G*是通常推薦的做振動分析的級別,但是有些體系,可能在B3LYP/6-31G*下面優化的結果很糟,諸如弱相互作用顯著的情況,進而可能由此導致算出的頻率和熱力學校正量都很差(即便已經經過校正)。此時就應該選擇更合適的方法做振動分析以及對應的優化。例如對于弱相互作用對結構影響顯著的時候可以用M06-2X、wB97XD。一個有點麻煩的問題是諸如M06-2X、wB97XD這樣的稍微較新的泛函缺乏ΔH和熵的校正因子,不過從JPC,100,16502的數據來看,雜化泛函的這兩類校正因子都比較接近于1.0,所以不做校正也沒什么問題。另外,也可以在B3LYP優化和振動分析時加上DFT-D校正解決這個問題,看《亂談DFT-D》(http://www.shanxitv.org/83)、《DFT-D色散校正的使用》(http://www.shanxitv.org/210)、《談談“計算時是否需要加DFT-D3色散校正?”》(http://www.shanxitv.org/413)。
各類校正因子對于所用基組都不敏感(除非是和很小的基組相比),所以如果你用的基組沒有對應的校正因子,可以直接用此方法在相近的基組下的校正因子。特別值得一提的是彌散函數對校正因子的影響微乎其微,所以假設你的基組用的是6-31+G*,但是找不到相應的校正因子,卻能找到6-31G*下的校正因子,那么就可以放心地用6-31G*的校正因子。另外,對于同一種方法,無論是校正前還是校正后,都不代表越大的基組下結果越準確,而通常是差不多(除非基組爛到3-21G這種層次),因此根本沒必要用雙zeta以上的基組、帶彌散函數做振動分析。經常有人問用混合基組時怎么選擇校正因子,實際上此時沒有完美的做法,一種做法是索性根本不考慮校正因子,還一種也有道理的做法是使用占原子數最多的基組的頻率校正因子。而如果你只對某個振動頻率感興趣,而這個振動模式涉及的原子基本上用的都是基組A,那可以使用A基組對應的頻率校正因子。
要提醒的是,振動分析所用的級別和優化幾何結構的級別必須嚴格一致,而得到的熱力學校正量是可以加到其它更高水平方法算的電子能量上的。比如說,要算一個體系某溫度下的氣相下的自由能,那么可以用B3LYP/6-31G*去優化體系然后在考慮校正因子的情況下做振動分析,得到ZPE、ΔH和熵,再用很好的方法諸如PWPB95-D3/def2-QZVP去獲得這個結構下的電子能量E。最后體系的自由能G=E+ZPE+ΔH-T*S。
實際上,自行擬合頻率校正因子不是難事,見《自行擬合基頻校正因子與ZPE校正因子的簡單方法》(http://www.shanxitv.org/391),用自行擬合的顯然比湊合借用一個更有信服力。
注意校正因子并非總能很好解決振動問題,盡管通常可以。用校正因子后只是整體統計結果變好了,并不代表對每種體系都能很好處理,比如很多級別下對于氟分子的振動頻率算得都不咋樣。并且對于非諧振效應很強、很低頻率的情況,用了校正因子后結果肯定也還是不好,只能直接去做非諧振計算才行。更多相關討論見《基頻頻率校正因子實際效果測試》(http://www.shanxitv.org/390)。
4 在Gaussian中使用頻率校正因子
2020-May-17注:下文介紹的freqchk程序已經沒有任何存在意義了!因為使用《使用Shermo結合量子化學程序計算分子的各種熱力學數據示例》(http://www.shanxitv.org/552)文中介紹的Shermo程序計算熱力學量省事方便得多,且功能強大得多,還可以同時指定不同的校正因子,方便至極,再也不用像下文一樣折騰了!
下面簡單說一下在Gaussian中使用校正因子的方式。在freq任務過程中使用scale=X就代表校正因子為X,例如
#P HF/6-31G* freq scale=0.9135 temperature=323
注意輸出信息中只是在計算熱力學校正量的時候對所用頻率乘上了校正因子,而輸出的頻率還是原始的頻率!所以要想得到校正后的頻率,得自己手動把頻率乘上校正因子。Gaussian顯然是不會自動用校正因子的,即默認scale=1.0。
如果已經做過一次頻率計算并且保存了chk文件,想查看另一個校正因子下的熱力學數據,最好的方法是用Gaussian自帶的freqchk工具,它會從chk文件中讀取儲存的Hessian矩陣并作分析。啟動它之后,輸入chk文件路徑,然后依次輸入
N //不輸出Hyperchem文件
323 //溫度。如果輸入0則為默認的298.15K
1 //壓力。如果輸入0則為默認的1atm
0.8905 //校正因子
[回車]
[回車]
然后立刻就看到了原始頻率以及等同于scale=0.8905 temperature=323時freq任務給出的熱力學數據,和Gaussian自身輸出的一樣。
利用freqchk,就可以反復利用不同校正因子得到較準確的ZPE、ΔH和熵,并獲得自由能G=H-T*S,過程在下面說一下。以下是freq任務會輸出的四個量,為了省事用字母表示,使用不同校正因子時它們都會有變化。
Zero-point correction= A
Thermal correction to Energy= B
Thermal correction to Enthalpy= C 注意這一項不是ΔH,而是ZPE+ΔH
Thermal correction to Gibbs Free Energy= D
首先,在ZPE校正因子下得到ZPE(即A)。然后在熵校正因子下得到-T*S(即D-C)。之后在ΔH校正因子下得到ΔH(即C-A)。最終將這三個量加在一起,再加上高精度方法算的電子能量E,就得到了自由能G=E+ZPE+ΔH-T*S。
在作IR、Raman、VCD這些振動光譜的譜圖時,要先把頻率乘上校正因子后再做圖。手動一一去改Gaussian輸出文件里的頻率值顯然太費事了。Multiwfn (http://www.shanxitv.org/multiwfn)的繪制光譜圖的功能可以方便地設定校正因子。啟動Multiwfn后,載入Gaussian的freq任務的輸出文件,選11進入光譜繪制功能,選IR或Raman,之后選0就可以顯示光譜圖。如果選14 Multiply the vibrational frequencies by a factor,就可以令不同編號范圍的頻率乘以指定的校正因子,重新作圖后就能看到效果了。用戶因此不僅可以對所有頻率乘上統一的校正因子,還可以根據一些文章的思想,對于不同范圍的頻率乘上不同的校正因子以得到更準確結果。關于Multiwfn繪制光譜的更多信息看《使用Multiwfn繪制紅外、拉曼、UV-Vis、ECD、VCD和ROA光譜圖》(http://www.shanxitv.org/224)。
最后順帶一提的是,除了對頻率進行校正,還有人提出了直接對力常數進行校正的方法,不過這樣的方法涉及到太多參數,難搞,所以沒流行起來。對于NMR、UV/Vis等其它譜的計算也有人提出過不同級別下的各種校正方法,不過都遠沒有像頻率校正因子這樣被廣泛使用。