密度泛函计算中基组的收敛现象:强成键体系CO分子

黄莺，雷志丹，张晓青，李龙

（湖南中医药大学药学院，中国长沙 410208）

密度泛函计算中基组的收敛现象:强成键体系CO分子

黄莺*，雷志丹，张晓青，李龙

（湖南中医药大学药学院，中国长沙 410208）

探讨了用密度泛函理论中7个常见交换－相关密度泛函及8种常见Dunning调和相关基组处理CO分子的三种性质（平衡键长、谐振基频和离解能）时表现出的基组收敛现象.用Dunning的单指数插补法，求得了上述性质的全基组极限值.用扩展的调和相关基组，发现键长的系统误差在0.002 9（LDA）到0.007 3 Å（BP86）之间，基频误差为17（LDA）至56（B3BP86）cm－1，离解能误差为12.6（BLYP）至186.6（LDA）kJ/mol.这些结果对系统地认识波函数的基组收敛现象将提供有益启示.

密度泛函理论;基组收敛;调和相关基组

近年来，以密度泛函理论（DFT）为基础的计算方法的运用得到了广泛的普及和承认［1－6］.DFT的理论基础源于Hohenberg和Kohn的两个著名定理，以电子密度作为多电子体系量化处理的基本参量，使之包含原子与分子基态性质的所有信息［7］.

几乎在所有情况下，DFT的具体计算都涉及到Kohn－Sham轨道方程的求解［8］

从方程（1）中可以看出，Kohn－Sham方程引入了轨道概念.这样，DFT方法与以波函数为基础的常规从头计算的方法存在同样的缺陷，即用非完全基组对波函数的展开.对全基组（CBS）极限的研究将有利于了解DFT计算中给定泛函的系统偏差.Dunning等人［9－10］曾发现分子性质对基组选择呈良好的单指数函数收敛行为.据此，人们便可估计出全基组极限及相应的计算性质［11－12］的系统偏差.

Geerlings等人［13］曾就电荷分布，化学活性指数，及相关性质的全基组收敛行为作了探讨.本文运用7种常见交换相关泛函组合，以CO分子为例，对强成键体系CO分子的全基组收敛性质及不同泛函的表现进行了系统研究.

1 方法

DFT方法中一给定泛函的特定精确度可用其所求得一系列分子性质的全基组极限值来评估.本文所用性质为分子平衡键长re，解离能D0，和谐振基频ωe.要求得全基组极限值，必须知道基组扩充的有效收敛规律.Dunning的扩展的调和相关价层双ZETA（（aug）－ccp－VDZ），三ZETA（（aug）－ccp－VTZ），四ZETA（（aug）－ccp－VQZ），和五ZETA（（aug）－ccp－V5Z）基组即属此类基组［14］.对此类基组，人们发现［15］某种给定性质的计算值Q（X）能用如下单指数函数很好地表达:

其中，X为所用基组的ZETA数，Q（∞）为全基组极限值［9－10］.一般说来，理论与试验的误差，ΔQapp，与基组大小的关系有两种不同类型.第一种随基组大小的增加而逐渐减小，故计算值随基组的不断增大而越来越靠近实验值.第二种类型则有一极小点.在该极小点，计算值与实验值非常接近，但仍与Q（∞）有很大差别.此类型常会被误认为所用泛函的质量很好.

交换－相关泛函可表达为

其中，EX［ρ］与EC［ρ］分别是交换和相关能量密度泛函.如果它们仅是密度（（r）本身的函数，则称它们是局域的，否则是非局域的.本文将采用7种常用的交换－相关泛函.对交换泛函，我们采用Slater（S）［16］，梯度修正的Becke 88（B）［17］，及Hartee－Fock交换式［18］，而对相关能泛函则采用Vosko－Wilk－Nusair（VWN）［19］，梯度修正的Perdew 86（P86）［20］，Perdew－Wang91（PW91）［21］及Lee－Yang－Parr（LYP）［22］公式.二者组合可得如下7种常见方法:SVWN，BLYP，B3LYP，BP86，B3P86，BPW91和B3PW91.

