统计方法校正O3LYP方法计算的生成热

王秀军 龙 汨

(华南理工大学化学与化工学院,广东省燃料电池技术重点实验室,广州510640)

统计方法校正O3LYP方法计算的生成热

王秀军*龙 汨

(华南理工大学化学与化工学院,广东省燃料电池技术重点实验室,广州510640)

由于引入各种内在近似,密度泛函理论存在固有误差.本文采用O3LYP/6-311+G(3df,2p)//O3LYP/ 6-31G(d)计算了220个中小型有机分子的生成热(ΔfH⊖calc),随后应用神经网络(ANN)和多元线性回归(MLR)方法对ΔfH⊖calc进行校正.采用计算得到的生成热、零点能、分子中原子总数、氢原子个数、双中心成键电子数、双中心反键电子数、单中心价层孤对电子数、单中心内层电子数作为ANN和MLR的描述符.以180个分子作为训练集构造ANN或MLR模型,并对40个独立测试集分子的ΔfH⊖calc进行了预测.结果表明:经过ANN和MLR校正后,训练集分子生成热的理论计算值和实验值间的均方根偏差(RMSD)从24.7 kJ·mol-1分别降低到11.8、13.0 kJ· mol-1;独立测试集分子的RMSD从21.3 kJ·mol-1分别降低到10.4、12.1 kJ·mol-1.因此ANN模型的拟合和预测能力要明显优于MLR模型.

O3LYP;神经网络;多元线性回归;生成热;均方根偏差

1 引言

分子的能量是分子最基本的物理量.在实际研究中无法测定体系精确的能量,因此我们更关心那些能够表征体系能量的热力学函数,其中生成热(ΔfH⊖)就是一个最基本的表征化合物在标准状态下稳定性的热力学量,对判断有机反应的方向和限度也十分重要,并且实验上可以测得其精确值,1-3因此生成热的应用十分广泛.但是实验上要精确测定一个化合物的生成热也是很困难的,不仅需要一套合适的装置,还需要花费大量的时间,并且还需要测定者有很丰富的测定经验,而且不同的计算方法间也会有偏差,4-6有时新测定的结果还会否定旧的实验结果,7因此近年测定新的化合物生成热数据的报导很少,于是发展一个简单的并能准确预测生成热的计算模型是非常必要的.

随着计算机技术的飞速发展,科学工作者对计算化学越来越依赖,从而计算化学已成为对实验结果的解释和预测不可缺少的工具.8与经典从头算方法相比,密度泛函理论(DFT)具有计算花费少,而且可以计算更大的分子,它在计算量子化学领域广泛应用,已经成为化学、凝聚态物理、材料科学和分子生物学中重要的研究工具,但是DFT计算仍然存在一些计算偏差.例如目前最流行的杂化密度泛函B3LYP计算存在一些缺陷:随着系统增大偏差增大;9,10它所计算键的离解焓偏低;11,12得到不准确的异构体能量差别;12低估了反应能垒;13,14不能准确计算分子间的范德华作用力.15Cohen和Handy16在B3LYP的基础上发展了一种新的泛函—O3LYP,他们引入一个交换泛函OPTX.近几年的研究发现O3LYP的计算结果要优于B3LYP.17-23

