呂志濤,劉玉存,聶江穩(wěn),賈康輝,李 強

(1.中北大學 環(huán)境與安全工程學院,山西 太原 030051;2.山西興新安全生產(chǎn)技術(shù)服務(wù)有限公司,山西 太原 030024)

引 言

含能材料是武器系統(tǒng)實現(xiàn)推進、毀傷的化學能源,對其研究包括分子、晶體、配方和應(yīng)用等方面,其中含能分子的設(shè)計和合成是含能材料研究和發(fā)展的基礎(chǔ)。為了適應(yīng)現(xiàn)代武器高安全、高可靠和高效毀傷的發(fā)展要求,鈍感及高能量密度成為含能材料追求的主要目標。氮、氧雜環(huán)含能材料由于具有高正生成焓,密度大,燃燒產(chǎn)物多為無毒氣體等優(yōu)點[1],是設(shè)計和合成鈍感高能含能材料理想的結(jié)構(gòu)單元,其中雜環(huán)結(jié)構(gòu)單元包括吡唑、咪唑、三唑、四唑、呋咱、吡啶、三嗪、四嗪環(huán)等,通過引入取代基如:硝基、氨基、硝酸酯基、硝胺基、疊氮基、肼基等含能基團,另有偶氮基等連接單元,可以構(gòu)建出大量含能材料分子。如果盲目合成這類含能材料再進行篩選則非常困難,預(yù)先進行能量特性預(yù)測再進行目標化合物合成可大大減小合成工作量。

密度是含能材料重要的能量特性參數(shù),與其他參數(shù)如爆速、爆壓等有密切關(guān)系。因此,建立精準的密度預(yù)測模型對含能分子進行預(yù)測篩選,可以剔除性能不佳的含能材料的研制工作,避免不必要的人力、時間和資金的浪費。目前,含能材料密度的理論預(yù)測方法主要有等電子密度法、分子表面靜電勢法、基團加和法、定量構(gòu)效關(guān)系法等[2]。每種方法各有優(yōu)勢與局限性,運用時需綜合考慮各方面因素。

最早提出預(yù)測密度的方法是基團加和法[3],將原子和官能團的體積相加得到分子摩爾體積。此法簡單但是無法區(qū)分組成相同而分子構(gòu)型、構(gòu)象不同對密度的影響。等電子密度面法是利用量子化學方法計算出0.001a.u.等電子密度面包圍的體積作為分子體積,進而由摩爾質(zhì)量與分子體積的比值求出晶體密度。邱玲、肖鶴鳴等[4]采用此法,比較了45種硝銨類化合物在不同基組下的密度。Rice等[5]研究了180個中性分子和38個高氮分子,均方根誤差在4%以內(nèi)。然而在某些情況下誤差較大,可能是沒有考慮分子間相互作用。故Politzer等[6]將分子表面靜電勢引入,對等電子密度法進行校正,誤差有顯著的改進。Rice等[7]運用Politzer[6]提出的修正方法,用之前研究的化合物驗證得出該校正方法提高了預(yù)測精度。因此分子表面靜電勢法在分子設(shè)計與性能預(yù)測方面得到了廣泛的應(yīng)用。

最初構(gòu)建結(jié)構(gòu)與性能定量關(guān)系的計算方法是以多元線性原理為基礎(chǔ)[8]。馮雪艷[9]運用多元線性回歸對36種多硝基芳香族化合物,以官能團個數(shù)為自變量,對密度進行預(yù)測;來蔚鵬[10]選取最弱鍵長、旋轉(zhuǎn)鍵、氫鍵供體、中點勢、分子總能量等描述符對30種芳香系單質(zhì)炸藥利用多元線性回歸建立密度預(yù)測模型;Narges Zohari等[11]將唑類含能材料的元素組成、化合物不飽和度、硝基亞胺基以及一些結(jié)構(gòu)性參數(shù)的貢獻值建立了一種新的回歸模型。此方法原理較簡單,但需分析所選取的自變量與因變量的線性相關(guān)性,且各自變量間不存在極強的多重共線性等適用條件。近年來,機器學習逐漸被用來預(yù)測含能材料的性質(zhì),它可以對已有數(shù)據(jù)特征進行提取和建模,預(yù)測未知含能材料分子的性能。其中,神經(jīng)網(wǎng)絡(luò)算法因其自身的特點和優(yōu)越性,應(yīng)用愈發(fā)廣泛,可以用作分類、聚類、預(yù)測等。而BP神經(jīng)網(wǎng)絡(luò),是一種按誤差逆向傳播算法訓(xùn)練的多層前饋網(wǎng)絡(luò),結(jié)構(gòu)簡單,是目前應(yīng)用最廣泛的神經(jīng)網(wǎng)絡(luò)模型之一,可用于函數(shù)擬合、預(yù)測等方面。于國強等[12-16]運用神經(jīng)網(wǎng)絡(luò)分別建立了撞擊感度、炸藥臨界直徑、熱感度、生成焓、密度、爆速的預(yù)測模型。上述文獻研究,模型都達到了一定的預(yù)測精度,但是樣本數(shù)量較少,適用范圍有一定的局限性。因此,針對數(shù)據(jù)集小,一種特殊的機器學習——深度學習,越來越受到重視。陳杰[17]從劍橋晶體結(jié)構(gòu)數(shù)據(jù)庫中篩選出2000條只含有硝基的化合物數(shù)據(jù),分別基于隨機森林和圖神經(jīng)網(wǎng)絡(luò)建立了含能材料密度預(yù)測模型;陳超等[18]基于空間矩陣特征化的方法建立了精確的機器學習模型,降低了預(yù)測的平均絕對誤差;Alex等[19]采用分子三維電子結(jié)構(gòu)為描述符,基于卷積神經(jīng)網(wǎng)絡(luò)得到含能材料包括晶體密度的8種主要性質(zhì)的預(yù)測模型;Phan Nguyen等[20]將高能化合物化學結(jié)構(gòu)轉(zhuǎn)化為RDKit二維分子描述符作為輸入,以信息傳遞神經(jīng)網(wǎng)絡(luò)模型預(yù)測密度。這些深度學習算法都獲得了較高的預(yù)測精度,但是需要大量的數(shù)據(jù)集,經(jīng)過特征編碼以及矩陣運算,對專業(yè)知識及硬件要求非常高。