所有计算采用Gaussian′03 Window版［23］软件在微机上完成.

2 结果与讨论

2.1 平衡键长re

表1列出了用7种泛函组合与8种Dunning基组计算得到的平衡键长.实验值为1.128 3 Å.由表中可见，计算值很好的收敛到CBS极值.然而，各个近似泛函的收敛行为却不尽相同.所有梯度修正的泛函均表现为第一类收敛行为，而SVWN和以Hartee－Fock交换能为基础的混合泛函方法，B3LYP，B3P86，和B3PW91则表现为第二类收敛行为.这表明用后面4种泛函，人们可以用较小基组得到较精确的平衡键长re值，但这并不说明这些计算方法较其他更好.一个例子是SVWN/cc－pVTZ，其re值为1.128 1 Å，与实验值已非常接近，但该计算水平总的说来非常粗糙.

表1 用不同泛函和不同基组近似所得CO分子的平衡键长及各计算方法的全基组（CBS）极值和系统误差单位:Å

对给定基组下某性质计算值与全基组极限值的比较可用下式表达（ΔQbs/Q（∞））×100，其中ΔQbs为计算值与极限值之差.从表中可得，对（aug）－cc－pVDZ基组该值为±1%，（aug）－cc－pVTZ为±0.2%，（aug）－ccpVQZ为±0.02%，（aug）－cc－pV5Z为±0.005%.ZETA数越大，偏差率越小.另外，从表中可见，基组的扩展与否对平衡键长计算值的影响非常小，系统误差的差别小于1 mÅ.与CBS极值相比，SVWN方法下的小扩展调和相关基组具有最大的偏差，但随ZETA数增加而快速收敛至CBS极值.

从系统误差来看，不同泛函组合的差别并不大，以SVWN的结果最小，低估2.9 mÅ.混合泛函平均低估5.0 mÅ，而梯度修正泛函平均高估7.0 mÅ.

2.2 谐振基频ωeC

在表2中，我们列出了用上述7种方法、8种基组计算所得CO分子的谐振基频，实验值为217 0 cm－1［24］.对于扩展的调和相关基组，计算值很好地收敛到CBS极值，其中梯度修正泛函表现为第一类收敛行为，而SVWN和混合泛函表现为第二类收敛行为.对非扩展的调和相关基组，由于从cc－pVQZ到cc－pV5Z有的没有增加，有的反而降低，故用（4）式无法得到其CBS极值.

表2 用不同泛函和不同基组近似所得CO分子的谐振基频ωe及各计算方法的全基组（CBS）极值和系统误差单位:cm－1

如果考察ωe计算值的精度随基组大小变化的关系，我们发现，跟平衡键长类似，与CBS极值相比较，SVWN泛函在较小基组时给出较大的偏差.随基组ZETA数的增加，其值迅速收敛至CBS极值.

此外，同一基组不同泛函所得结果的差别并不大.对aug－cc－pVDZ基组，变化率为1.2%，aug－cc－pVTZ为0.2%，aug－cc－pVQZ为0.04%，aug－cc－pV5Z为0.008%.从系统误差来看，以SVWN的结果最小，高估17 cm－1，混合泛函平均高估50 cm－1，而梯度修正泛函平均则平均低估50 cm－1.

2.3 离解能D0

表3报道了对CO分子离解能的计算值及其比较，其中实验值为1 071.9 kJ/mol［24］.我们观察到了与谐振基频一样的收敛行为，即对扩展的调和相关基组，可以看到很好的收敛现象，而对另外的基组，由于从ccpVQZ到cc－pV5Z，计算值不增反降，用（4）式无法得到其CBS极值.对扩展的调和相关基组而言，混合泛函B3LYP及B3PW91表现为第一类收敛行为，其它则为第二类.各基组的偏差百分率分别为2.7%（aug－ccpVDZ），0.5%（aug－cc－pVTZ），0.06%（aug－cc－pVQZ），和0.01%（aug－cc－pV5Z）.与平衡键长及谐振基频一致，ZETA数越大，偏差值越小.从系统误差而言，SVWN泛函的最大，高估186.6 kJ/mol，而梯度修正及混合泛函均表现出良好的精度.