虽然Gn(n=1-4)方法在计算小分子生成热方面成效非常显著,9,24-29但是它们的计算花费不经济,科学家一般不利用此类方法计算大分子的生成热.另一方面,由于DFT计算的一些固有偏差,其局限性也越来越突出,尤其是在计算大分子体系时更是如此.因此,一些学者提出了一些校正模型对DFT计算结果进行校正,并取得了令人满意的结果. Chen等30,31首次采用ANN校正了密度泛函理论的计算结果,生成热的计算值与实验值之间的均方根偏差从89.5 kJ·mol-1(或50.2 kJ·mol-1)降为13.0 kJ· mol-1(或13.8 kJ·mol-1).Fan等32提出了用自然键轨道(NBO)分析并运用神经网络校正350个分子的计算生成热,减少了350个分子的生成热的计算偏差,随后他们运用多元线性回归算法校正DFT计算的生成热,也成功地减小了DFT的计算偏差.33Li等34利用支持向量机模型也成功地校正了DFT计算的生成热.这说明神经网络,多元线性回归,支持向量机等统计方法是适合于校正DFT计算结果的. Chen35及本课题组36也采用神经网络模型通过调节B3LYP杂化密度泛函中三个参数来减少密度泛函理论计算的分子热力学性质数值同实验值的偏差,并取得了成功.本文利用ANN和MLR的方法来校正O3LYP计算的220个有机分子的生成热,并对这两种校正方法的计算结果进行比较.

2 方法理论

2.1 计算方法

本文选用的是包括小型(重原子数目＜6)和中型(重原子数目≥6)的220个有机分子,所有的分子都是由H、C、N、O、F、S、Cl、Br这几种原子组成,每个分子中含有1-12个C原子,这些分子的生成热实验值都是从实验手册中获得.1-3所有分子在这三本手册中的ΔfH⊖的实验值偏差都小于4 kJ·mol-1.30在上述分子中,选取文献30中的180个分子作为训练集用于构造ANN或MLR模型,利用交叉验证的方法来检验ANN和MLR的校正结果.本文首先用O3LYP/ 6-31G(d)方法优化所有分子的结构,然后用O3LYP/ 6-311+G(3df,2p)方法计算其生成热,作为训练集的180个分子和测试集的40个分子计算结果与实验值的RMS偏差分别是24.7、21.3 kJ·mol-1.对所有分子进行了NBO分析.NBO分析得到的分子的结构参数用于构造ANN(MLR)模型的描述符(自变量).所有O3LYP计算均采用Gaussian 09程序包.37

在构造统计校正模型时,最重要的一点就是关于分子的物理描述符的选择,因为这些描述符被用作神经网络的输入层和多元线性回归的自变量.首先我们选择O3LYP计算的生成热()作为统计模型的第一个描述符.Chen等30在利用神经网络改进B3LYP计算的时,他们发现对神经网络的校正结果非常显著,因此可以作为第一个描述符.另外Gn(n=1-4)在计算中小型分子生成热时效果显著,这是因为在Gn方法中引入了一个高阶校正(HLC)项来抵消一系列计算后剩余能量偏差,26而HLC主要采用了分子和原子中的价层孤对电子和未成对电子对其能量进行最后拟合得到.Gn方法认为大部分的相关能是由占据在同一个分子轨道上的电子之间的相互作用引起的,因此电子越多,它们的相关能就越大,且成对电子对相关能的贡献和未成对电子是不一样的.而且G4对原子和分子采用不同的参数,说明电子在分子环境和原子环境中,它们对相关能的贡献也不一样.鉴于以上分析,本文也采用NBO计算分子的双中心成键电子数(BD)、双中心反键电子数(BD*)、单中心价层孤对电子数(LP)、单中心的内层电子数(CR)这四个参数作为描述符.考虑到分子的大小对的计算偏差有很大的影响,其偏差随着原子数的增加而递增,因而把原子总数(NT)和氢原子个数(NH)作为描述符.分子的零点能(ZPE)是计算的一个重要参数,因此将ZPE作为另一个重要的描述符.30

2.2 多元线性回归模型

采用SPSS11.01统计软件分析180个分子训练集的上述描述符参数与生成热实验值的关联度,并构建实验值的多元线性回归方程如下:

其中,ΔfH⊖为应用MLR方程校正后的生成热值,r为相关系数,RMSD为均方根偏差.

2.3 神经网络模型

