精確計算核磁共振屏蔽常數的擴展焦點分析方法

王康麗,孫　萌,吳安安

(廈門大學化學化工學院,福建 廈門 361005)

通過核磁共振波譜研究可以獲得分子幾何構型、分子中原子的成鍵情況以及相互作用等信息[1]，因此,核磁共振技術是測定分子幾何結構最重要的實驗手段之一[2-4],并且已經被廣泛應用到化學、物理以及醫(yī)學等領域．化學位移因其可以反映原子核周圍的電子環(huán)境而成為核磁共振波譜的重要參數之一．然而,在譜圖的解析過程中仍會出現很多錯誤的結果[5]．屏蔽常數的理論計算不僅可以提高譜圖解析的準確度,而且在復雜分子的結構解析與預測中發(fā)揮著重要的作用,因而受到了越來越多化學家的關注[6-8]．

目前,有很多量子化學計算方法應用到屏蔽常數的理論計算中．采用未考慮電子相關效應的HF(Hartree-Fock)方法雖然能夠獲得部分合理的屏蔽常數,但是越來越多的情況表明電子相關效應對于精確計算屏蔽常數是至關重要的[9-12]．因此,大量考慮電子相關效應的量子化學方法發(fā)展起來,并應用到核磁共振的理論計算中[10,13-15]．其中,二階MP2(M?ller-Plesset)方法雖然在HF方法的基礎上考慮了二階電子相關效應[10,14],但是仍無法很好地描述相關效應很強的體系;密度泛函理論(density functional theory,DFT)方法既考慮了電子相關效應,又能計算包含上百甚至上千個原子的大分子體系,但是現有的DFT方法在計算屏蔽常數時會表現出過強的去屏蔽效應[15];CCSD(T)方法因其能精確計算屏蔽常數而被認為是黃金方法,但是在大基組下,CCSD(T)只能計算不超過4個重原子的體系．因此,為了減少計算量和提高計算精度,人們提出了各種用于計算屏蔽常數的復合方法,如:Sun等[16]提出的焦點分析(focal point analysis for magnetic,FPA-M)方法,其計算精度可以達到CCSD(T)完全基組(complete basis set,CBS)的極限;Kupka等[17]提出的在CCSD(T)/aug-cc-pVTZ-J的基礎上,加上DFT/CBS與DFT/aug-cc-pVTZ-J的差值作為校正項來近似求得CCSD(T)/CBS的值．

雖然對于相對較小的分子體系(小于6個重原子的體系)采用FPA-M方法計算的屏蔽常數能達到CCSD(T)/CBS的外推值,但是應用FPA-M方法時需要計算MP2/cc-pV5Z和CCSD(T)/cc-pVTZ的屏蔽常數,其計算量對于大分子體系,特別是7個重原子以上的體系來說仍然不可承受．為解決這個問題,在FPA-M方法的基礎上,本研究提出了一種計算量更小的方法,并對42個小分子進行了計算．

1　研究方法及研究體系

1.1　研究方法

本研究中屏蔽常數的計算是在CCSD(T)和HF級別下完成的．為了解決規(guī)范原點問題,所有計算都采用GIAO(gauge including atomic orbital)方法[13,18]．基組則采用Dunning開發(fā)的cc-pVnZ基組 (correlation consistent-polarized valence basis set ofnple ζ,n=D、T、Q、5,簡寫為VnZ)[19]、cc-pCVnZ 基組(correlation consistent-polarized core-valence basis set ofnpleζ,n=D、T、Q、5,簡寫為CVnZ)[20]、aug-cc-pVnZ基組(augmented cc-pVnz，n=D、T、Q、5,簡寫為aVnZ)[21];Jensen開發(fā)的pc-n(polarization consistent basis sets,n=2～4)[22]、pcJ-n(n=2～4)[23]和pcS-n(n=2～4)[24]基組．在本研究中,采用指數形式的兩點[25]和指數函數/高斯函數混合的三點[26]外推至CBS,公式表達如下:

