文琳元，何曉凱，師進文，劉英哲，王伯周

(1.西安交通大學 能源與動力工程學院，陜西 西安 710049；2.西安近代化學研究所，陜西 西安 710065)

引 言

作為炸藥、發(fā)射藥、推進劑和火工品中高能量組分的含能材料，其能量密度水平一定程度上決定了武器發(fā)射、推進及毀傷能力，乃至國家整體軍事裝備的戰(zhàn)略威懾和打擊能力。追求更高能量密度，成為了含能材料領域科學家為之長期奮斗的目標。為此，研究者們提出了引入高能致爆基團構筑新型高能量密度含能材料的策略[1]。同時，由于新材料的合成通常需要大量的時間與金錢，因此準確預測含能材料性能對于開發(fā)具有優(yōu)良性能的新型含能材料至關重要。隨著計算科學與理論化學的快速發(fā)展，新型含能材料的理論設計技術日漸成熟，具體為通過量化計算預估新化合物的生成焓、密度、爆轟速度、爆轟壓力等參數(shù)，輔以性能指標要求篩選到最終的目標化合物。其中，因為K-J方程中爆轟速度與密度正相關且爆轟壓力與密度的平方成正比，所以密度被認為是影響含能材料爆轟性能的關鍵參數(shù)，如何準確預測新型含能材料晶體密度也就成為新型含能材料研發(fā)中一大關鍵問題。

當前含能材料晶體密度預測方法包含物理模型與統(tǒng)計模型兩種。其中物理模型可細分為基于分子體積[2-4]以及基于晶體體積[5-9]兩類，而基于統(tǒng)計模型中較為主流的方法則是機器學習[10]。目前基于分子體積方法用的較多的是分子表面靜電勢法，該方法是用分子摩爾質量(M)除以使用Monte-Carlo方法計算得到的分子0.001a.u.電子密度等值面包裹內(nèi)體積V0.001的值來表示晶體密度[11-12]，該方法未考慮分子間作用和晶體堆積，計算誤差與速度適中。Politzer 在分子表面靜電勢法的基礎上，針對CHON類材料提出了分子表面靜電勢參數(shù)修正的方法[13-14]，其預測密度與實驗值誤差小于0.05g/cm3，被廣泛運用于含能材料晶體密度預測。

基于晶體體積方法則是通過晶體結構預測，以獲得含能材料準確的晶體結構，從而得到該材料的晶體密度。該方法需要準確描述分子間相互作用，存在計算速度較慢的不足，是當前含能材料晶體密度研究的一大挑戰(zhàn)。而機器學習模型主要是通過各種回歸算法構建材料結構與密度之間的關系，其中大多數(shù)是黑箱模型，缺乏解釋性。近年來，不少研究者基于該方法在CHON類含能材料晶體密度預測上做了許多優(yōu)秀工作，其中張朝陽研究員針對超過2000個硝基化合物構建了均方根誤差(RMSE)=0.049g/cm3的圖神經(jīng)網(wǎng)絡密度預測模型[15]，T. Yong-Jin Han針對10251個含有至少一根氮氧鍵的CHON類化合物開發(fā)了RMSE=0.062g/cm3消息傳遞神經(jīng)網(wǎng)絡模型[16]，張慶華研究員對超過1000個CHON類含能材料利用核嶺回歸的方法構建了MAE=0.042g/cm3的密度預測模型[17]。

從上述研究可以發(fā)現(xiàn)，其樣本量均大于103，這是因為機器學習模型的計算精度十分依賴數(shù)據(jù)樣本，樣本量過少會使得模型過擬合。這其中，符號回歸作為一種回歸分析方法，在挖掘變量間內(nèi)在關系以及幫助理解物理機制并發(fā)現(xiàn)隱藏的數(shù)學表達式方面取得一定的成果[18-19]。

由于當前晶體密度預測方法大多針對CHON類含能材料，沒有考慮氟代偕二硝基(—CF(NO2)2)、二氟氨基(—C(NF2))以及氟代偕硝基(—CF(NO2))等含氟高能致爆基團對密度預測的影響，難以保證上述方法在預測含氟高能致爆材料晶體密度的精度。因此，有必要建立一個適用于含氟高能致爆材料的晶體密度預估方法。其中，基于晶體體積方法存在計算速度較慢、缺乏準確描述分子間相互作用的不足，因此本研究僅采用基于分子體積與基于機器學習的方法。先在Politzer基于CHON類含能材料密度預測公式基礎上，利用多元線性擬合方法得到適合含氟高能致爆材料體系的晶體密度預測公式?？紤]到含有上述含氟高能致爆基團材料數(shù)量較少的情況，本研究將機器學習模型中的符號回歸算法與基于分子體積方法相結合，從多種靜電勢參數(shù)中確定出更適合于該研究體系的變量表現(xiàn)形式，構造出了適用于含氟高能致爆材料體系的R2=0.77、RMSE=0.045g/cm3的晶體密度預測公式。所得公式能提高含氟高能致爆材料晶體密度預測精度，以期為含氟高能致爆材料的高效準確設計奠定基礎。