氮、氧雜環(huán)含能材料作為合成高能量密度材料的一種選擇,有必要對這類含能材料建立較為精準的密度預(yù)測模型,為新型雜環(huán)含能化合物的篩選提供理論依據(jù)。鑒于已合成的氮、氧雜環(huán)含能材料的數(shù)量有限,且結(jié)構(gòu)參數(shù)和性能參數(shù)的匱乏,不適合選用深度學習算法,因此本研究采用與文獻[6-7]相同的泛函以及M06-2X泛函,在原有參數(shù)下添加其他合適的結(jié)構(gòu)參數(shù),分別采用多元線性回歸(MLR)和BP神經(jīng)網(wǎng)絡(luò)建立密度預(yù)測模型,以比較不同泛函對密度預(yù)測精度的影響,并與文獻[6-7]方法進行比較,得到一種分子結(jié)構(gòu)與密度的關(guān)聯(lián)模型,提高預(yù)測的準確度和可靠性,為新一代高性能含能分子設(shè)計提供理論指導(dǎo)。

1 實 驗

1.1 樣本數(shù)據(jù)計算與收集

本研究雜環(huán)含能材料樣本數(shù)據(jù)通過查詢《含能材料及相關(guān)物手冊》《火炸藥手冊》《四唑類富氮化合物研究》等書籍及與含能材料相關(guān)的期刊文獻,收集到雜環(huán)含能材料共112種,其中包括10種吡唑類、6種咪唑類、18種三唑類、24種四唑類、29種呋咱類、9種三嗪類、8種四嗪類和8種吡啶類,并從手冊和CCDC數(shù)據(jù)庫中查詢各物質(zhì)標準狀態(tài)下穩(wěn)定晶型的密度實驗值,將其記錄并保存。

1.2 研究方法

計算和查詢到的共10種結(jié)構(gòu)參數(shù)在建立模型時需進行變量的篩選來確定引入模型的變量。有的參數(shù)可能對密度的影響不是很大,而且各參數(shù)之間不完全相互獨立,可能有互相作用的關(guān)系。因此,這里利用逐步回歸方法篩選并剔除引起多重共線性的變量?；舅枷胧菍⒆兞恐饌€引入模型,并對引入的變量進行F檢驗,對已經(jīng)選入的變量進行t檢驗,當原來引入的變量由于后面變量的引入變得不再顯著時,則將其刪除。如此反復(fù),直到既沒有顯著的變量被選入,也沒有不顯著的變量被剔除為止。這樣建立的回歸模型預(yù)測效果會更好。同時嘗試用篩選出的參數(shù)建立BP神經(jīng)網(wǎng)絡(luò)模型。為了驗證模型的可靠性,采用平均絕對誤差(MAE)、平均絕對百分比誤差(MAPE)、均方誤差(MSE)、均方根誤差(RMSE)、相關(guān)系數(shù)(R)、決定系數(shù)(R2)評價。其中,平均絕對誤差、平均絕對百分比誤差、均方誤差和均方根誤差可以從整體上評價模型預(yù)測的誤差大小,數(shù)值越小說明模型預(yù)測的精度越高。相關(guān)系數(shù)和決定系數(shù)越接近于1,表明所建立的模型效果越好。

2 結(jié)果與討論

