結合神經網絡方法和擴大訓練基組構建新B3LYP泛函

張家虎 王秀軍

(華南理工大學化學與化工學院應用化學系,廣州 510640)

結合神經網絡方法和擴大訓練基組構建新B3LYP泛函

張家虎 王秀軍*

(華南理工大學化學與化工學院應用化學系,廣州 510640)

神經網絡方法成功地應用于修正密度泛函理論B3LYP方法中的三個參數(a0、ax和ac)以構建新B3LYP交換相關泛函.本文采用包含輸入層、隱藏層和輸出層的三層式神經網絡結構.總電子數、多重度、偶極矩、動能、四極矩和零點能被選為物理描述符.296個能量數據被隨機地分成兩組,246個能量數據作為訓練集以確定神經網絡的最優結構和最優突觸權重,50個能量數據作為測試集以測試神經網絡的預測能力.修正后的三個參數?0、?x、?c從輸出層處得到,并用于計算體系的熱化學性質如原子化能(AE)、電離勢(IP)、質子親合能(PA)、總原子能(TAE)和標準生成熱(ΔfH?).修正后的計算結果優于傳統B3LYP/6-311+G(3df,2p)方法的計算結果.經過神經網絡修正后,296個物種的總體均方根偏差從41.0 kJ·mol-1減少到14.2 kJ·mol-1.

B3LYP泛函;神經網絡;描述符;訓練基組

密度泛函理論(DFT)以其良好的計算性能和較低的計算成本,成為化學家和物理學家理解原子、分子、固體及其相關電子結構的一個非常寶貴的工具[1],并且成為化學、凝聚態物理、材料科學和分子生物學中重要的研究工具[2,3].盡管DFT在實際應用中取得了很大的成功,但是它通常并不能足夠準確地預測分子的實驗值,尤其是對于大分子體系,偏差產生的原因就在于DFT中的固有近似[4].DFT計算結果的準確度是由其精確交換相關泛函XC泛函決定的[2],而精確交換相關泛函是不知道的.所有的DFT計算均使用近似交換相關泛函,這就進一步擴大了計算結果的偏差.所以尋找更加精確的交換相關泛函對提高DFT計算精度有著重要的意義.

最初的近似交換相關泛函是基于局域梯度近似(LDA)[5],后來廣義梯度近似(GGA)[6-9]和雜化泛函[10]的發展對交換相關泛函的發展起了很大的促進作用.交換相關泛函包括局域和非局域貢獻,除此之外,高階項相關能對交換相關泛函的貢獻也不可忽視.Chen等[11]討論了高階項相關能對B3LYP泛函的貢獻.近年來交換相關泛函的一個研究趨勢就是利用準確的實驗數據確定交換相關泛函表達式中的參數值.例如:Becke使用最小二乘法擬合G2[12,13]基組中原子和分子能量數據,確定了B3LYP泛函表達式中的三參數值[10].

B3LYP泛函在過去的十幾年中得到了重大的發展,尤其是在使用B3LYP預測小型和中型分子體系的熱化學性質方面[14-16].但是B3LYP計算結果與實驗結果之間也存在著較大的偏差.因此,許多學者提出了一些修正模型,如線性回歸模型[17,18]、支持向量機模型[19]和神經網絡模型[20-25],對B3LYP方法進行修正,極大地改進了計算結果的偏差.這些模型只是對DFT計算結果進行修正,而文獻[11]提出了用神經網絡來構建B3LYP泛函.本文在此基礎上將數據庫從116個數據擴展到了296個(116個來自文獻[10],180個來自文獻[20],文獻[10,20]已經考證這些數據的實驗值是可靠的),大大地提高了此方法的適用性.

1 計算方法

雜化B3LYP泛函可以表示為:高階項相關能是系統相關的,并且高階項相關能在式(1)中各項均有體現,因此,式(1)中的三參數也應該是系統相關的.三個參數可以通過神經網絡得到.

