結(jié)合神經(jīng)網(wǎng)絡(luò)方法和擴大訓(xùn)練基組構(gòu)建新B3LYP泛函

張家虎 王秀軍

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

結(jié)合神經(jīng)網(wǎng)絡(luò)方法和擴大訓(xùn)練基組構(gòu)建新B3LYP泛函

張家虎 王秀軍*

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

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

B3LYP泛函;神經(jīng)網(wǎng)絡(luò);描述符;訓(xùn)練基組

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

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

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

1 計算方法

雜化B3LYP泛函可以表示為:高階項相關(guān)能是系統(tǒng)相關(guān)的,并且高階項相關(guān)能在式(1)中各項均有體現(xiàn),因此,式(1)中的三參數(shù)也應(yīng)該是系統(tǒng)相關(guān)的.三個參數(shù)可以通過神經(jīng)網(wǎng)絡(luò)得到.

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

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

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

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

其中,

圖1 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖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.

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

2 結(jié)果與討論

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

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

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

表1 神經(jīng)網(wǎng)絡(luò)修正前后的均方根偏差Table 1 RMS errors before and after the neural network correction

表2 最優(yōu)化的突觸權(quán)重Wij和Table 2 Optimized synaptic weightages Wijand jk

表2 最優(yōu)化的突觸權(quán)重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 三個參數(shù)對各個描述符(xi)的導(dǎo)數(shù)Table 3 Derivatives ofwith respect to each physical descriptor(xi)

表3 三個參數(shù)對各個描述符(xi)的導(dǎo)數(shù)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

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

3 結(jié) 論

本文中,我們使用神經(jīng)網(wǎng)絡(luò)方法對B3LYP三參數(shù)進行校正,并借助神經(jīng)網(wǎng)絡(luò)構(gòu)建新B3LYP泛函.采用gs、Ne、D、Ek、Q和ZPE作為神經(jīng)網(wǎng)絡(luò)的描述符,并且以246個能量數(shù)據(jù)作為訓(xùn)練集以確定神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)和突觸權(quán)重值.新的突觸權(quán)重用于構(gòu)建新B3LYP泛函,新B3LYP泛函再用于預(yù)測分子的性質(zhì).通過采用B3LYP/6-311+G(3df,2p)和B3LYP/ 6-311+G(3df,2p)-NN方法分別計算56個AE、42個IP、8個PA、10個TAE以及180個ΔfH?,計算結(jié)果與實驗結(jié)果間的均方根偏差分別為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.此神經(jīng)網(wǎng)絡(luò)對AE和PA的結(jié)果不能有效修正,原因可能是數(shù)據(jù)庫中生成熱數(shù)據(jù)較多以及生成熱數(shù)值原始偏差較大,但對IP、TAE和ΔfH?起到很好的修正作用,尤其是對于ΔfH?的修正效果更為明顯.應(yīng)用此神經(jīng)網(wǎng)絡(luò)可以有效地預(yù)測分子的ΔfH?,同時,此神經(jīng)網(wǎng)絡(luò)將數(shù)據(jù)庫從116個擴大到296個,增大了此神經(jīng)網(wǎng)絡(luò)的應(yīng)用范圍.

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)資助項目