2.1 多元線性回歸模型的建立

多元線性回歸模型考慮的是多個因素對同一個因素的影響,即多個自變量與因變量的關(guān)系,可表示為:

yi=b0+b1xi1+b2xi2+…+bpxip

(1)

式中:yi為預(yù)測值,x1、x2……xp為p個自變量,b0為常數(shù)項,b1、b2……bp為回歸系數(shù),采用矩陣形式求出每個回歸系數(shù)。

為了驗證所建立的模型有良好的外部預(yù)測能力,將112個樣本數(shù)據(jù)集隨機分為訓(xùn)練集92個,用于建立模型。剩余20個化合物作為測試集,驗證該模型的預(yù)測能力。首先選用B3PW91/6-31G(d,p)下計算所得參數(shù),采用SPSS軟件進行多元逐步回歸,結(jié)果見表1。

表1 多元逐步回歸結(jié)果Table 1 Results of multiple stepwise regression

表2 模型系數(shù)Table 2 Parameters for models

分別將訓(xùn)練集與測試集中的化合物各參數(shù)代入上述模型,計算得到的預(yù)測值與實測值作對比,結(jié)果見圖1和圖2。由圖1可知,3種計算理論水平下,數(shù)據(jù)較集中的分布在對角線附近兩側(cè),說明該模型有一定的預(yù)測能力。由圖2可知,隨著預(yù)測值的變化,殘差均在一定的水平帶狀范圍內(nèi)波動,說明模型效果較好,由于少數(shù)誤差較大的點使得帶狀區(qū)域較寬,預(yù)測精度不是很高。經(jīng)比較,在B3PW91泛函下計算所得參數(shù)建立的MLR模型對于訓(xùn)練集和測試集的平均絕對誤差、平均絕對百分比誤差、均方誤差均低于B3LYP泛函和M06-2X泛函計算結(jié)果,而相關(guān)系數(shù)和決定系數(shù)高于B3LYP泛函和M06-2X泛函。M06-2X泛函計算結(jié)果稍好于B3LYP泛函,具體各評價指標見表3。但是從整體分析得出,3種理論水平下建立的MLR預(yù)測模型精度都不高,預(yù)測結(jié)果不是很好,由此考慮采用其他方法建立模型。因BP神經(jīng)網(wǎng)絡(luò)能學習和存儲大量的輸入與輸出之間的映射關(guān)系,無需事前提供描述這種映射關(guān)系的數(shù)學方程[13]。因此,嘗試采用相同的訓(xùn)練集和測試集,建立BP神經(jīng)網(wǎng)絡(luò)模型預(yù)測雜環(huán)含能材料的晶體密度。

圖1 MLR模型預(yù)測值與實測值比較Fig.1 Comparison of predicted value obtained by MLR model and experimental value

圖2 MLR模型預(yù)測值與殘差的關(guān)系Fig.2 The relationship between predicted value obtained by MLR model and residual

表3 MLR預(yù)測模型的評價指標Table 3 Evaluation metrics for MLR model

2.2 BP神經(jīng)網(wǎng)絡(luò)模型的建立

由BP神經(jīng)網(wǎng)絡(luò)模型的拓撲結(jié)構(gòu)可知,BP神經(jīng)網(wǎng)絡(luò)的構(gòu)建需確定輸入層、隱含層和輸出層。由逐步回歸方法篩選出的參數(shù)可知,有4個輸入?yún)?shù)和1個輸出參數(shù),但隱含層數(shù)和節(jié)點數(shù)是不確定的。為了提高網(wǎng)絡(luò)訓(xùn)練精度,一般采取增加節(jié)點數(shù)和隱含層數(shù)的方法。但是在網(wǎng)絡(luò)結(jié)構(gòu)實現(xiàn)上,增加節(jié)點數(shù)要比增加隱含層數(shù)簡單,因此選擇單隱含層。而隱含層節(jié)點數(shù)的確定,目前在理論上還沒有作出明確的規(guī)定。一般根據(jù)經(jīng)驗公式,在一定范圍內(nèi),采用嘗試法,通過對不同節(jié)點數(shù)進行訓(xùn)練對比,選取最合適的網(wǎng)絡(luò)[25]。

BP神經(jīng)網(wǎng)絡(luò)不同的傳遞函數(shù)、訓(xùn)練函數(shù)會影響網(wǎng)絡(luò)預(yù)測結(jié)果。具體參數(shù)設(shè)定為:學習速率0.01,訓(xùn)練次數(shù)1000次,訓(xùn)練目標最小誤差0.0001。在MATLAB中編寫程序完成算法過程。分別在不同的傳遞函數(shù)、訓(xùn)練函數(shù)和隱含層節(jié)點個數(shù)下進行訓(xùn)練,經(jīng)過多次訓(xùn)練得到測試集預(yù)測的精度最好時的網(wǎng)絡(luò)結(jié)構(gòu)為4-7-1。

BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)測值與實測值比較結(jié)果如圖3和圖4所示;BP預(yù)測模型的評價指標見表4。由圖3和圖4可以看出,對于訓(xùn)練集和測試集,數(shù)據(jù)整體上較為緊密地分布在對角線附近,殘差點也分布在較窄的帶狀區(qū)域內(nèi),且B3PW91泛函計算結(jié)果優(yōu)于B3LYP泛函和M06-2X泛函,B3LYP泛函和M06-2X泛函計算結(jié)果相差不大。同時由表4可知,和多元回歸模型比較,平均絕對誤差、平均絕對百分比誤差、均方誤差均有所降低,而相關(guān)系數(shù)、決定系數(shù)均有所提高,說明BP神經(jīng)網(wǎng)絡(luò)模型較多元回歸模型的預(yù)測效果好,預(yù)測精度高。

圖3 BP模型預(yù)測值與實測值比較Fig.3 Comparison of predicted value obtained by BP model and experimental value

圖4 BP模型預(yù)測值與殘差的關(guān)系Fig.4 The relationship between predicted value obtained by BP model and residual

表4 BP預(yù)測模型的評價指標Table 4 Evaluation metrics for BP model

2.3 預(yù)測模型評估與比較

為了評估模型的可靠性和穩(wěn)定性,除了應(yīng)用測試集對模型進行驗證外,同時采用留出法將樣本數(shù)據(jù)進行多次建模后取平均值作為評估結(jié)果,并與Politzer等[6]和Rice等[7]提出的預(yù)測公式進行比較,相關(guān)評價指標見表5。

表5 所建模型與現(xiàn)有模型的比較Table 5 Comparison of the established models and the existing models

由表5評估結(jié)果可知,3種泛函下計算所得參數(shù)建立的模型,B3PW91泛函結(jié)果優(yōu)于其他兩種。在B3LYP/6-31G(d,p)水平下,平均絕對誤差、均方誤差、平均絕對百分比誤差,MLR模型預(yù)測結(jié)果較Rice等[7]的計算結(jié)果分別降低了0.003g/cm3、0.002(g/cm3)2、0.002,而BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型分別降低了0.004g/cm3、0.004(g/cm3)2、0.002,結(jié)果較好于MLR預(yù)測結(jié)果。同樣,在B3PW91/6-31G(d,p)水平下,MLR模型預(yù)測結(jié)果較Politzer等[6]的計算結(jié)果分別降低了0.01g/cm3、0.001(g/cm3)2、0.006,而BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型分別降低了0.028g/cm3、0.006(g/cm3)2、0.018,結(jié)果仍然優(yōu)于MLR預(yù)測結(jié)果。總體可看出,對于氮/氧雜環(huán)含能化合物在原有分子靜電勢法的基礎(chǔ)上添加了氫鍵供體數(shù)和氫鍵受體數(shù)兩個參數(shù)后,預(yù)測精度提高,且在B3PW91/6-31G(d,p)水平下,BP神經(jīng)網(wǎng)絡(luò)建立的密度預(yù)測模型精度最高,但是沒有達到最理想的效果,可能與量子化學計算方法所得參數(shù)的精度及樣本數(shù)量有關(guān)。BP神經(jīng)網(wǎng)絡(luò)由于具有較強的映射能力,能夠通過學習自動提取輸入、輸出數(shù)據(jù)間的“合理規(guī)則”,通過反向傳播不斷調(diào)整網(wǎng)絡(luò)的權(quán)值和閾值,使網(wǎng)絡(luò)的誤差平方和最小,得到比多元回歸模型較好的預(yù)測結(jié)果。

3 結(jié) 論

(1)運用量子化學計算方法,分別選擇B3LYP泛函、B3PW91泛函和M06-2X泛函,在6-31G(d,p)基組水平下對112種氮、氧雜環(huán)含能材料的分子結(jié)構(gòu)進行優(yōu)化,分析其振動頻率后得到最穩(wěn)定分子結(jié)構(gòu),檢索出7類量化結(jié)構(gòu)參數(shù),并從PubChem數(shù)據(jù)庫查詢到3類結(jié)構(gòu)參數(shù),為建立模型統(tǒng)一了結(jié)構(gòu)參數(shù)的來源。

(2)通過逐步多元回歸方法篩選出對密度影響較大的參數(shù)并分別建立回歸模型和BP神經(jīng)網(wǎng)絡(luò)模型,結(jié)果表明,在3種泛函計算水平下,均為BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)測誤差較小、精度較高。

(3)添加氫鍵供體數(shù)和氫鍵受體數(shù)兩個參數(shù)后,與已有模型比較,MLR和BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)測精度均有所提高,且在B3PW91/6-31G(d,p)水平下,BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)測精度最高。