采用神經網絡進行計算時,首先要構建有效的神經網絡結構,而構建有效的神經網絡結構就要選擇能充分反映體系特征性質的物理描述符.由于密度泛函理論是基于體系電子密度的理論,而交換相關泛函也只是電子密度的泛函,因此選擇的物理描述符要能反映體系的電子特性.同時,選擇的物理描述符要能夠反映體系電子的分布,因為電子的分布和電子密度緊密相關.在選取了物理描述符之后,我們就隨機選用246個能量數據作為訓練集,50個能量數據作為測試集,用訓練集確定神經網絡的結構,并用此神經網絡確定三參數值的校正值.同時,得到新的基于神經網絡的B3LYP交換相關泛函.

Perdew等[28,29]在GGA基礎上提出了meta-GGA,在meta-GGA中,交換相關泛函不僅是電子密度和密度梯度的泛函,同時也是動能密度的泛函,因此動能Ek被選取為關鍵物理描述符.電子密度分布ρ(r)唯一地決定著交換相關泛函并且可以從多極矩的角度被擴展,作為零階項的擴展,總電子個數Ne也被選作是物理描述符.體系的偶極矩和四極矩也是如此,我們選擇作為偶極描述符,Q≡作為四極描述符,其中di(i=x,y,z)是偶極矩向量的分量,Qii(i=x,y,z)是四極矩張量的對角分量.因為交換泛函同自旋電子之間的交換作用相關,所以多重度gs也被選作是物理描述符.零點能(ZPE)對能量的校正起著重要的作用,因此了也被選作關鍵物理描述符.

我們采用三層式的神經網絡結構:包括輸入層、隱藏層和輸出層(見圖1).相鄰層之間通過突觸權重(Wij或)連接.物理描述符gs、Ne、D、Ek、Q和 ZPE的值從輸入層輸入神經網絡.對于每一個原子或分子而言,三參數a0、ax和ac修正后(修正后的三參數表示為?0,?x和?c)的值可以從輸出層獲得.我們使用246個能量實驗值作為訓練集,并用它來確定神經網絡的結構和突觸權重值.我們采用B3LYP/6-311+G(3df,2p)進行計算,每個分子的幾何構型直接用B3LYP/6-311+G(d,p)方法優化得到.D、Ek、Q和ZPE的值也由B3LYP/6-311+G(3df,2p)方法計算得出.除了物理描述符gs、Ne、D、Ek、Q和ZPE的值輸入輸入層之外,偏置量也被引入了輸入層和隱藏層,并且偏置量在整個訓練過程中保持不變.同時,采用誤差反向傳播學習方法優化突觸權重.對于每一個分子或原子而言,系統相關的三參數?0、?x和?c被用來修正相應的B3LYP泛函,修正后的B3LYP泛函依次用來計算相應的原子化能(AE)、電離勢(IP)、質子親合能(PA)、總原子能(TAE)和標準生成熱(ΔfH?)的值.修正后的能量值同相應的實驗值比較,誤差反饋給神經網絡,同時,所有的突觸權重值獲得相應地調整.這個步驟重復進行直到整個訓練集能量的計算值和實驗值足夠的接近.然后,整個神經網絡就可被認為是收斂的,突觸權重也相應最終確定.

我們將296個能量數據(56個用于計算AE、42個用于計算IP、8個用于計算PA、10個用于計算TAE以及180個用于計算ΔfH?)隨機地分成大小相近的六組.其中五個用于訓練神經網絡的權重,稱之為訓練組.第六組用于檢驗神經網絡的預測效果,稱之為驗證組.如此重復六次以優化神經網絡的結構.我們從1到10改變隱藏神經元的數目,結果表明,8個隱藏神經元效果最好.即:總體均方根偏差以及訓練組與驗證組的均方根偏差最小.預測組預測結果確保了神經網絡的預測效果.因此,我們采用了7-9-3結構(見圖1).為了使輸出數據同目標數據易于比較,需對輸入層描述符(xi)數值進行預處理.除了偏置量之外,每一個xi的值都被換算到區間(0.1, 0.9)內.經神經網絡修正的?0、?x和?c值在輸出層處獲得,并且它們同{xi}之間的關系如下:

其中,