表1　研究的4個分子集

Tab.1　Four sets of molecules under investigation

分子集分子結構基準1C2H2,C2H4,CH2N2,CH2NN,CH2NN,CH4,CO2,CO2,CO,CO,F2,HF,H2CO,H2CO,H2O,HCN,HCN,N2,NNO,NNO,NNO,NH3實驗構型CCSD(T)/CBS-345-VnZ2CH2CH+,CH2CH+,CF2,CF2,cis-N2F2,cis-N2F2,trans-N2F2,trans-N2F2,N2O-cyc,N2O-cyc優(yōu)化構型CCSD(T)/CBS-345-VnZ3CH3COCH3,CH3COCH3,CH3COCH3,C6H6,C2H6,CH2CCH2,CH2CCH2,CH3CH2CH3,CH3CH2CH3,CF4,CF4,CH3CN,CH3CN,CH3CN,CH3F,CH3F,CH3NH2,CH3NH2,CH3OH,CH3OH,CHF3,CHF3實驗構型FPA-M/CBS-345-VnZ4CO2-cyc,CO2-cyc,NO2-,NO2-,C6H5NN+,C6H5NN+,C6H5NN+,優(yōu)化構型FPA-M/CBS-345-VnZ

注：下劃線表示本文中研究的核;優(yōu)化構型是在MP2/VQZ級別下優(yōu)化,實驗構型及優(yōu)化構型的坐標見文獻[17]．

σe(n)=σe(CBS)+ae1-X，

(1)

σe(n)=σe(CBS)+ae1-X+be-(1-X)2，

(2)

其中,σe(n)為計算所得的屏蔽常數,X為上述基組的指數值,σe(CBS)為CBS下屏蔽常數的外推值．σe(CBS)、a、b這3個參數由非線性最小二乘法擬合得到．為了更好地區(qū)別上述兩種外推方法得到的屏蔽常數,用CBS-34、CBS-45、CBS-345分別表示將X=3～4帶入式(1),將X=4～5帶入式(1),將X=3～5帶入式(2)得到的值．對于VnZ、CVnZ和aVnZ基組,X=3對應TZ基組,X=4對應QZ基組,以此類推;對于pc-n、pcJ-n和pcS-n基組,X=3對應pc-2、pcJ-2、pcS-2,X=4對應pc-3、pcJ-3、pcS-3,以此類推．

FPA-M方法是在MP2/CBS的極限值上,加上CCSD(T)/VTZ與MP2/VTZ的差值作為高階相關校正項,其計算公式如下:

σe(CCSD(T)/CBS)≈σe(MP2/CBS)+

Δσe(CCSD(T)),

(3)

Δσe(CCSD(T))=σe(CCSD(T)/VTZ)-

σe(MP2/VTZ).

(4)

本研究在FPA-M方法的基礎上,提出FPA-M-HF方法,即外推HF屏蔽常數至CBS,并加上較小基組下CCSD(T)與HF計算屏蔽常數的差值作為高階相關校正項,以獲得近似CCSD(T)/CBS的屏蔽常數,其計算公式如下:

σe(CCSD(T)/CBS)≈σe(HF/CBS)+

Δσe(CCSD(T)),

(5)

Δσe(CCSD(T))=σe(CCSD(T)/small)-

σe(HF/small),

(6)

其中,small表示小基組,σe(HF/CBS)表示在HF級別下用外推公式得到的屏蔽常數．在本研究中,VTZ、CVTZ、aVTZ、pc-2、pcJ-2 和 pcS-2被選作小基組．

本研究中CCSD(T)計算采用CFOUR程序完成,HF計算采用GAUSSIAN 09程序完成．

1.2　研究體系