表3 用不同泛函和不同基组近似所得CO分子的离解能D0及各计算方法的全基组（CBS）极值和系统误差单位:kJ/mol

2.4 讨论

在现代密度泛函理论中，分子结构与性质的计算精度由两个彼此独立的因素决定，即交换－相关能的近似密度泛函形式及基组的相对大小（也就是由Kohn－Sham轨道组成的近似波函数的精确性）.由于前者相对于后者更难精确地近似，DFT过去几十年的发展均以前者为研究核心，代表性工作有LDA，GGA，meta－GGA，hybrid－GGA，double－hybrid GGA等.而对基组的作用特别是基组在计算结构和性质时表现出的收敛现象及其规律，仍然不见系统研究，而且在文献中也没有引起足够的重视.本文以强成键体系CO分子为例，对3种常见物理化学性质，即平衡键长、谐振基频和离解能，在计算时表现出的全基组收敛现象进行了系统研究，同时采用不同的泛函近似形式探讨了不同的全基组收敛模式.发现不同的性质不同的泛函表现出不同的全基组收敛方式，反映了全基组收敛的复杂性.更有意思的是，从全基组收敛的角度来看，波函数计算的精确性与近似泛函的精确性并没有相关性.这些结果表明，上述决定分子结构和性质计算精确度的两个因素是彼此独立的，没有相关性.

3 结论

本文用密度泛函理论中7种常见泛函组合探讨了用8种扩展和非扩展的调和相关基组来计算CO分子3种物理化学性质（即平衡键长，谐振基频和离解能）时表现出的全基组收敛现象.并利用Dunning等人的单指数插补函数确定了这些性质的CBS极值.结果发现，不同的性质可能表现出截然不同的收敛行为.对平衡键长与谐振基频而言，定域密度近似（LDA）和混合泛函方法表现为第二类收敛行为，而梯度修正则为第一类.对离解能，B3LYP和B3PW91两个混合泛函方法表现出第一类，而其它则为第二类.令人吃惊的是，LDA方法给出了该分子平衡键长与谐振基频的最好结果，而梯度修正（GGA）与混合泛函（hybrid GGA）方法的结果相当.对其它成键体系的研究正在进行中.

［1］PARR R G，YANG W.Density functional theory of atoms and molecules［M］.New York:Oxford University Press，1989.

［2］DREIZLER R M，GROSS E K U.Density functional theory［M］.New York:Springer－Verlag，Berlin Heidelberg，1990.

［3］POLITZER P，SEMINARIO J.Modern density functional theory;a tool for chemistry［M］.Amsterdam:Elsevier，1995.

［4］PARR R G，YANG W.Density－functional theory of the electronic structure of molecules［J］.Ann Rev Phys Chem，1995，46: 701－728.

［5］KOHN W，BECKE A D，PARR R G.Density functional theory of electronic structure［J］.J Phys Chem，1996，100:12974－12980.

［6］LIU S B.Conceptual density functional theory and some recent developments［J］.Acta Phys－Chim Sin，2009，25:590－600.

［7］HOHENBERG P，KOHN W.Inhomogeneous electron gas［J］.Phys Rev B，1964，136:864－871.

［8］KOHN W，SHAM L J.Self－consistent equations including exchange and correlation effects［J］.Phys Rev A，1965，140:1133－1138.

［9］WOON D E，DUNNING JR T H.Benchmark calculations with correlated molecular wave functions.Ⅰ.Multireference configuration interaction calculations for the second row diatomic hydrides［J］.J Chem Phys，1993，99:1914－1929.

［10］PETERSON K A，KENDALL R A，DUNNING JR T H.Benchmark calculations with correlated molecular wave functions.Ⅱ. Configuration interaction calculations on first row diatomic hydrides［J］.J Chem Phys，1993，99:1930－1944.