圖1 神經網絡結構圖Fig.1 Architectural graph of the neural networkxi(i=1-6)are spin multiplicity(gs),the total number of electrons(Ne),dipole moment(D),kinetic energy(Ek),quadrupole moment(Q),and zero point energy(ZPE),respectively.yj(j=1-8)are the number of hidden neuron.x0and y0are bias in input layer and in hidden layer,respectively.{Wij}are the synaptic weights connecting the input layer(xi)and the hidden neurons(yj),and{Wj′k}connect the hidden neurons and the output.?0,?xand ?care the modified three parameters in B3LYP method.AE,IP,PA,TAE,and ΔfH?are atomic energy,ionization potential,proton affinity,total atomic energy, and standard heat of formation calculated with B3LYP,respectively.

參數α和γ分別控制S形函數Siga(v)和Sigb(v)的陡度,β是控制Sigb(v)彎曲程度的參數.Wij是連接輸入層和隱藏層的突觸權重,是連接隱藏層和輸出層的突觸權重(i=0-6,j=0-8,0代表輸入層和隱藏層的偏置量項).

2 結果與討論

采用B3LYP/6-311+G(3df,2p)計算56個物種的AE、42個物種的IP、8個物種的PA、10個物種的TAE以及180個物種的ΔfH?,計算結果與實驗結果間的均方根偏差分別為12.6、20.5、6.7、43.1和50.2 kJ·mol-1(見表1).

訓練集中每一個分子或原子的物理描述符gs、Ne、D、Ek、Q和ZPE的值從輸入層輸入神經網絡,輸出層得到修正的三參數新的三參數用于構建新B3LYP泛函,然后再將新B3LYP泛函用于計算AE、IP、PA、TAE和ΔfH?.這些值再同訓練集中各物種的能量值進行比較以調整突觸權重{Wij}和.最優化的突觸權重值顯示在表2中.

在xi=0.5(i=1-6)和x0=1(偏置量)時對xi(i=0-6)的導數值列于表3中,其值的大小反映了相應物理描述符對確定值的貢獻,該數值越大,相應的物理描述符對確定?0、?x和?c的值越重要.我們發現gs、Ne、D、Ek和Q對?0、?x和?c的貢獻都很大,這與文獻[11]中選用的描述符相同,說明我們選擇的描述符比較合適.因此,gs、Ne、D、Ek和Q被認為是五個最重要的物理描述符以確定?0、?x和?c.訓練集中所有的分子、原子和相關離子?0、?x和?c的平均值是0.728、0.685和0.852,這不同于使用有限基組得到的初始B3LYP三參數值.對于訓練組中每一個分子或原子而言,?x＜?0＜?c始終成立,更重要的是,?0、?x和?c的值彼此之間略有差異.即:相應的B3LYP泛函是系統相關的.

表1 神經網絡修正前后的均方根偏差Table 1 RMS errors before and after the neural network correction

表2 最優化的突觸權重Wij和Table 2 Optimized synaptic weightages Wijand jk

表2 最優化的突觸權重Wij和Table 2 Optimized synaptic weightages Wijand jk

i=0-6:for bias and input neurons;j=0-8:for bias and hidden neurons; k=1-3:for modified three parameters in output layer

j=1 j=2 j=3 j=4 j=5 j=6 j=7 j=8 j=0 i=1 0.778 3.361 0.533-0.530 1.132-0.855 3.167 0.205 -i=2 1.135 0.166 0.812-0.882-0.803-0.080-0.377 0.772 -i=3-0.840-0.588 0.970 2.317-1.474-3.453 1.805-1.888 -i=4 1.785-2.250 1.891 0.153 0.871 0.200-0.632-1.669 -i=5 1.023-0.161-0.137-2.072 2.247 4.274 3.214-0.008 -i=6-0.151 0.484-3.460-3.929-2.167 0.175-1.060-0.630 -i=0-0.598-1.737-0.639 0.202-0.220-1.002-1.628 0.475 -k=1-0.400 0.131 0.232 0.096-0.143 0.262-0.049 0.416-0.131 k=2 0.265 0.070 0.067 0.376 0.290-0.054-0.043 0.070-0.661 k=3 0.467 0.358 0.180-0.186 0.136-0.351 0.397 0.020 0.020