本文中研究的分子分為4個分子集,見表1．分子集1和3采用的是實驗構型,分子集2和4采用的是優(yōu)化構型．在給出屏蔽常數的誤差時,分子集1和2以CCSD(T)/CBS-345-VnZ方法計算的屏蔽常數為基準；而分子集3和4中分子包含的原子數比較多,因此很難計算出CCSD(T)/V5Z級別的屏蔽常數．從文獻[16]可知FPA-M方法的計算精度達到CCSD(T)/CBS,所以對于分子集3和4,以FPA-M/CBS-345-VnZ方法計算的屏蔽常數為基準．

2　結果與討論

2.1　分子集1和2的計算結果

采用FPA-M-HF/CBS-34方法下的不同基組計算分子集1和2的屏蔽常數,見附錄(jxmu.xmu.edu.cn/upload/html20180202.html)表S1．從表2可看出，在FPA-M-HF/CBS-34方法下相較于其他基組,采用pcJ-n和CVnZ基組計算分子集1的屏蔽常數能夠獲得很好的計算結果,其平均絕對偏差分別為0.6×10-6和1.1×10-6,最大誤差分別為2.5×10-6和3.3×10-6,尤其是pcJ-n的平均絕對偏差只有CCSD(T)/V5Z的3/8(CCSD(T)/VnZ的計算結果見附錄表S2)．采用FPA-M-HF/CBS-34方法下VnZ和aVnZ基組計算屏蔽常數的最大誤差都來自F2中的F原子,分別為9.0×10-6和9.8×10-6,這是因為CCSD(T)和HF方法計算的屏蔽常數之間的差值受基組的影響較大,即差值的系統性比較差．以VnZ基組為例,當n=T、Q、5時,采用CCSD(T)方法計算F2中F原子的屏蔽常數分別為-165.9×10-6,-178.1×10-6和-185.4×10-6,而采用HF方法計算得到的屏蔽常數分別為-151.2×10-6,-160.6×10-6和-167.0×10-6,兩者之間的差值分別為-14.7×10-6,-17.5×10-6和-18.4×10-6,即差值隨基組變化比較明顯,這表明用[σe(CCSD(T)/VTZ)-σe(HF/VTZ)]作為高階校正項的不準確性．當用更高階差值[σe(CCSD(T)/VQZ)-σe(HF/VQZ)]作為校正項時,可將F2的誤差縮小至6.1×10-6,但是其計算量會大大地增加．而FPA-M-HF/CBS-34-CVnZ、pc-n和pcS-n的最大誤差分別來自H2CO，分別為C原子的3.3×10-6、O原子的14.4×10-6和O原子的12.4×10-6,這主要是因為含有雜原子多重鍵的體系需要更高階的相關效應來描述．進一步研究表明,在FPA-M-HF/CBS-34方法下,pcJ-n和CVnZ基組的平均絕對偏差和最大誤差都小于CCSD(T)/V5Z的相應值(平均絕對偏差為1.6×10-6,最大誤差為4.3×10-6),說明這兩種方法的精度比CCSD(T)/V5Z高;盡管VnZ、aVnZ和pcS-n基組計算的屏蔽常數的精度(平均絕對偏差分別為3.2×10-6,2.9×10-6和3.8×10-6)沒有達到CCSD(T)/V5Z的精度,但是都遠遠超過了CCSD(T)/VQZ(平均絕對偏差為4.1×10-6)．

從表2還可以看出對于分子集2,FPA-M-HF/CBS-34-pcJ-n和CVnZ的平均絕對偏差分別為1.6×10-6和1.5×10-6,均比CCSD(T)/V5Z的(2.2×10-6)小;而VnZ、aVnZ和pcS-n基組的平均絕對偏差則比CCSD(T)/VQZ的(5.9×10-6)小．

表2　采用FPA-M-HF和CCSD(T)方法計算分子集1和2的屏蔽常數

注：MAD表示平均絕對偏差,MAX表示最大誤差．

