氮/氧雜環含能材料晶體密度的預測研究

呂志濤,劉玉存,聶江穩,賈康輝,李 強

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

引 言

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

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

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

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

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

1 實 驗

1.1 樣本數據計算與收集

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

1.2 研究方法

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

2 結果與討論

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

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

yi=b0+b1xi1+b2xi2+…+bpxip

(1)

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

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

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

表2 模型系數Table 2 Parameters for models

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

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

圖2 MLR模型預測值與殘差的關系Fig.2 The relationship between predicted value obtained by MLR model and residual

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

2.2 BP神經網絡模型的建立

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

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

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

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

圖4 BP模型預測值與殘差的關系Fig.4 The relationship between predicted value obtained by BP model and residual

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

2.3 預測模型評估與比較

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

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

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

3 結 論

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

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

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