表3 三個參數對各個描述符(xi)的導數Table 3 Derivatives ofwith respect to each physical descriptor(xi)

表3 三個參數對各個描述符(xi)的導數Table 3 Derivatives ofwith respect to each physical descriptor(xi)

i=0:for bias;i=1-6:for input neurons(gs,Ne,D,Ek,Q,and ZPE).

i=1 i=2 i=3 i=4 i=5 i=6 i=0??0/?xi 0.012 0.043 -0.014 -0.031 0.025 -0.023 -0.070??x/?xi 0.065 0.054 0.027 0.056 0.100 0.007 -0.011??c/?xi 0.091 0.103 0.146 0.135 0.084 0.039 0.015

為了檢驗神經網絡的性能,我們從原有數據庫中隨機抽取了50個數據作為測試組,用以對確定的神經網絡進行測試.測試組不同于訓練組中的分子或原子,它主要用于對神經網絡的性能進行檢測.在測試階段神經網絡可以測出任何新分子的?0、?x和?c,因為神經網絡結構和所有的突觸權重已由前面的訓練過程所確定.新分子物理描述符的值可以通過傳統的B3LYP計算得到,然后輸入最優神經網絡的輸入層,經神經網絡修正的?0、?x和?c在輸出層給出.修正后的三參數用于構建新B3LYP泛函以計算體系的特定性質.測試組50個數據同實驗值的均方根偏差為18.4 kJ·mol-1.

3 結 論

本文中,我們使用神經網絡方法對B3LYP三參數進行校正,并借助神經網絡構建新B3LYP泛函.采用gs、Ne、D、Ek、Q和ZPE作為神經網絡的描述符,并且以246個能量數據作為訓練集以確定神經網絡的結構和突觸權重值.新的突觸權重用于構建新B3LYP泛函,新B3LYP泛函再用于預測分子的性質.通過采用B3LYP/6-311+G(3df,2p)和B3LYP/ 6-311+G(3df,2p)-NN方法分別計算56個AE、42個IP、8個PA、10個TAE以及180個ΔfH?,計算結果與實驗結果間的均方根偏差分別為12.6、20.5、6.7、43.1和50.2 kJ·mol-1和15.9、16.7、23.4、25.5和11.3 kJ·mol-1,總的均方根偏差從41.0 kJ·mol-1減小到14.2 kJ·mol-1.此神經網絡對AE和PA的結果不能有效修正,原因可能是數據庫中生成熱數據較多以及生成熱數值原始偏差較大,但對IP、TAE和ΔfH?起到很好的修正作用,尤其是對于ΔfH?的修正效果更為明顯.應用此神經網絡可以有效地預測分子的ΔfH?,同時,此神經網絡將數據庫從116個擴大到296個,增大了此神經網絡的應用范圍.

1 Koch,W.;Holthausen,M.C.A chemist′s guide to density functional theory.2nd ed.Weinheim:Wiley-VCH,2001:117-188

2 Parr,R.G.;Yang,W.Density-functional theory of atoms and molecules.New York:Oxford University Press,1989:285-317

3 Cramer,C.J.Essentials of computational chemistry:theories and models.2nd ed.West Sussex:John Wiley&Sons Ltd.,2004:355-424

4 Hu,L.H.;Wang,X.J.;Wong,L.H.;Chen,G.H.J.Chem.Phys., 2003,119:11501

5 Ceperley,D.M.;Alder,B.J.Phys.Rev.Lett.,1980,45:566

6 Becke,A.D.Phys.Rev.A,1988,38:3098

7 Lee,C.;Yang,W.;Parr,R.G.Phys.Rev.B,1988,37:785

8 Perdew,J.P.;Wang,Y.Phys.Rev.B,1992,45:13244

9 Perdew,J.P.;Chevary,J.A.;Vosko,S.H.;Jackson,K.A.; Pederson,M.R.;Singh,D.J.;Fiolhais,C.Phys.Rev.B,1992,46: 6671

10 Becke,A.D.J.Chem.Phys.,1993,98:5648

11 Zheng,X.;Hu,L.H.;Wang,X.J.;Chen,G.H.Chem.Phys.Lett., 2004,390:186