神经网络是一个由许多简单的并行工作的处理单元组成的系统,其功能取决于网络的结构、连接强度以及各单元的处理方式.所采用的神经网络模型包含有输入层,隐藏层和输出层三层.本文中的神经网络结构见图1,其中输入层(X1,X2,X3,…)是由8个描述符和一个偏置量(bias)组成,隐藏层包含有3个隐藏神经元(y1,y2,y3),一个输出层是输出校正后的分子生成热.{Wxij}和{Wyj}表示突触权重,其中{Wxij}是连接输入神经元和隐藏神经元的突触权重,{Wyj}是连接隐藏神经元和输出神经元的突触权重.在神经网络模型中所用的转换函数是sigmoid函数.输出层Z与输入层Xi的关系由公式(2)和(3)描述:

其中,

这里Z为经ANN校正后的生成热值,Xi对应图1中给出的8个描述符和一个偏置量,Sig(v)为sigmoid函数.

2.4 模型的验证

本文采用六重交叉验证的方法来验证MLR和ANN模型的校正结果,即将180个分子分成两组,其中一组共有150个分子作为训练集,另一组30个分子作为测试集,这个过程循环重复了6次.同时计算评价模型拟合能力的标准回归系数(q2)和RMSD. Tropsha等38将q2定义为:

其中,Zi和分别表示物种的实验值和预测值,Zi表示物种训练值的平均值,因此本文采用的计算公式如下:

图1 神经网络的结构Fig.1 Structure of the neural networkXi(i=1-9)are ΔfH⊖calc,the total quantities of atoms(NT),hydrogen atoms(NH),2-center bond(BD),2-center antibond(BD*),1-center valence lone pair(LP),and 1-center core pair(CR),the zero point energy(ZPE)and bias,respectively.yj(j=1-3)are the number of hidden neuron.{Wxij}are the synaptic weights connecting the input layer(Xi)and the hidden neurons(yj)and{Wyj}connect the hidden neurons and the output(Z).Z is the ΔfH⊖corrected byANN.

3 结果与讨论

3.1 MLR校正生成热

其中,系数a-g对每次回归分析都有不同的数值,与方程(1)相比,经过交叉验证后得到的r与r2几乎不变,MLR的结果分析见表1.由表1可知,MLR经交叉验证后的q2值是-0.53,训练集的平均RMSD是12.5 kJ·mol-1,但是测试集的平均RMSD高达15.1 kJ·mol-1,总体RMSD是13.0 kJ·mol-1.Cramer等39曾指出q2为负值是因为测试集的平均偏差要大于训练集的平均偏差,说明模型预测能力较差.Golbraikh和Tropsha40曾提出统计参数q2＞0.5,r2＞0.6的模型有效,同时Tropsha等也认为模型的预测能力比拟合能力更为重要,即q2的大小更为关键.由上述实验结果可知,利用MLR来预测中小型分子的生成热是不可靠的.MLR虽然在训练集中可以降低生成热的RMSD,但是其模型的预测能力较差.

3.2 神经网络校正生成热

图2是O3LYP计算结果、MLR校正结果(O3LYP-MLR)、ANN校正结果(O3LYP-ANN)和实验值的对比以及偏差直方图.由图2(a,b)可知,与实验值相比,原始的O3LYP计算结果存在很大的系统偏差,偏差在±4 kJ·mol-1范围内的分子个数只有24个,经过神经网络校正后(见图2(e,f)),生成热的校正结果明显更加接近于实验值,偏差分布非常集中,所有分子的偏差范围都集中在±32 kJ·mol-1的范围内,而且偏差在±4 kJ·mol-1范围内的分子个数有63个,明显好于O3LYP的计算结果.MLR尽管也能在一定程度上对O3LYP的计算结果进行校正,但其模型预测能力较差.

表3 ΔfH⊖的实验值和不同方法的计算值和实验值间的偏差aTable 3 Experimental ΔfH⊖and the deviation between the calculated and experimental values with different methodsa

continued Table 3

表4 测试集ΔfH⊖的实验值和不同方法的计算值和实验值间的偏差aTable 4 Experimental ΔfH⊖of testing set and the deviation between the calculated and experimental values with different methodsa