1 流程及原理

從劍橋晶體結構數(shù)據(jù)庫(CCDC)中，按照實驗密度大于1.6g/cm3、晶體結構為單晶、分子中存在含氟高能致爆基團等標準，篩選出56個化合物。再利用Gaussian16[20]軟件的B3PW91/6-31g(d,p)方法，對上述化合物的分子進行幾何結構優(yōu)化并得到能量最優(yōu)的構型，經(jīng)頻率分析確定無虛頻。

根據(jù)分子表面靜電勢理論，利用Multiwfn[21]計算得到各參數(shù)，表1列出幾個典型靜電勢參數(shù)?；诙嘣€性擬合方法，對Politzer所提的密度擬合公式進行修正，并在靜電勢參數(shù)基礎上采用符號回歸方法得到最優(yōu)的公式形式。

表1 56個含氟高能致爆材料的部分靜電勢參數(shù)及實驗晶體密度Table 1 Electrostatic potential parameters and experimental crystal densities for 56 fluorine-containing high energetic explosophoric materials

Continued

其中，選取了Rice在2007年所提的公式(1)[12]、Politzer在2009年所提的兩個公式(2)和(3)[13]作為對比，并計算了各公式在預測含氟高能致爆材料時的均方根誤差(RMSE)。

(1)

(2)

(3)

2 結果與討論

根據(jù)公式(1)，其實驗晶體密度與預測密度關系如圖1所示。

從圖1可以清楚地發(fā)現(xiàn)，其預測密度幾乎都大于對應的實驗晶體密度，這首先說明氟元素加入改變了其分子間相互作用強弱，使得0.001a.u.電子密度等值面包裹的體積Vm小于等效的晶體體積；同時，按照Goh的標準[22]，公式(1)計算中僅有17.9%的化合物預測密度可視為準確(誤差小于 0.03g/cm3)，25.0%的化合物預測密度可視為有用(誤差小于0.05g/cm3)，而公式(1)預測密度的均方根偏差高達0.106g/cm3，不足以達到晶體密度預測的要求。

圖公式下預測密度與實驗晶體密度關系Fig.2 The relationship between predicted density and experimental crystal density under formula

從圖1和圖2可以看出，公式(2)得到預測密度比公式(1)預測密度小，這主要在于公式(2)中主要影響因素(M/Vm)前面的系數(shù)α小于公式(1)中的系數(shù)α=1。該公式使得預測結果更偏向于實際晶體密度，但從圖2看，其仍呈現(xiàn)出高估了晶體密度的現(xiàn)象。相較沒考慮分子間相互作用的公式(1)，公式(2)計算得到的晶體密度與實際晶體密度的偏差減小，均方根偏差從0.106g/cm3降到0.075g/cm3，在56種含氟高能致爆材料中，預測晶體密度可視為準確的比例提高到26.8%，而視為有用的比例提高到44.6%，幾乎是成倍地增長。這說明，加入靜電勢修正項，是符合物理本質且對于含氟晶體密度預測有用的。但是，其中仍有21.4%的結果偏差是大于0.100g/cm3，說明這一套僅考慮C、H、O、N元素的靜電勢修正擬合參數(shù)，不足以滿足現(xiàn)有的需求。

圖3 Politer-Π公式下預測密度與實驗晶體密度關系Fig.3 The relationship between predicted density and experimental crystal density under Politer-Π formula

由于上述使用Politer所提的分子表面靜電勢修正公式在預測含氟高能致爆材料晶體密度時誤差仍較大，因此為了更為準確地描述氟元素加入后分子間相互作用對預測公式的影響。綜合考慮到(M/Vm)在公式(1)、(2)、(3)中關鍵作用，在公式(2)、(3)的基礎采用多元線性擬合方法修正密度預測公式,得到如下式(4)、(5)：

(4)

(5)