12 Curtiss,L.A.;Raghavachari,K.;Trucks,G.W.;Pople,J.A. J.Chem.Phys.,1991,94:7221

13 Curtiss,L.A.;Raghavachari,K.;Redfern,P.C.;Pople,J.A. J.Chem.Phys.,1997,106:1063

14 Jung,D.;Chen,C.J.;Bozzelli,J.W.J.Phys.Chem.A,2000,104: 9581

15 Chen,C.C.;Lay,T.H.;Bozzelli,J.W.J.Phys.Chem.A,2003, 107:6451

16 Sabzyan,H.;Omrani,A.J.Mol.Struct.-Theochem,2005,713:43

17 Duan,X.M.;Song,G.L.;Li,Z.H.;Wang,X.J.;Chen,G.H.;Fan, K.N.J.Chem.Phys.,2004,121:7086

18 Duan,X.M.;Li,Z.H.;Hu,H.R.;Song,G.L.;Wang,W.N.;Chen, G.H.;Fan,K.N.Chem.Phys.Lett.,2005,409:315

19 Yan,G.K.;Li,J.J.;Li,B.R.;Hu,J.;Guo,W.P.J.Theor.Comput. Chem.,2007,6:495

20 Hu,L.H.;Wang,X.J.;Wong,L.H.;Chen,G.H.J.Chem.Phys., 2003,119:11501

21 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

22 Duan,X.M.;Li,Z.H.;Song,G.L.;Wang,W.N.;Chen,G.H.; Fan,K.N.Chem.Phys.Lett.,2005,410:125

23 Zhao,J.;Cheng,X.L.;He,B.;Yang,X.D.Struct.Chem.,2006, 17:501

24 Li,H.;Shi,L.L.;Zhang,M.;Su,Z.;Wang,X.J.;Hu,L.H.;Chen, G.H.J.Chem.Phys.,2007,126:144101

25 Wu,J.;Xu,X.J.Chem.Phys.,2007,127:214105

26 Kohn,W.;Sham,L.J.Phys.Rev.A,1965,140:1133

27 Vosko,S.H.;Wilk,L.;Nusair,M.Can.J.Phys.,1980,58:1200

28 Perdew,J.P.;Kurth,S.;Zupan,A.;Blaha,P.Phys.Rev.Lett., 1999,82:2544

29 Perdew,J.P.;Kurth,S.;Zupan,A.;Blaha,P.Phys.Rev.Lett., 1999,82:5179

August 4,2009;Revised:October 17,2009;Published on Web:November 24,2009.

Neural Network Approach for a New B3LYP Functional with an Enlarged Training Set

ZHANG Jia-Hu WANG Xiu-Jun*
(Department of Applied Chemistry,College of Chemistry and Chemical Engineering,South China University of Technology, Guangzhou 510640,P.R.China)

A neural network approach was used to correct three parameters(a0,ax,and ac)in the B3LYP method for constructing a new B3LYP exchange correlation functional.A three-layer architecture which consisted of an input layer,a hidden layer,and an output layer,was adopted in the neural network.The total number of electrons,spin multiplicity,dipole moment,kinetic energy,quadrupole moment,and zero point energy were chosen as the most important physical descriptors.In this work,296 energy values were randomly divided into two subsets:246 energy values were used as the training set to determine the optimized structure of the neural network and the optimized synaptic weights;50 energy values were used as a testing set to test the prediction capability of our neural network. Three modified parameters ?0,?x,and ?cthat were obtained from the output layer were used to calculate thermochemical data such as the atomic energy(AE),ionization potential(IP),proton affinity(PA),total atomic energy (TAE),and standard heat of formation(ΔfH?).The new results obtained,based on the neural network approach,are much better than the results calculated using the conventional B3LYP/6-311+G(3df,2p)method.Upon the neural network correction,the overall root-mean-square(RMS)error for the 296 species decreased from 41.0 to 14.2 kJ·mol-1.

B3LYP functional;Neural network; Descriptor; Training set

O641

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

The project was supported by the National Natural Science Foundation of China(20975040).

國家自然科學基金(20975040)資助項目