采用FPA-M-HF/CBS-45方法計算分子集1和2的屏蔽常數,見附錄表S3．從表2可以看出對于FPA-M-HF/CBS-45計算的分子集1,pcJ-n和CVnZ基組仍然能夠給出很精確的計算結果,其平均絕對偏差分別為0.8×10-6和1.1×10-6,最大誤差分別為2.5×10-6和3.2×10-6；VnZ和aVnZ基組的最大誤差來自F2中的F原子,分別為4.3×10-6和7.2×10-6,而CVnZ、pc-n和pcS-n的最大誤差則來源于H2CO,分別為C原子的3.2×10-6、O原子的16.4×10-6和O原子的7.7×10-6,其原因已經在前文論述過,這里不再贅述．此外,pcJ-n基組的最大誤差來自CH2N2中的Nt(頂端N原子),這是因為CH2N2屬于雜原子多重鍵體系,對于這種體系需要高階相關效應來描述．同時，在FPA-M-HF/CBS-45方法中,pcJ-n和CVnZ基組的精度遠高于CCSD(T)/V5Z的精度,尤其是pcJ-n的平均絕對偏差只有CCSD(T)/V5Z的1/2;VnZ和aVnZ的精度(平均絕對偏差分別為1.6×10-6和1.7×10-6)與CCSD(T)/V5Z相當;盡管pcS-n基組計算的屏蔽常數的精度(平均絕對偏差為3.0×10-6)低于CCSD(T)/V5Z,但卻比CCSD(T)/VQZ的精度要高．對于分子集2,FPA-M-HF/CBS-45-pcJ-n、VnZ、CVnZ的平均絕對偏差只有1.7×10-6,1.7×10-6和1.8×10-6,比CCSD(T)/V5Z的(平均絕對偏差為2.2×10-6)要?。梢园l(fā)現在FPA-M-HF方法中,VnZ、aVnZ和pcS-n基組下CBS-45的精度要比CBS-34的高許多,CVnZ和pcJ-n基組下的差異不大,而pc-n基組下則是CBS-45的精度低于CBS-34．進一步分析發(fā)現,HF級別下pc-n基組的CBS與其他5種基組的相差比較大,尤其是CBS-345和CBS-45,所以計算的分子集1和2的平均絕對偏差也比較大．

1.CCSD(T)/VTZ；2.FPA-M-HF/CBS-45-pc-n；3.FPA-M-HF/CBS-34-pc-n;4.CCSD(T)/VQZ；5.FPA-M-HF/CBS-34-pcS-n；6.FPA-M-HF/CBS-45-pcS-n；7.FPA-M-HF/CBS-34-aVnZ；8.FPA-M-HF/CBS-34-VnZ；9.FPA-M-HF/CBS-45-aVnZ；10.FPA-M/CBS-34-VnZ；11.CCSD(T)/V5Z；12.FPA-M-HF/CBS-45-VnZ；13.FPA-M-HF/CBS-45-CVnZ；14.FPA-M-HF/CBS-34-CVnZ；15.FPA-M-HF/CBS-45-pcJ-n；16.FPA-M-HF/CBS-34-pcJ-n.圖1　分子集1和2中4種元素的平均絕對偏差Fig.1 Mean absolute deviations (MADs) for four elements in sets 1 and 2

由表2可以看出采用FPA-M-HF/CBS-345和CBS-45(附錄表S4)計算的屏蔽常數幾乎相同,這是因為隨著基組的增大,基組的收斂趨勢越來越平滑、明確．因為二者的結果一樣,這里不再贅述．

為了更好地分析FPA-M-HF方法對于不同元素的計算誤差,給出分子集1和2中的4種元素(13C、15N、17O、19F)的平均絕對偏差,并按照4種元素總的平均絕對偏差從大到小排序(圖1)．因為對于FPA-M-HF來說,采用CBS-345和CBS-45得到的屏蔽常數幾乎相同,所以這里僅列出了CBS-34和CBS-45的計算結果．從圖1可以看出,除了FPA-M-HF/CBS-45-pc-n和34-pc-n方法外,其余方法的平均絕對偏差都小于CCSD(T)/VQZ．對于FPA-M-HF/pcS-n、aVnZ、34-VnZ方法,總的平均絕對偏差雖然低于CCSD(T)/VQZ,但對于15N和19F來說,其平均絕對偏差比較大．考慮到計算量和4種元素的平均絕對偏差,FPA-M-HF/CBS-34-CVnZ方法可以看作是最好的計算方法,其計算量要比FPA-M-HF/CBS-34-pcJ-n和CCSD(T)/VQZ的計算量要小,且對于4種元素計算的平均絕對偏差都在1.5×10-6以內．