180个训练集分子经ANN校正后的结果列于表3.为了对比,表3中也列出了O3LYP的计算值与实验值间的偏差,MLR得到的方程(1)拟合的计算值与实验值间的偏差及全部分子的ΔfH⊖实验值.

综上所述,不管是模型的拟合还是模型的预测,利用ANN校正生成热的结果要明显优于MLR模型.原因可能是ANN模型采用了Sigmod函数来校正理论生成热,而MLR模型不是利用非线性函数.在ANN模型里,每个输入的参数都与隐藏层的神经元有一定的关系,而本文所采用的ANN模型的输入层中包含一个偏置量,而且隐藏层包含三个神经元,因此ANN在拟合的过程中产生的变量数目要比MLR多,而在MLR方法中,变量的数目与输入的参数数目相同.

3.3 独立测试集的预测结果

为了验证MLR和ANN模型的预测能力,我们取出40个分子建立了一个独立测试集,经过MLR和ANN模型校正后,测试集分子生成热的RMSD从21.3 kJ·mol-1分别降低到12.1、10.4 kJ·mol-1.因此,ANN模型的预测能力要优于MLR模型.表4列出了独立测试集分子经MLR和ANN校正后的结果,同时也列出了生成热的O3LYP计算值与实验值间的偏差.

4 结论

采用O3LYP/6-311+G(3df,2p)//O3LYP/6-31G(d)计算220个有机分子的生成热,然后采用ΔfH⊖calc、双中心成键电子数、双中心反键电子数、单中心价层孤对电子数、单中心的内层电子数、原子总数和氢原子个数、分子的零点能作为描述符参数,运用神经网络方法和多元线性回归方法对生成热进行校正.研究结果表明:经过神经网络校正后,ΔfH⊖的理论计算值与实验值间的RMSD从24.7 kJ·mol-1降低到11.8 kJ·mol-1,ANN模型的标准回归系数是0.96, MLR虽然在一定程度上降低了生成热的RMSD,但是它的标准回归系数很小,说明MLR的预测能力很差.采用ANN和MLR模型校正40个独立测试集分子的生成热,结果表明其RMSD从21.3 kJ·mol-1分别降低到10.4、12.1 kJ·mol-1.因此ANN模型的拟合和预测能力要明显优于MLR.在利用统计方法校正能量时,描述符的选取是问题的关键,因此有必要选择合适的描述符,从而进一步改善统计方法对理论计算获得的能量的校正结果.

(1) Pedley,J.B.;Naylor,R.D.;Kirby,S.P.Thermochemical Data of Organic Compounds;Chapman and Hall:New York,1986.

(2)Yaws,C.L.Chemical Properties Handbook;McGraw-Hill: New York,1999.

(3) Lide,D.R.CRC Handbook of Chemistry and Physics,3rd. electronic ed.;BocaRaton:FL,2000.

(4) Wu,J.;Xu,X.J.Chem.Phys.2007,127,214105.doi:10.1063/ 1.2800018

(5) Curtiss,L.A.;Raghavachari,K.;Redfern,P.C.;Pople,J.A. Chem.Phys.Lett.1997,270,419.doi:10.1016/S0009-2614(97) 00399-0

(6) Schmitz,L.R.;Chen,K.H.;Labanowski,J.;Allinger,N.L. J.Phys.Org.Chem.2001,14,90.doi:10.1002/1099-1395 (200102)14:2＜90::AID-POC330＞3.0.CO;2-O

(7) Curtiss,L.A.;Raghavachari,K.;Trucks,G.W.;Pople,J.A. J.Chem.Phys.1991,94,7221.doi:10.1063/1.460205

(8) Lado-Tourino,I.;Tsobnang,F.Comp.Mater.Sci.1998,11,181. doi:10.1016/S0927-0256(98)80004-9

(9) Curtiss,L.A.;Raghavachari,K.;Redfern,P.C.;Pople,J.A. J.Chem.Phys.2000,112,7374.doi:10.1063/1.481336

