文琳元，何曉凱，師進(jìn)文，劉英哲，王伯周

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

引 言

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

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

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

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

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

1 流程及原理

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

根據(jù)分子表面靜電勢(shì)理論，利用Multiwfn[21]計(jì)算得到各參數(shù)，表1列出幾個(gè)典型靜電勢(shì)參數(shù)。基于多元線性擬合方法，對(duì)Politzer所提的密度擬合公式進(jìn)行修正，并在靜電勢(shì)參數(shù)基礎(chǔ)上采用符號(hào)回歸方法得到最優(yōu)的公式形式。

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

Continued

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

(1)

(2)

(3)

2 結(jié)果與討論

根據(jù)公式(1)，其實(shí)驗(yàn)晶體密度與預(yù)測(cè)密度關(guān)系如圖1所示。

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

圖公式下預(yù)測(cè)密度與實(shí)驗(yàn)晶體密度關(guān)系Fig.2 The relationship between predicted density and experimental crystal density under formula

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

圖3 Politer-Π公式下預(yù)測(cè)密度與實(shí)驗(yàn)晶體密度關(guān)系Fig.3 The relationship between predicted density and experimental crystal density under Politer-Π formula

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

(4)

(5)

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

圖4 多元線性擬合改進(jìn)后密度預(yù)測(cè)公式中預(yù)測(cè)密度與實(shí)驗(yàn)晶體密度關(guān)系Fig.4 The relationship between predicted density and experimental crystal density in the improved density prediction formula by multiple linear regression

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

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

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

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

(6)

(7)

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

圖6 符號(hào)回歸方法改進(jìn)后密度預(yù)測(cè)公式中預(yù)測(cè)密度與實(shí)驗(yàn)晶體密度關(guān)系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 各公式預(yù)測(cè)精度對(duì)比Table 2 Comparison of the prediction accuracy of different formulas

3 結(jié) 論

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

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

(3)符號(hào)回歸構(gòu)建的預(yù)測(cè)公式提高了新型含氟高能致爆材料晶體密度預(yù)測(cè)精度，有助于加速含氟高能致爆材料的設(shè)計(jì)。與此同時(shí)，符號(hào)回歸方法還可以運(yùn)用于其他研究體系，以破譯材料性能間高維復(fù)雜的關(guān)系，構(gòu)建出可理解的數(shù)學(xué)表達(dá)式，為材料設(shè)計(jì)提供思路。