精确计算核磁共振屏蔽常数的扩展焦点分析方法

王康丽,孙　萌,吴安安

(厦门大学化学化工学院,福建 厦门 361005)

通过核磁共振波谱研究可以获得分子几何构型、分子中原子的成键情况以及相互作用等信息[1]，因此,核磁共振技术是测定分子几何结构最重要的实验手段之一[2-4],并且已经被广泛应用到化学、物理以及医学等领域．化学位移因其可以反映原子核周围的电子环境而成为核磁共振波谱的重要参数之一．然而,在谱图的解析过程中仍会出现很多错误的结果[5]．屏蔽常数的理论计算不仅可以提高谱图解析的准确度,而且在复杂分子的结构解析与预测中发挥着重要的作用,因而受到了越来越多化学家的关注[6-8]．

目前,有很多量子化学计算方法应用到屏蔽常数的理论计算中．采用未考虑电子相关效应的HF(Hartree-Fock)方法虽然能够获得部分合理的屏蔽常数,但是越来越多的情况表明电子相关效应对于精确计算屏蔽常数是至关重要的[9-12]．因此,大量考虑电子相关效应的量子化学方法发展起来,并应用到核磁共振的理论计算中[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级别下完成的．为了解决规范原点问题,所有计算都采用GIAO(gauge including atomic orbital)方法[13,18]．基组则采用Dunning开发的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开发的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优化构型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+,优化构型FPA-M/CBS-345-VnZ

注：下划线表示本文中研究的核;优化构型是在MP2/VQZ级别下优化,实验构型及优化构型的坐标见文献[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个参数由非线性最小二乘法拟合得到．为了更好地区别上述两种外推方法得到的屏蔽常数,用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采用的是优化构型．在给出屏蔽常数的误差时,分子集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)要小．可以发现在FPA-M-HF方法中,VnZ、aVnZ和pcS-n基组下CBS-45的精度要比CBS-34的高许多,CVnZ和pcJ-n基组下的差异不大,而pc-n基组下则是CBS-45的精度低于CBS-34．进一步分析发现,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方法的优势所在．

表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原子连在芳香环上且带正电荷;N2表示N原子仅与N相连;C7H9+中的C原子标号见表1．

3　结　论

本文中系统地研究了在FPA-M方法基础上发展起来的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.