(10) Wodrich,M.D.;Corminboeuf,C.;Schleyer,P.v.R.Org.Lett. 2006,8,3631.doi:10.1021/ol061016i

(11) Check,C.E.;Gilbert,T.M.J.Org.Chem.2005,70,9828.doi: 10.1021/jo051545k

(12) Izgorodina,E.I.;Coote,M.L.;Radom,L.J.Phys.Chem.A 2005,109,7558.doi:10.1021/jp052021r

(13) Schreiner,P.R.;Fokin,A.A.;Pascal,R.A.,Jr.;de Meijere,A. Org.Lett.2006,8,3635.doi:10.1021/ol0610486

(14)Zhao,Y.;Gonzalez-Garcia,N.;Truhlar,D.G.J.Phys.Chem.A 2005,109,2012.doi:10.1021/jp045141s

(15) Zhang,I.Y.;Luo,Y.;Xu,X.J.Chem.Phys.2010,132,194105. doi:10.1063/1.3424845

(16) Cohen,A.J.;Handy,N.C.Mol.Phys.2001,99,607.doi: 10.1080/00268970010023435

(17)Yang,K.;Peverati,R.;Truhlar,D.G.;Valero,R.J.Chem.Phys. 2011,135,044118.doi:10.1063/1.3607312

(18) Heerdt,G.;Morgon,N.H.Quimica Nova 2011,34,868.doi: 10.1590/S0100-40422011000500024

(19) Bochevarov,A.D.;Friesner,R.A.;Lippard,S.J.J.Chem. Theory Comput.2010,6,3735.doi:10.1021/ct100398m

(20)Dobado,J.A.;Gomez-Tamayo,J.C.;Calvo-Flores,F.G.; Martinez-Garcia,H.;Cardona,W.;Weiss-Lopez,B.;Ramirez-Rodriguez,O.;Pessoa-Mahana,H.;Araya-Maturana,R.Magn. Reson.Chem.2011,49,358.doi:10.1002/mrc.2745

(21) Qian,Z.S.;Feng,H.;He,L.N.;Yang,W.J.;Bi,S.P.J.Phys. Chem.A 2009,113,5138.doi:10.1021/jp810632f

(22) Strassner,T.;Taige,M.A.J.Chem.Theory Comput.2005,1, 848.doi:10.1021/ct049846+

(23) Baker,J.;Pulay,P.J.Comput.Chem.2003,24,1184.doi: 10.1002/jcc.10280

(24) Curtiss,L.A.;Jones,C.;Trucks,G.W.;Raghavachari,K.; Pople,J.A.J.Chem.Phys.1990,93,2537.doi:10.1063/ 1.458892

(25) Curtiss,L.A.;Raghavachari,K.;Redfern,P.C.;Pople,J.A. J.Chem.Phys.1997,106,1063.doi:10.1063/1.473182

(26) Curtiss,L.A.;Raghavachari,K.;Redfern,P.C.;Rassolov,V.; Pople,J.A.J.Chem.Phys.1998,109,7764.doi:10.1063/ 1.477422

(27) Curtiss,L.A.;Raghavachari,K.;Redfern,P.C.;Pople,J.A. J.Chem.Phys.2000,112,7374.doi:10.1063/1.481336

(28) Curtiss,L.A.;Redfern,P.C.;Raghavachari,K.Chem.Phys. Lett.2010,499,168.doi:10.1016/j.cplett.2010.09.012

(29) Dorofeeva,O.V.;Kolesnikova,I.N.;Marochkin,I.I.;Ryzhova, O.N.J.Struct.Chem.2011,22,1303.doi:10.1007/s11224-011-9827-7

(30)Hu,L.H.;Wang,X.J.;Wong,L.H.;Chen,G.H.J.Chem. Phys.2003,119,11501.doi:10.1063/1.1630951

(31)Wang,X.J.;Wong,L.H.;Hu,L.H.;Chan,C.Y.;Su,Z.M.; Chen,G.H.J.Phys.Chem.A 2004,108,8514.doi:10.1021/ jp047263q