［11］黄雪，刘文龙，钟爱国，等.双核镉配聚体及其衍生物荧光光谱发光机理的计算化学研究［J］.湖南师范大学自然科学学报，2010，33（4）:69－74.

［12］HUANG Y，ZHONG A G，LIU S B.Prediciting pKa values for singly and multiply substituted benzoic acids with density functional reactivity theory［J］.湖南师范大学自然科学学报，2011，34（1）:52－55.

［13］GEERLINGS P，DE PROFT F，LANGENAEKER W.Density functional theory:a source of chemical concepts and a costeffective methodology for their calculation［J］.Adv Quant Chem，1999，33:303－328.

［14］DUNNING JR T H.Gaussian basis sets for use in correlated molecular calculations.Ⅰ.The atoms boron through neon and hydrogen［J］.J Chem Phys，1989，90:1007－1023.

［15］MARTIN J M L.Energetics of stable molecules and reactive intermediates［M］.NATO ASI Symposium Volume，M E Minas da Piedade and K K Irikura（Ed.），Dordrecht:Kluwer Academic Publishers，Dordrecht，1999.

［16］SLATER J C.A simplification of the hartree－Fock method［J］.Phys Rev，1951，81:385－391.

［17］BECKE A D.Density－functional exchange－energy approximation with correct asymptotic behavior［J］.Phys Rev A，1988，38: 3098－3100.

［18］BECKE A D.Density－functional thermochemistry.Ⅲ.The role of exact exchange［J］.J Chem Phys，1993，98:5648－5652.

［19］VOSKO S H，WILK L，NUSAIR M.Accurate spin－dependent electron liquid correlation energies for local spin density calculations:a critical analysis［J］.Can J Phys，1980，58:1200－1211.

［20］PERDEW J P.Density－functional approximation for the correlation energy of the inhomogeneous electron gas［J］.Phys Rev B，1986，33:8822－8824.

［21］PERDEW J P.Electronic structure of solids［M］.P Ziesche and H Eschrig（Ed.）Berlin:Akademie Verlag，1991.

［22］LEE C，YANG W，PARR R G.Development of the colle－salvetti correlation－energy formula into a functional of the electron density［J］.Phys Rev B，1988，37:785－789.

［23］FRISCH M J，TRUCKS G W，SCHLEGEL H B，et al.Gaussian 03，Revision E01，Gaussian，Inc.，Pittsburgh PA，2003.

［24］HUBER K P，HERZBERG G.Molecular spectra and molecular structureⅣ.Constants of diatomic molecules［M］.Princeton:Van Nostrand，1979.

Basis Set Convergence Phenomenon in Density Functional Theory: CO as an Example of Strongly Bonded Systems

HUANG Ying*，LEI Zhi-dan，ZHANG Xiao-qing，LI Long

In this work we systematically investigate the basis set convergence phenomenon with seven commonly used approximate exchange－correlation energy density functionals and eight Dunning correlated basis sets for three physicochemical properties，equilibrium bond distance，first harmonic frequency，and dissociation energy，of the CO molecule.By employing Dunning′s single－exponent extrapolation approach，we have obtained the complete basis set limit for these properties.With the augmented correlation basis sets employed in this study，we found that the systematic error in the equilibrium bond distance is between 0.002 9（LDA）to 0.007 3 Å（BP86），that in the first harmonic frequency lies between 17（LDA）and 56（B3BP86）cm－1，and that in the dissociation energy falls between 12.6（BLYP）and 186.6 kJ/mol.These results shed lights on our systematic understanding of the basis set convergence phenomenon in the wave function.

density functional theory;basis set convergence;correlated basis set

O641.12+1

A

1000－2537（2011）05－0050－05

2011－05－20

湖南省自然科学基金资助项目（11JJ5065）;湖南省大学生研究性学习和创新性实验计划项目（湘教通［2011］272号）

*通讯作者，E－mail:huanglingying@126.com

（编辑杨春明）