為了進一步比較FPA-M-HF和FPA-M方法,給出了FPA-M-HF/CBS-34-CVnZ與FPA-M/CBS-345-VnZ計算甲烷中13C屏蔽常數的時間(分別為11 min 16 s和21 min 18 s)以及誤差(分別為0.2×10-6和0)．雖然FPA-M-HF/CBS-34-CVnZ計算的誤差達到了0.2×10-6,但是卻將計算時間縮減至約FPA-M/CBS-345-VnZ的一半,這說明FPA-M-HF/CBS-34-CVnZ在保證精度的前提下,極大地減少了計算量．

2.2　分子集3和4的計算結果

表3列出了以FPA-M/CBS-345-VnZ計算的屏蔽常數為基準,用FPA-M-HF/CBS-34-CVnZ方法計算所得分子集3和4的屏蔽常數，其平均絕對偏差只有1.2×10-6,這說明對于絕大多數體系,該方法能夠得到很高的計算精度．最大誤差來源于丙酮中的O原子(7.3×10-6),這是因為丙酮是含有多重鍵的體系,需要更高階相關校正項．在用FPA-M方法計算C6H5N2+和C7H9+的屏蔽常數時,用的是MP2關于VDZ、VTZ和VQZ的外推值,而FPA-M-HF/CBS-34-CVnZ計算這二者的誤差都不超過1.0×10-6．對于像C6H5N2+和C7H9+這樣的體系,無法用CCSD(T)/V5Z來計算．但是FPA-M-HF卻能在較小計算量的基礎上得到高精度計算結果,這也是FPA-M-HF方法的優(yōu)勢所在．

表3　HF/CVnZ、HF/CBS-34、CCSD(T)/VTZ、FPA-M-HF/CBS-34-CVnZ計算分子集3和4的屏蔽常數

注：CO表示C原子連在O上;CN表示C原子連在N上;Cc表示中心C原子;Oc表示中心O原子;C1表示C原子與N相連;C2、C3、C4分別表示鄰、間、對位的C原子;N1表示N原子連在芳香環(huán)上且?guī)д姾?N2表示N原子僅與N相連;C7H9+中的C原子標號見表1．

3　結　論

本文中系統地研究了在FPA-M方法基礎上發(fā)展起來的FPA-M-HF方法,并對比了該方法下VnZ、 CVnZ、 aVnZ、pc-n、 pcJ-n和pcS-n6種基組的計算結果．FPA-M-HF方法能夠精確計算分子的屏蔽常數,尤其是用CVnZ和pcJ-n基組,其精度要高于CCSD(T)/V5Z,而計算量要比相應基組下的FPA-M方法小．相比于CCSD(T)/V5Z只能計算不超過4個重原子的體系,FPA-M-HF方法能計算包含7個甚至更多個重原子的體系．綜合考慮精確度和計算量,FPA-M-HF/CBS-34-CVnZ方法更好,其計算量要比CCSD(T)/VQZ小很多．這為進一步研究更大分子體系提供了指導．

參考文獻：

[1] 萬堅.核磁共振光譜參數的理論研究[D].武漢:華中師范大學,2001:1-2.

[2]WUTHRICH K.NMR of proteins and nucleic acids[M].New York:Wiley-Interscience,1986.

[3]CORNILESCU G,DELAGLIO F,BAX A.Protein backbone angle restraints from searching a database for chemi-cal shift and sequence homology[J].Journal of Biomole-cular NMR,1999,13(3):289-302.