(32)Duan,X.M.;Li,Z.H.;Song,G.L.;Wang,W.N.;Chen,G.H.; Fan,K.N.Chem.Phys.Lett.2005,410,125.doi:10.1016/ j.cplett.2005.05.046

(33)Duan,X.M.;Song,G.L.;Li,Z.H.;Wang,X.J.;Chen,G.H.; Fan,K.N.J.Chem.Phys.2004,121,7086.doi:10.1063/ 1.1786582

(34) Yan,G.K.;Li,J.J.;Li,B.R.;Hu,J.;Guo,W.P.J.Theor. Comput.Chem.2007,6,495.doi:10.1142/ S0219633607003118

(35) Zheng,X.;Hu,L.;Wang,X.;Chen,G.Chem.Phys.Lett.2004, 390,186.doi:10.1016/j.cplett.2004.04.020

(36) Zhang,J.H.;Wang,X.J.Acta Phys.-Chim.Sin.2010,26, 188. [张家虎,王秀军.物理化学学报,2010,26,188.] doi:10.3866/PKU.WHXB20100116

(37) Frisch,M.J.;Trucks,G.W.;Schlegel,H.B.;et al.Gaussian 09, Revision B.01;Gaussian Inc.:Wallingford,CT,2010.

(38) Zheng,W.F.;Tropsha,A.J.Chem.Inf.Comput.Sci.2000,40, 185.doi:10.1021/ci980033m

(39) Cramer,R.D.,III;Patterson,D.E.;Bunce,J.D.J.Am.Chem. Soc.1988,110,5959.doi:10.1021/ja00226a005

(40) Golbraikh,A.;Tropsha,A.J.Chem.Inf.Comput.Sci.2003,43, 144.doi:10.1021/ci025516b

April 9,2012;Revised:July 16,2012;Published on Web:July 17,2012.

Statistical Correction of Heat of Formation Calculated by the O3LYP Method

WANG Xiu-Jun*LONG Mi
(Key Laboratory of Fuel Cell Techology of Guangdong,College of Chemistry and Chemical Engineering, South China University of Technology,Guangzhou 510640,P.R.China)

The results of density functional theory calculations are known to contain inherent numerical errors caused by various intrinsic approximations.In this paper,O3LYP/6-311+G(3df,2p)//O3LYP/6-31G(d) calculations were used to derive the heats of formation(ΔfH⊖calc)of 220 small to medium-sized organic molecules,followed by the application of artificial neural network(ANN)and multiple linear regression (MLR)analyses to correct the values.The physical descriptors chosen were ΔfH⊖calcand zero point energy as well as the total quantities of atoms,hydrogen atoms,2-center bonds,2-center antibonds,1-center valence lone pairs and 1-center core pairs.The ANN and MLR systems were initially constructed using a 180 training set.The trained ANN and MLR systems were subsequently used to predict values of ΔfH⊖calcfor a 40 individual testing set.The results demonstrated that the root mean square(RMS)deviations between the calculated and experimental ΔfH⊖values in the training set were reduced from 24.7 to 11.8 and 13.0 kJ· mol-1after ANN and MLR corrections,respectively.For the individual testing set,the deviations(RMSD) were reduced from 21.3 to 10.4 and 12.1 kJ·mol-1,respectively.Based on these results,it can be concluded thatANN exhibits superior fitting and predictive abilities compared with MLR.

O3LYP;Neural network;Multiple linear regression;Heat of formation;Root mean square deviation

10.3866/PKU.WHXB201207172

∗Corresponding author.Email:xjwangcn@scut.edu.cn;Tel:+86-20-87112943.

The project was supported by the National Natural Science Foundation of China(20975040)and Natural Science Foundation of Guangdong Province,China(10351064101000000).

国家自然科学基金(20975040)及广东省自然科学基金(10351064101000000)资助项目

O641