發(fā)現(xiàn)修正公式(4)和(5)的決定系數(shù)與均方根誤差都為0.65和0.055g/cm3；而公式(4)在密度預測誤差小于0.05g/cm3和小于0.03g/cm3上的表現(xiàn)優(yōu)于公式(5)，因此在后續(xù)僅討論公式(4)的效果。其預測密度與實際晶體密度如圖4所示，相較于前面提及的3種方法，改進后預測結果的相對誤差都在5%以內(nèi)，而且均方根偏差也降至0.055g/cm3。預測晶體密度可視為準確的比例為33.9%，而視為有用的比例提高到60.7%。不僅如此，預測晶體密度的偏差大于0.100g/cm3的比例也降到1.8%。

圖4 多元線性擬合改進后密度預測公式中預測密度與實驗晶體密度關系Fig.4 The relationship between predicted density and experimental crystal density in the improved density prediction formula by multiple linear regression

雖然公式(4)中靜電勢方差參數(shù)在含氟高能致爆材料晶體密度預估公式中作用不明顯，但將公式(4)中靜電勢方差替換為其他靜電勢參數(shù)的組合，可能會更好地體現(xiàn)含氟高能致爆材料體系內(nèi)相互作用。由此課題組基于符號回歸概念，開展含氟高能致爆材料的晶體密度預測研究。在本研究中，考慮到公式(2)以及(3)中靜電勢參數(shù)主要的運算方式為“加、減、乘、除”(分別對應于“add、sub、mul、div”等4個算符)，

因此，在符號回歸中基于上述4種運算符構建模型。如圖 5所示，利用符號回歸模型后得到的含氟高能致爆材料晶體密度預測公式的二叉樹表示，其中每個葉節(jié)點都是變量或者常數(shù)，而內(nèi)部節(jié)點是選定的運算符。

圖5 含氟高能致爆材料晶體密度預測公式-二叉樹表示Fig.5 Density prediction formula of fluorine-containing high energetic explosophoric materials-binary tree representation

根據(jù)圖5的二叉樹，可以得到其表達式，見公式(6)：

(6)

(7)

符號回歸得到的密度預測公式，其R2=0.77，RMSE=0.045g/cm3，優(yōu)于利用多元線性擬合方法修正的公式(4)。預測密度與實驗晶體密度的分布如圖 6所示，可以發(fā)現(xiàn)該公式彌補之前數(shù)個公式在預測實驗密度低于1.75g/cm3的材料時高估密度以及預測實驗密度高于2.00g/cm3的材料時低估密度的問題，表現(xiàn)在圖中就是頭尾的數(shù)據(jù)更加趨近于藍色點線。如表2所示，相比于多元線性擬合方法，符號回歸在降低預測均方根誤差的同時，還提高了預測公式的綜合表現(xiàn)。具體而言，一半以上材料的密度預測誤差控制在0.03g/cm3以下，超過70%材料的預測誤差低于0.05g/cm3。結合之前對于公式(2)、(3)和(4)的分析，公式(7)等號右側第二項是由靜電勢方差、靜電勢偏差、0.001a.u電子密度等值面面積以及分子摩爾質量組合的新變量，該參數(shù)大致體現(xiàn)含氟高能致爆材料復雜的分子間相互作用并以數(shù)學表達式呈現(xiàn)，使得該公式能夠彌補之前談到的M/Vm對密度的高估問題，從而也就解釋了為何該公式在密度預測時較優(yōu)的表現(xiàn)。

圖6 符號回歸方法改進后密度預測公式中預測密度與實驗晶體密度關系Fig.6 The relationship between the predicted density and the experimental crystal density in the density prediction formula after the improvement of the symbolic regression method

表2 各公式預測精度對比Table 2 Comparison of the prediction accuracy of different formulas

3 結 論

(1)基于分子體積方法，在Politzer的CHON晶體密度預測公式基礎上，運用多元線性擬合得到了均方根誤差為0.055g/cm3的修正公式。

(2)結合符號回歸方法構建了適合含氟高能致爆材料晶體密度預測的公式，其均方根誤差為0.045g/cm3，極大地提高了預測精度，其中有超過50%的化合物預測誤差低于0.03g/cm3。

(3)符號回歸構建的預測公式提高了新型含氟高能致爆材料晶體密度預測精度，有助于加速含氟高能致爆材料的設計。與此同時，符號回歸方法還可以運用于其他研究體系，以破譯材料性能間高維復雜的關系，構建出可理解的數(shù)學表達式，為材料設計提供思路。