[4]SANDERS J K M,HUNTER B K.Modern NMR spectroscopy:a guide for chemists[M].Oxford:Oxford University Press,1987.

[5]NICOLAOU K C,SNYDER S A.Chasing molecules that were never there:misassigned natural products and the role of chemical synthesis in modern structure elucidation[J].Angewandte Chemie International Edition,2005,44(7):1012-1044.

[7]BIFULCO G,DAMBRUOSO P,GOMEZ-PALOMA L,et al.Determination of relative configuration in organic compounds by NMR spectroscopy and computational methods[J].Chemical Reviews,2007,107(9):3744-3779.

[8]CAVALLI A,SALVATELLA X,DOBSON C M,et al.Protein structure determination from NMR chemical shifts[J].Proceedings of the National Academy of Sciences of the United States of America,2007,104(23):9615-9620.

[9]PRICE D R,STANTON J F.Computational study of [10] annulene NMR spectra[J].Organic Letters,2002,4(17):2809-2811.

[10]GAUSS J.Calculation of NMR chemical shifts at second-order many-body perturbation theory using gauge-including atomic orbitals[J].Chemical Physics Letters,1992,191(6):614-620.

[11]GAUSS J.Effects of electron correlation in the calculation of nuclear magnetic resonance chemical shifts[J].The Journal of Chemical Physics,1993,99(5):3629-3643.

[12]BüHL M,GAUSS J,HOFMANN M,et al.Decisive electron correlation effects on computed11B and13C NMR chemical shifts.Application of the GIAO-MP2 method to boranes and carbaboranes[J].Journal of the American Chemical Society,1993,115(26):12385-12390.

[14]AUER A A.Quantitative prediction of gas-phase17O nuclear magnetic shielding constants[J].The Journal of Chemical Physics,2009,131(2):024116.

[15]孫萌.核磁共振屏蔽波譜的理論計算[D].廈門:廈門大學,2013:19-51.

[16]SUN M,ZHANG I Y,WU A A,et al.Accurate prediction of nuclear magnetic resonance shielding constants:towards the accuracy of CCSD (T) complete basis set limit[J].The Journal of Chemical Physics,2013,138(12):124113.

[18]DITCHFIELD R.Molecular orbital theory of magnetic shielding and magnetic susceptibility[J].The Journal of Chemical Physics,1972,56(11):5688-5691.

[19]DUNNING T H,JR.Gaussian basis sets for use in correlated molecular calculations.Ⅰ.The atoms boron through neon and hydrogen[J].The Journal of Chemical Physics,1989,90(2):1007-1023.

[20]WOON D E,DUNNING T H,JR.Gaussian basis sets for use in correlated molecular calculations.Ⅴ.Core-valence basis sets for boron through neon[J].The Journal of Chemical Physics,1995,103(11):4572-4585.

[21]KENDALL R A,DUNNING T H,JR,HARRISON R J.Electron affinities of the first-row atoms revisited.Systematic basis sets and wave functions[J].The Journal of Chemical Physics,1992,96(9):6796-6806.

[22]JENSEN F.Polarization consistent basis sets:principles[J].The Journal of Chemical Physics,2001,115(20):9113-9125.

[23]JENSEN F.The basis set convergence of spin-spin coup-ling constants calculated by density functional methods[J].Journal of Chemical Theory and Computation,2006,2(5):1360-1369.

[24]JENSEN F.Basis set convergence of nuclear magnetic shielding constants calculated by density functional methods[J].Journal of Chemical Theory and Computation,2008,4(5):719-727.

[25]FELLER D,SORDO J A.A CCSDT study of the effects of higher order correlation on spectroscopic constants.Ⅰ.First row diatomic hydrides[J].The Journal of Chemical Physics,2000,112(13):5604-5610.

[26]PETERSON K A,WOON D E,DUNNING T H,JR.Benchmark calculations with correlated molecular wave functions.Ⅳ.The classical barrier height of the H+H2→H2+H reaction[J].The Journal of Chemical Physics,1994,100(10):7410-7415.