常曉雅,文明杰,張 迪,王永錦,初慶釗,朱 通,陳東平

(1.北京理工大學(xué) 爆炸科學(xué)與技術(shù)重點實驗室,北京 100081;2.上海市分子治療與新藥開發(fā)工程研究中心, 華東師范大學(xué) 化學(xué)與分子工程學(xué)院,上海 200062)

引 言

含能材料通常指在受到外界刺激時能快速分解并釋放巨大能量的化合物或混合物[1-2],被廣泛應(yīng)用于軍事燃燒及爆破、火箭導(dǎo)彈發(fā)射、推進劑等領(lǐng)域。在推進劑領(lǐng)域,含能材料以燃燒的形式釋放能量,實現(xiàn)了從化學(xué)能向動能的有效轉(zhuǎn)化,成為先進導(dǎo)彈武器發(fā)動機的動力來源,是推動武器裝備系統(tǒng)發(fā)展的支撐和制約技術(shù)之一[3-4];在火炸藥領(lǐng)域,含能材料在極短的時間內(nèi)釋放大量的熱能和氣體,以爆轟形式作功形成沖擊波,是武器裝備完成發(fā)射和高效毀傷的能源材料。作為現(xiàn)代武器系統(tǒng)的主要能量來源,含能材料燃燒與爆炸性能直接影響武器裝備的打擊范圍和毀傷能力,是先進武器裝備優(yōu)化設(shè)計的關(guān)鍵發(fā)展方向。

含能材料釋放的能量主要來源于原子核外層電子在反應(yīng)過程中轉(zhuǎn)移所釋放的能量[5],其宏觀特性本質(zhì)上是由其微觀分子結(jié)構(gòu)所決定的。例如,含能材料爆轟反應(yīng)放熱量取決于產(chǎn)物分子的生成熱;含能材料的感度與分子結(jié)構(gòu)、分子間相互作用以及反應(yīng)活化能等性質(zhì)有關(guān);各組分之間的相容性與分子間的相互作用有關(guān)[6-7]。因此,在微觀結(jié)構(gòu)層面上研究含能分子結(jié)構(gòu)燃燒特性和力學(xué)穩(wěn)定性對于新型含能材料的設(shè)計至關(guān)重要。分子模擬作為原子尺度的研究方法,可以用來預(yù)測含能材料的力學(xué)行為和化學(xué)反應(yīng)過程,揭示分子結(jié)構(gòu)與宏觀燃燒性能之間的關(guān)系,從而為含能材料的理論設(shè)計提供科學(xué)依據(jù)。在分子模擬計算中,描述微觀體系中不同分子構(gòu)型與能量的映射關(guān)系被稱為勢能面。分子模擬的實質(zhì)是在被模擬材料分子體系的勢能面上進行動力學(xué)研究,例如搜尋分子的最低能量構(gòu)形或者計算化學(xué)反應(yīng)速率。其中構(gòu)建勢能面的計算方法主要包括第一性原理計算和分子動力學(xué)。第一性原理計算是通過求解薛定諤方程來描述粒子的電子結(jié)構(gòu),在不需要任何先驗參數(shù)情況下,準確計算分子、原子的理化性質(zhì),計算精度較高。在空間尺度較小且精度要求高的場景,比如少量分子體系的化學(xué)反應(yīng)路徑研究、動力學(xué)參數(shù)計算等,第一性原理計算表現(xiàn)出獨特的優(yōu)勢。然而高精度的計算水平伴隨著高昂的計算成本,時間復(fù)雜度往往是O(N3)甚至更高,現(xiàn)有的計算能力無法承載超過1 000個原子的分子體系。

分子動力學(xué)不考慮電子運動的自由度,電子效應(yīng)體現(xiàn)在原子間相互作用的勢函數(shù)中,原子的運動遵循牛頓第二定律,并據(jù)此實現(xiàn)目標體系在給定初始條件下隨時間的動態(tài)演化。分子模擬的精度取決于勢函數(shù)對真實勢能面重建的精度。勢函數(shù)的數(shù)學(xué)形式和分子各部分對能量貢獻的可拓展性等均對勢能面重建精度有影響[8]。非反應(yīng)力場假設(shè)的函數(shù)形式簡單,如施加兩體或三體截斷,對于材料的力學(xué)性質(zhì)描述較為準確,其研究尺度可達上百萬甚至千萬原子的體系[9]。但這樣的經(jīng)驗力場并不適用于研究電子轉(zhuǎn)移、原子變價的化學(xué)反應(yīng)過程,無法用于模擬含能材料燃燒與爆炸等復(fù)雜反應(yīng)過程。反應(yīng)力場是基于鍵級來描述體系的鍵合作用,鍵級隨著化學(xué)反應(yīng)的發(fā)生自發(fā)進行調(diào)整[10],包括Airebo[11]、COMB3[12]和ReaxFF等力場[13-14]。其中,應(yīng)用最廣泛的是Van Duin等[13-14]提出的ReaxFF力場,通過對鍵能、范德華能等多個能量項進行參數(shù)擬合來構(gòu)建勢能面,目前在含能材料點火、燃燒、沖擊等方面均有大量研究[15-20]。利用反應(yīng)力場進行分子模擬,結(jié)合模擬軌跡分析,可得到反應(yīng)中間體、反應(yīng)產(chǎn)物等隨時間的變化,從而提取其主要的基元反應(yīng)、構(gòu)建詳細反應(yīng)機理。精確的力場模型是探究含能材料反應(yīng)機理的基礎(chǔ)。此外,并不是所有含能材料體系均已開發(fā)ReaxFF力場,以AP為例,盡管已經(jīng)廣泛驗證其作為氧化劑可以增加含能材料的能量釋放和反應(yīng)活性,但是并未有前人開發(fā)相應(yīng)的ReaxFF力場,因此目前并沒有關(guān)于AP燃燒過程的分子模擬研究。

第一性原理計算精度高,但計算成本大;分子動力學(xué)模擬效率較高,但是精度較低。兩種方法均無法兼顧精度與效率。如何高效率、高精度地重建勢能面是研究人員面臨已久的難題。近年來,基于機器學(xué)習的勢函數(shù)模型有望解決傳統(tǒng)計算方法的“卡脖子”問題。該方法繞過理論模型精確解釋的難題,通過完全數(shù)值化對抽象的物理規(guī)律和復(fù)雜的化學(xué)過程之間的潛在關(guān)系進行推斷,從而準確高效地構(gòu)建反應(yīng)體系的勢能面,實現(xiàn)多種復(fù)合含能體系的理論計算,建立分子結(jié)構(gòu)與燃燒性能的聯(lián)系。

本文歸納了機器學(xué)習勢函數(shù)模型搭建方案及訓(xùn)練集搭建策略,通過對高精度的量子化學(xué)勢能面進行近似,成功應(yīng)用于含能材料燃燒爆炸方面的研究,包括常用的RDX、高氮炸藥CL-20、新型硝胺炸藥ICM-102、強氧化劑AP、高能顆粒(鋁和硼)等物質(zhì),并綜述了機器學(xué)習勢在碳氫燃料燃燒中的研究進展。最后總結(jié)了機器學(xué)習精度和效率方面的優(yōu)勢,并對其發(fā)展趨勢提出展望。

1 機器學(xué)習勢

1.1 機器學(xué)習勢模型

精確構(gòu)造含能材料的勢能面是開發(fā)機器學(xué)習勢的基本要求,也是研究其微觀燃燒機制的基礎(chǔ)。機器學(xué)習勢從第一性原理數(shù)據(jù)出發(fā),用高維函數(shù)來描述分子結(jié)構(gòu)與能量之間的關(guān)系,進而擬合勢能面,以較低的計算成本實現(xiàn)對目標體系的分子模擬計算。一個合格的機器學(xué)習勢模型需要至少滿足廣延性、對稱不變性和光滑性[9]。體系的廣延性是指體系總能量和原子數(shù)成正比,故機器學(xué)習勢的能量分解形式為:

(1)

(2)

為了達到上述要求,通常引入描述符D(Ri)的概念對原子三維坐標進行映射。對于進行等價對稱操作的原子坐標,其描述符相同,即:

D(Ri)=D(URi),?U∈u

(3)

從而將原子能量貢獻改寫為:

(4)

其中F為擬合方法。將原子坐標映射得到的結(jié)構(gòu)描述符作為輸入,通過不同的擬合方法得到分子結(jié)構(gòu)和能量的對應(yīng)關(guān)系,即機器學(xué)習勢模型。模型的光滑性由擬合方法的光滑性和描述符共同保證。本節(jié)簡要介紹機器學(xué)習勢的擬合方法和描述符構(gòu)造方法。

目前主流的機器學(xué)習擬合方法包括核方法和神經(jīng)網(wǎng)絡(luò)方法。其中,核方法通過測量描述符與數(shù)據(jù)集中描述符的相似程度來獲得系統(tǒng)的能量或力。2010年,Csanyi等[21]基于核方法提出高斯近似勢(Gaussian Approximation Potential, GAP)模型,通過對原子化學(xué)環(huán)境相似性進行度量,采用高斯回歸過程對勢能面進行內(nèi)插,從而得到原子能量、受力的最優(yōu)估值[22]。分子構(gòu)型越相似,在勢能面上的位置就越接近,構(gòu)型的勢能也應(yīng)當越相近[23]。這種方法并未預(yù)設(shè)函數(shù)形式,訓(xùn)練過程人工干預(yù)較少,且可解釋性較強。目前在材料領(lǐng)域已經(jīng)發(fā)表了大量的相關(guān)研究工作,如鈣鈦礦[24]、過渡態(tài)金屬[25-26]、硅烯[27]、半導(dǎo)體[28]、合金[29]等。2018年以來,基于核方法的(對稱)梯度域機器學(xué)習((Symmetric)Gradient Domain Machine Learning, GDML)[30]方法快速發(fā)展。與高斯近似勢不同,GDML模型在初始數(shù)據(jù)集中僅加入第一性原理計算的原子受力,通過對原子受力進行積分得到體系總能量[31]。引入受力項的模型比基于能量的模型自由度更高、參數(shù)更多、形式更靈活。

神經(jīng)網(wǎng)絡(luò)勢函數(shù)(Neural Network Potential, NNP)模型是用一個多參數(shù)的非線性函數(shù)來擬合勢能面。早在1995年,Blank等[32]利用前饋神經(jīng)網(wǎng)絡(luò)擬合勢能面,研究CO在Ni表面以及H2在Si表面的擴散問題。該模型的輸入坐標并不滿足平移、旋轉(zhuǎn)和交換對稱不變性,故前饋神經(jīng)網(wǎng)絡(luò)僅僅應(yīng)用于低維體系。2007年,Behler和Parrinello等[33]提出了高維神經(jīng)網(wǎng)絡(luò)(High-Dimensional Neural Network, HDNN)模型, 采用以原子為中心的對稱函數(shù)矢量來描述原子的鄰域環(huán)境,每個原子用一個單獨的前饋神經(jīng)網(wǎng)絡(luò)表示,從而預(yù)測每個原子對體系全局能量的貢獻。這樣的神經(jīng)網(wǎng)絡(luò)預(yù)測能力更強,但需要反復(fù)調(diào)節(jié)和擬合大量的網(wǎng)絡(luò)參數(shù), 且需要充足的訓(xùn)練樣本。一般來說,神經(jīng)網(wǎng)絡(luò)模型的深度比其寬度更受青睞,即隱藏層的層數(shù)比每層的節(jié)點數(shù)更重要[34]。這是因為更深的神經(jīng)網(wǎng)絡(luò)具有更好的非線性表達能力,可學(xué)習更復(fù)雜的變化,從而擬合更復(fù)雜的特征。換言之,包含多個連續(xù)隱藏層的深度神經(jīng)網(wǎng)絡(luò)的高維非線性函數(shù)擬合能力更為強大。2018年,鄂維南、張林峰等[35-36]基于深度神經(jīng)網(wǎng)絡(luò)提出了深度勢(Deep Potential, DP)模型。同年,該團隊開發(fā)了開源工具DeePMD-kit[37]來支持DP模型的開發(fā)工作。目前DP模型已被用于材料、力學(xué)、化學(xué)、生物等各個領(lǐng)域的研究,如電池[38-39]、催化[40-41]、蛋白質(zhì)[42-43]、含能材料[44-46]等。以DP模型為代表的機器學(xué)習勢模型,通過機器學(xué)習算法加速和大規(guī)模并行優(yōu)化,可在保持第一性原理計算精度的同時實現(xiàn)上億原子的分子動力學(xué)模擬,將精確的物理建模帶入了更大尺度的材料分子模擬[47]。

機器學(xué)習勢通過構(gòu)造結(jié)構(gòu)和能量之間的映射關(guān)系來重建勢能面,因此采用何種描述符來提取分子結(jié)構(gòu)特征進而用于訓(xùn)練機器學(xué)習模型對于力場開發(fā)的精度和效率至關(guān)重要。Behler和Parrinello等[33]構(gòu)造的原子中心對稱函數(shù)(Atom-Centered Symmetry Function, ACSF)滿足物理對稱性以及相對于移動原子的可微性,被廣泛集成于多個軟件中[48-49]。該方法以原子為中心,用徑向函數(shù)和角向函數(shù)將鄰域原子的化學(xué)信息編碼為復(fù)雜的“多體”形式,以完整表示原子的局域環(huán)境,并引入一個平滑函數(shù),使其在截斷半徑處平滑過渡到零,以保證原子運動的平滑性和連續(xù)性。值得注意的是,ACSF方法給出的是構(gòu)型的絕對描述,非常適用于在高維神經(jīng)網(wǎng)絡(luò)勢中進行原子化學(xué)環(huán)境描述。

Bartok等[50]并未試圖建立局域構(gòu)型的描述,而是直接建立兩個局域構(gòu)型之間距離的對稱性度量,即平滑重疊原子位置(Smooth Overlap of Atomic Position, SOAP)。這種方法的主要依據(jù)是局域原子密度,將兩個局域構(gòu)型之間的“距離”定義為兩個局域密度之間的交疊程度,并將核原子與其余原子的相對位置表示為徑向項和角向項的乘積。對于原子數(shù)量較多的團簇,SOAP描述符具有獨特的優(yōu)勢,通常能得到更準確和穩(wěn)定的勢能面。然而,SOAP給出的并非是某個構(gòu)型的絕對描述,而是構(gòu)型之間的距離,因此并不適合用于深度神經(jīng)網(wǎng)絡(luò),而非常適用于核方法[9]。同樣地,庫侖矩陣(Coulomb matrix)[51]描述并非針對局域原子環(huán)境,而是針對體系中所有的原子,將體系中N個原子的三維坐標轉(zhuǎn)換為一個N×N的矩陣,多用于開放邊界的體系,而非周期性體系的描述。庫倫矩陣保證了坐標的平移和旋轉(zhuǎn)不變性,對庫倫矩陣排序后可以保證交換對稱性。

描述符構(gòu)造是開發(fā)機器學(xué)習勢的基礎(chǔ),只有將描述符和擬合方法有效結(jié)合,才能最大程度地發(fā)揮機器學(xué)習勢的優(yōu)勢。需要注意的是,關(guān)于何種分子描述符在機器學(xué)習勢模型中表現(xiàn)最優(yōu)、從原子坐標中提取哪種特征,目前學(xué)界并無定論[8]。

1.2 數(shù)據(jù)集搭建方法

1.2.1 采樣方法

機器學(xué)習勢的開發(fā)過程中需要首先利用高精度的量子化學(xué)計算方法來構(gòu)建體系的勢能面。從效率和精度的角度出發(fā),通常采用領(lǐng)域中廣泛使用的密度泛函(DFT)方法。由此可見機器學(xué)習勢的精度上限是量子化學(xué)數(shù)據(jù)集的精度。考慮到DFT的計算成本,完全覆蓋復(fù)雜體系高維全局勢能面也是不切實際的,因此開發(fā)機器學(xué)習勢的核心是搭建合理的數(shù)據(jù)集。數(shù)據(jù)集的充分性、多樣性對于機器學(xué)習勢的精度至關(guān)重要。訓(xùn)練集需要盡可能地覆蓋目標體系的勢能面,對不同含能材料組分、多種原子局域環(huán)境進行批量采樣,以得到精度高、拓展性強、普適性高的分子動力學(xué)機器學(xué)習勢函數(shù)模型。常用的勢能面采樣方法(見圖1)包括結(jié)構(gòu)微擾法、分子動力學(xué)(MD)模擬方法和增強采樣方法。其中結(jié)構(gòu)微擾法是對原子坐標施加結(jié)構(gòu)縮放、錯位和隨機擾動,從而得到大量隨機分子結(jié)構(gòu);MD采樣方法采用現(xiàn)有的反應(yīng)力場、經(jīng)驗力場等進行分子動力學(xué)計算,得到分子的運動軌跡和反應(yīng)路徑;增強采樣法則是通過人為施加約束,降低化學(xué)反應(yīng)能壘,提高反應(yīng)發(fā)生的概率。

圖1 常見的采樣方法:結(jié)構(gòu)微擾法、MD采樣法、增強采樣法和虛擬現(xiàn)實MD采樣法[45]

含能材料組分復(fù)雜,通常涉及到3種以上的化學(xué)組分,5種以上化學(xué)元素;且多應(yīng)用于高溫(>3 000 K)、高壓(>10 GPa)等極端場景,其反應(yīng)過程涉及到氣、液、固3個相態(tài)及兩相甚至多相之間的非均相反應(yīng),這些條件決定了含能材料機器學(xué)習勢的開發(fā)難度遠大于其他材料,往往需要更大規(guī)模的數(shù)據(jù)集來覆蓋可能存在的熱力學(xué)條件。因此在進行勢能面采樣時,需要有針對性地構(gòu)造數(shù)據(jù)集。例如,在高溫條件下,分子化學(xué)鍵斷裂生成大量自由基,原子局域環(huán)境的復(fù)雜度顯著增加,因此初始數(shù)據(jù)集中需包括已知的反應(yīng)路徑、過渡態(tài)結(jié)構(gòu)和反應(yīng)產(chǎn)物。為了描述含能材料的高壓構(gòu)象空間,需要在數(shù)據(jù)集中引入不同壓力狀態(tài)的結(jié)構(gòu)構(gòu)象,即材料的狀態(tài)方程(壓力—體積曲線),以描述含能材料在極端場景的力學(xué)特性。

事實上,在傳統(tǒng)分子動力學(xué)模擬的時間尺度內(nèi),體系很難跨越較高的反應(yīng)能壘,在相空間內(nèi)的采樣效率具有明顯的局限性,使用增強采樣方法能夠更高效地對稀有事件(化學(xué)反應(yīng))實現(xiàn)采樣。然而,經(jīng)典的增強采樣方法需要預(yù)先選定反應(yīng)坐標,通過施加偏置勢或偏置力使得體系沿著既定的方向發(fā)生反應(yīng),考慮到含能材料復(fù)雜的分子結(jié)構(gòu)和苛刻的熱力學(xué)條件,潛在反應(yīng)坐標數(shù)量指數(shù)性增長,因此經(jīng)典的增強采樣方法會存在反應(yīng)坐標無法合理確定的困境,以至于對潛在分解反應(yīng)路徑、過渡態(tài)結(jié)構(gòu)等采樣不足。近年,本課題組[52]開發(fā)了一款虛擬現(xiàn)實(virtual reality, VR)分子動力學(xué)模擬軟件MANTA,為含能材料數(shù)據(jù)集的搭建提供了新的采樣工具。研究人員可以借助VR技術(shù)并憑借化學(xué)直覺對目標反應(yīng)路徑開展實時交互,快速地對特定的反應(yīng)過程勢能面進行采樣,將反應(yīng)構(gòu)型實時加入到訓(xùn)練集中,是一種極具應(yīng)用潛力的輔助采樣工具?？傊?勢能面采樣是開發(fā)機器學(xué)習勢的關(guān)鍵步驟,在訓(xùn)練過程中應(yīng)結(jié)合上述手段,盡可能完備地描述化學(xué)反應(yīng)的勢能面空間,以提高勢函數(shù)的準確性。

1.2.2 主動學(xué)習

主動學(xué)習是指在機器學(xué)習勢訓(xùn)練過程中,對目標體系的相空間進行不斷地探索和補充,實現(xiàn)自動拓展數(shù)據(jù)集的一種訓(xùn)練方法[53]。對于復(fù)雜體系,主動強化學(xué)習是一種十分合適的數(shù)據(jù)挖掘、構(gòu)型探索、勢能面拓展方法。以DP-GEN[54]為例,這一工具是在DeePMD-kit的基礎(chǔ)上開發(fā)的同步學(xué)習工作流,其工作流程如圖2所示。

圖2 DP-GEN工作流程

針對初始數(shù)據(jù)集,調(diào)用DeePMD-kit工具,利用不同的神經(jīng)網(wǎng)絡(luò)初始化設(shè)置,訓(xùn)練多個機器學(xué)習勢函數(shù)模型。隨機選取其中一個勢函數(shù)模型進行分子動力學(xué)模擬(調(diào)用LAMMPS等[55]軟件),對運動軌跡中的構(gòu)型能量和受力進行評估。若該模型預(yù)測構(gòu)型的能量或受力與其他幾個機器學(xué)習勢函數(shù)模型有明顯偏差,則表示當前模型對于該構(gòu)型的預(yù)測精度較差,需要補充到初始數(shù)據(jù)集中從而進一步優(yōu)化勢函數(shù)模型。DP-GEN引入模型偏差來判斷多個機器學(xué)習勢間的誤差,通過設(shè)定誤差上下限來挑選候選構(gòu)型。隨機挑選其中部分構(gòu)型進行DFT計算得到新的數(shù)據(jù),并加入到現(xiàn)有的數(shù)據(jù)集中重新訓(xùn)練,從而得到精度比之前更高的機器學(xué)習勢模型。上述流程對數(shù)據(jù)集進行補充挖掘并得到新的機器學(xué)習勢模型,被稱為完成了一輪迭代。

采用DP-GEN方案進行多次迭代,每輪迭代可以改變系綜設(shè)置,完成在目標溫度、壓力或體積條件下的批量采樣和定量篩選,實現(xiàn)主動學(xué)習,旨在有效覆蓋所需樣本空間。高效的迭代訓(xùn)練過程設(shè)計能讓勢能面覆蓋區(qū)域從局域空間平滑穩(wěn)定地拓展到整個構(gòu)象空間,例如從短的模擬時間逐步增加為長時間模擬,從低溫模擬逐步擴展為高溫模擬,從低壓模擬逐步增加為高壓模擬,從非反應(yīng)體系過渡至反應(yīng)體系。

2 基于機器學(xué)習勢的含能材料分子模擬

依據(jù)材料功能,含能材料組分一般可分為五大類,如圖3所示。

圖3 含能材料常見組分

其中,硝胺類是指以碳、氫、氧、氮等元素構(gòu)建的具有爆炸性基團的化合物,是含能材料的主體部分[5]。硝胺類組分以環(huán)狀和籠狀的硝胺硝基化合物為主,如環(huán)三亞甲基三硝胺(RDX)、環(huán)四亞甲基四硝胺(HMX)、六硝基六氮雜異伍茲烷(CL-20)等。這類分子受限于反應(yīng)產(chǎn)物的生成熱,導(dǎo)致其能量密度偏低。高能顆粒(鋁、硼、鎂等)憑借較高的燃燒焓被廣泛用做還原劑以提高含能材料配方體系的能量密度,調(diào)節(jié)其能量輸出特性,是強有力的動力來源和毀傷手段[56-58]。然而,高能顆粒的氧化層往往會限制其能量釋放速率,導(dǎo)致較長的點火延遲時間。在含能復(fù)合材料中加入聚偏二氟乙烯(PVDF)、高氯酸銨(AP)等強氧化劑有助于完全發(fā)揮納米顆粒的增益作用,其燃燒和熱分解性能極大改變了傳統(tǒng)含能材料的釋能機制[59-60]。為了降低含能材料感度,提高其力學(xué)性能,組分配方中通常包含少量黏結(jié)劑作為輔料。黏結(jié)劑組分多為高分子材料,如端羥基聚丁二烯(HTPB)、聚疊氮縮水甘油醚(GAP),這些高分子材料同樣包含一定含能基團以提高整體配方的能量密度。此外,近年來一部分高能量特性、熱穩(wěn)定性、機械感度及力學(xué)性能優(yōu)異的催化劑在含能材料領(lǐng)域中得到了廣泛的應(yīng)用,如氧化石墨烯(GO)、銅(Cu)等[61-62]。

本章將介紹機器學(xué)習勢在含能材料分子模擬方面的研究進展,涉及到硝胺類化學(xué)物(RDX/CL-20/ICM-102)、氧化劑(AP)、金屬添加劑(Al/B)等組分的反應(yīng)動力學(xué)研究。同時,包含了機器學(xué)習勢在碳氫燃料燃燒分子模擬中的進展。上述研究所采用的機器學(xué)習模型均為Zhang等[37]提出的DP模型,簡記為NNP。

2.1 硝胺類含能材料

2.1.1 RDX

黑索金(RDX)是一種典型的含能材料,穩(wěn)定性好、爆炸威力大,被廣泛應(yīng)用于火炸藥領(lǐng)域。Chu等[44]通過CP2K計算[63]將PBE/DZVP級別的第一性原理數(shù)據(jù)作為訓(xùn)練集,基于DP策略開發(fā)機器學(xué)習勢模型對RDX熱解的反應(yīng)動力學(xué)進行研究,并與ReaxFF[16]的結(jié)果進行對比。兩種方法得到的不同溫度下RDX分子數(shù)量變化趨勢相同,即體系分解速率均隨溫度升高而增加;但變化快慢卻大不相同,NNP模型預(yù)測的分解速率與第一性原理計算的結(jié)果一致,約為ReaxFF力場結(jié)果的兩倍。這表明開發(fā)的機器學(xué)習勢函數(shù)能準確描述RDX的熱解反應(yīng),而ReaxFF反應(yīng)力場低估了分解反應(yīng)的反應(yīng)速率。此外,NNP模型和ReaxFF力場通過阿倫尼烏斯公式得到的RDX熱解反應(yīng)活化能分別為105.6 kJ/mol和108.8 kJ/mol,與文獻值119.7 kJ/mol接近,即兩種力場模型均能較為準確地預(yù)測RDX分子發(fā)生單分子分解反應(yīng)的活化能。

RDX熱分解的中間體和平衡產(chǎn)物主要包括NO2、NO、HNO2、N2、H2O、CO2和CO等。機器學(xué)習勢函數(shù)以第一性原理的精度成功復(fù)現(xiàn)文獻結(jié)果,即RDX受熱時N—NO2鍵會首先發(fā)生均相斷裂而生成NO2,中間產(chǎn)物NO2進一步發(fā)生反應(yīng)生成NO和H2O分子。當達到平衡濃度時,體系中N2含量最多,其次是H2O和CO2。與ReaxFF模型的結(jié)果相比,機器學(xué)習勢預(yù)測到更多的NO分子生成量、更快的H2O生成率以及更高的N2平衡濃度,與文獻結(jié)果一致[64]。通過對分子動力學(xué)模擬的軌跡進行分析可以得到RDX熱解的詳細反應(yīng)路徑,如圖4所示。RDX的分解機制可分為兩個階段:(1)N—N鍵斷裂生成NO2分子和RDR中間體;(2)中間體RDR進一步發(fā)生N—N鍵斷裂、開環(huán)反應(yīng)、或氫加成反應(yīng)。其中,關(guān)鍵中間體NO2發(fā)生氫加成反應(yīng)會生成HONO,進而生成NO、N2、OH和H2O等物質(zhì)。

圖4 RDX分解的中間產(chǎn)物及詳細反應(yīng)路徑[44]

近鄰分子間的相互作用在很大程度上會影響甚至改變含能材料的分解路徑[6]。不同密度的RDX分子和RDR中間體發(fā)生N—N均裂、氫提取和開環(huán)反應(yīng)的通量并不相同。在高密度情況下,由于空間位阻效應(yīng),近鄰分子會阻礙RDX分子及其中間體的開環(huán)反應(yīng);隨著密度降低,近鄰分子的影響減弱,開環(huán)反應(yīng)通量增加,氫提取反應(yīng)通量減少,RDX的熱解機制與孤立分子的第一性原理結(jié)果相近。上述研究結(jié)果在保持第一性原理計算精度的同時,從原子尺度闡明了RDX的反應(yīng)路徑,建立了RDX熱分解詳細反應(yīng)動力學(xué)機制。

2.1.2 CL-20

第三代高爆炸藥CL-20是能量和密度最高的單質(zhì)炸藥之一,在氧平衡、爆速、爆壓和堆積密度等方面均優(yōu)于其他含能材料。然而,對機械刺激的高敏感性限制了其應(yīng)用前景。近年來,共晶技術(shù)在含能材料領(lǐng)域已成為一種廣泛應(yīng)用的策略,以更好地平衡感度和能量輸出之間的關(guān)系?，F(xiàn)有的實驗和理論工作表明,CL-20與其他含能材料共晶結(jié)合可以提高沖擊穩(wěn)定性和熱穩(wěn)定性。盡管已經(jīng)有許多關(guān)于CL-20的共晶材料被設(shè)計出來,但提高穩(wěn)定性的具體微觀機理仍有待進一步揭示。

Cao等[46]構(gòu)建了β-CL-20和CL-20/TNT共晶體系PBE/DZVP級別的神經(jīng)網(wǎng)絡(luò)勢能函數(shù),并在此基礎(chǔ)上開展分子模擬研究以揭示CL-20/TNT共晶穩(wěn)定性增強背后的分子機制。為了深入探索β-CL-20和CL-20/TNT共晶的熱解過程,基于機器學(xué)習勢進行了分子動力學(xué)升溫模擬。對于β-CL-20體系,CL-20分子在溫度升高到大約1 200 K時開始分解;而在CL-20/TNT共晶體系中,CL-20分子發(fā)生分解的溫度高達1 400 K左右。通過比較起爆溫度,可以發(fā)現(xiàn)共晶CL-20/TNT的熱穩(wěn)定性明顯提高。在不同的升溫速率下,也能觀察到相似的結(jié)果,即在CL-20/TNT共晶體系中的起爆溫度總比在β-CL-20體系中的要高。這是由于共晶體系會穩(wěn)定生成分子間氫鍵,而在單質(zhì)體系中并沒有氫鍵的出現(xiàn)。分子間的氫鍵相互作用可以有效延遲CL-20分子的裂解,降低活性物種的反應(yīng)速率,這一結(jié)果與張朝陽等[65]對含能材料晶體結(jié)構(gòu)的研究結(jié)果一致。通過對分子動力學(xué)模擬軌跡進行分析,得到單質(zhì)CL-20分子主要的熱分解路徑,如圖5所示。與之前的研究一致,CL-20分子熱解以N—NO2鍵的斷裂為主,脫去NO2后的中間體進一步發(fā)生N—N斷鍵和C—C斷鍵,生成產(chǎn)物NO2和CxNy小分子。在共晶體系中,也觀察到了類似的初始熱解路徑,但是速率較低。TNT中的自由基與游離的NO2分子會發(fā)生進一步反應(yīng),如TNT中的氫自由基與游離的NO2分子結(jié)合生成HONO,TNT中的甲基自由基與NO2結(jié)合生成甲醛(HCHO)和NO。由于TNT中的C—N鍵和六元環(huán)相對穩(wěn)定,因此在模擬過程中很難發(fā)生NO2的分離和開環(huán)反應(yīng)。值得注意的是,從CL-20中分離出的NO2與TNT的反應(yīng)進一步減緩了初始階段NO2的生成速度,阻礙了CL-20分子的進一步分解。因此,TNT可以被視為緩沖碎片,以降低活性中間體與CL-20分子碰撞的概率。上述工作從原子尺度揭示了CL-20/TNT共晶體穩(wěn)定性增強背后的物理化學(xué)機制。

圖5 CL-20熱解的初始反應(yīng)路徑[46]

2.1.3 ICM-102

新型含能材料的設(shè)計是提升我國先進武器裝備綜合性能的有力支撐。張慶華等[66]將“材料基因組工程”的研發(fā)模式引入到了含能材料領(lǐng)域,針對含能材料能量與安全匹配性差的問題,運用計算機輔助分子高通量設(shè)計與篩選的手段,成功研發(fā)了一種高能鈍感炸藥分子ICM-102(C4H6N6O4)。其能量水平高于低感炸藥LLM-105,感度與鈍感炸藥TATB相當[67],實現(xiàn)了炸藥能量與安全性的平衡。Chu等[45]基于虛擬現(xiàn)實增強采樣(VRMD)方案構(gòu)建了PBE/DZVP級別的數(shù)據(jù)集,并開發(fā)機器學(xué)習力場來描述ICM-102的反應(yīng)動力學(xué)。其中,第一性原理計算采樣的構(gòu)型與初始構(gòu)型接近,隨溫度的升高,構(gòu)型在勢能面上的分布區(qū)域擴大,但主要集中于低能區(qū)域,導(dǎo)致反應(yīng)構(gòu)型數(shù)量有限,不能完整描述整個反應(yīng)過程。而VRMD方案通過手柄施加外力,并依靠化學(xué)直覺對經(jīng)驗的反應(yīng)位點、反應(yīng)路徑進行實時交互采樣(見圖2)來構(gòu)建數(shù)據(jù)集,可以直接對斷鍵反應(yīng)等稀有事件進行直接采樣,有效提高了高能壘反應(yīng)過程的采樣效率。

數(shù)據(jù)集的質(zhì)量對于機器學(xué)習模型的性能至關(guān)重要。通過主成分分析可提取采樣結(jié)構(gòu)的主要特征并進行數(shù)據(jù)降維,從而將采樣到的結(jié)構(gòu)劃分為晶體構(gòu)型、部分分解構(gòu)型和氣態(tài)產(chǎn)物,如圖6所示。沿著總反應(yīng)發(fā)生的方向,即主成分PC1增加的方向,兩種采樣方案得到的構(gòu)型質(zhì)量均逐漸減少,但是對應(yīng)的數(shù)密度卻明顯不同。由于反應(yīng)能壘的存在,直接分子動力學(xué)采樣的數(shù)據(jù)集范圍明顯集中在低能勢阱,在高能區(qū)域的采樣較少,即得到的構(gòu)型大多并未涉及到反應(yīng)的發(fā)生。而VRMD采樣方案可針對能壘較高的反應(yīng)過程進行采樣,作為AIMD采樣方式的補充,極大地拓展了勢能面的采樣空間。此外,單一的AIMD采樣方案得到的模型能量精度為0.014 eV/atom,兩種方法相結(jié)合得到的模型偏差僅有0.006 eV/atom。這表明VRMD采樣方案可以作為補充手段,來構(gòu)造更為均衡、更豐富的訓(xùn)練集,從而得到更精準、更有效的勢函數(shù)模型。

圖6 AIMD和VRMD采樣方案的主成分分析[45]

通過高溫分子模擬,Chu等[45]發(fā)現(xiàn)ICM-102熱分解的主要中間體和產(chǎn)物包括C6H5N6O4、H2O、NO、N2、CHNO和CO2。利用聚類分析的方法可以將所有出現(xiàn)過的分子構(gòu)型分為八類:前兩類分子表示ICM-102及其失去一個小分子基團后的產(chǎn)物;第三、四類分子表示ICM-102分子經(jīng)過開環(huán)反應(yīng)的產(chǎn)物;第五、六、七類分子表示開環(huán)反應(yīng)后進一步生成的含碳或含氮的分子;第八類則是反應(yīng)過程中的水、一氧化氮和氮氣等小分子氣體。以熱解反應(yīng)網(wǎng)絡(luò)為基礎(chǔ),可直接依據(jù)通量分析推演任意中間物質(zhì)的詳細反應(yīng)網(wǎng)絡(luò)(見圖7)。

2.2 AP

AP的有效氧含量高、機械感度低、相容性好,爆炸過程時產(chǎn)生的氣體多、無固體殘渣、成本低,是固體推進劑和火炸藥等最重要的強氧化劑[68-69]。AP的熱分解是一個復(fù)雜的物理化學(xué)過程,涉及到吸熱相變和放熱反應(yīng)。目前AP的熱分解研究以傳統(tǒng)實驗為主,如使用熱重分析(TGA)、差示掃描量熱法(DSC)和質(zhì)譜分析(MS)等來確定分解速率和中間產(chǎn)物。傳統(tǒng)的實驗手段并不能在極小的時間、空間尺度上直接觀測到AP復(fù)雜的反應(yīng)機制。而AP反應(yīng)力場的缺失制約了研究者對其熱分解機理的微觀研究,是火炸藥領(lǐng)域中的“卡脖子問題”。Chu等[70]基于第一性原理計算(PBE/DZVP級別)首次建立了AP的機器學(xué)習勢,對其熱分解過程進行分子動力學(xué)模擬,從原子尺度成功提取AP分解的動力學(xué)演化過程,填補了目前含能材料燃燒機理缺失的重要部分。通過分子模擬復(fù)現(xiàn)TGA/DSC實驗,將AP晶體從常溫加熱至3 000 K,所開發(fā)的機器學(xué)習勢可描述AP熱分解過程中的質(zhì)量變化和熱流速率,即從正交相到立方相的吸熱相變,以及隨后的放熱分解反應(yīng)。

AP分子由NH4和ClO4基團構(gòu)成,通過直接分子動力學(xué)模擬,機器學(xué)習勢模型預(yù)測其分解過程中生成的中間體和反應(yīng)產(chǎn)物多達三萬種。其中,反應(yīng)的中間體主要是NH3、ClO2、ClO3、HClO4。隨著反應(yīng)的進行,最后的平衡產(chǎn)物以H2O、NO、O2、Cl2、N2和HCl為主。通過提取熱分解過程中的主要物質(zhì)和生成路徑得到AP反應(yīng)網(wǎng)絡(luò),如圖8所示。

圖8 AP分解的產(chǎn)物數(shù)量變化和主要反應(yīng)路徑[70]

考慮到AP中的強氧化性元素較多,據(jù)此可以將整個分解動力學(xué)劃分為AP的初始分解、氫化學(xué)、氮化學(xué)和氯化學(xué)。初始分解過程以質(zhì)子轉(zhuǎn)移為主,NH4中的氫自由基會轉(zhuǎn)移到HClO4中并生成NH3,或與ClO4中的O自由基結(jié)合生成ClO3和OH。在初級分解之后,氮化學(xué)、氫化學(xué)和氯化學(xué)反應(yīng)同時進行。特別指出的是,上述分子模擬預(yù)測到一些實驗未曾報道過的雙分子反應(yīng),例如NH4+NH2→2NH3和ClO4+NH3→HClO4+NH2。這些氣相物種和AP之間的雙分子反應(yīng)可能會影響AP的分解過程。上述研究以第一性原理計算的精度首次揭示了AP的詳細反應(yīng)機制,克服了AP反應(yīng)力場缺失對基礎(chǔ)研究的制約,為建立復(fù)合含能體系的反應(yīng)動力學(xué)模型提供了新的機遇。

2.3 高能顆粒(Al、B)

金屬顆粒能量密度大、燃燒熱值高,是極具吸引力的含能材料組分,能顯著提高固沖發(fā)動機的性能,增大導(dǎo)彈武器的射程,提高戰(zhàn)術(shù)武器的作戰(zhàn)效能和毀傷效益[40, 71-72]。金屬作為燃料和燃料添加劑可以減少燃燒不穩(wěn)定性、提高燃燒速率、降低爆震敏感性[57, 73]。其中應(yīng)用最廣的金屬顆粒是鋁,其燃燒熱值高、能量密度大且儲量大,是重要的含能材料添加劑。近期本課題組基于PBE/DZVP級別的第一性原理數(shù)據(jù)建立機器學(xué)習勢模型,訓(xùn)練集中包括Al、Al2O3、界面以及氧氣擴散等數(shù)據(jù)。所開發(fā)的模型可以復(fù)現(xiàn)晶體鋁的狀態(tài)方程曲線(見圖9),計算結(jié)果與DFT預(yù)測基本一致。模型在極端的壓縮和拉伸情況下仍然保持著良好的性能表現(xiàn),因此可以用于研究鋁在沖擊波下的力學(xué)響應(yīng)行為。值得注意的是,盡管訓(xùn)練集中僅包含了單位原子體積在10～25 ?3/atom的數(shù)據(jù)(圖中虛線區(qū)域的數(shù)據(jù)),機器學(xué)習勢模型可以準確地預(yù)測7～32 ?3/atom范圍內(nèi)的狀態(tài)方程,這證明了機器學(xué)習勢函數(shù)具有一定的外推性能。

圖9 鋁的狀態(tài)方程曲線

在所有金屬和類金屬中,除有毒性的鈹外,硼的體積熱值和質(zhì)量熱值最高,分別為138 kJ/cm3和59 kJ/g[74-75]。盡管硼的理論燃燒性能優(yōu)異,但是硼的熔點(2 348 K)和沸點(4 043 K)都很高,導(dǎo)致顆粒難以點火,限制了硼粉的實際應(yīng)用。此前Liu等[76]、Wang等[77]采用分子模擬的方式研究了硼顆粒的微觀燃燒機理,然而這些工作使用的勢函數(shù)是基于氨硼烷[78]開發(fā)的,并未針對硼燃燒進行優(yōu)化。這里我們提出一個問題:一個力場模型的開發(fā)數(shù)據(jù)集中并不包括硼單質(zhì)的數(shù)據(jù),這樣的模型可直接用于硼燃燒的模型嗎?勢函數(shù)模型是針對目標體系特定的應(yīng)用場景而開發(fā)的。一般而言,訓(xùn)練集只是全局勢能面上的局部空間,在勢能面上的外推性能非常局限。例如,同樣含有C/H/O/N元素的反應(yīng)力場可以描述RDX的熱解反應(yīng),卻不能直接用于CL-20的模擬,會導(dǎo)致預(yù)測誤差迅速升高。因此,在進行分子動力學(xué)模擬之前,需要對勢函數(shù)開發(fā)的訓(xùn)練集進行分析,確保數(shù)據(jù)集覆蓋了關(guān)注的勢能面區(qū)域,避免超出訓(xùn)練集的外推空間,并進行物理化學(xué)參數(shù)驗證,如晶格參數(shù)、狀態(tài)方程、融點、比熱容等。

Chang等[79]利用CP2K構(gòu)建了PBE/DZVP級別的量子化學(xué)數(shù)據(jù)集,開發(fā)了用于描述硼和硼氧化合物的DP模型。開發(fā)的DP模型對能量、受力、狀態(tài)方程、晶格參數(shù)和徑向分布函數(shù)等進行驗證,機器學(xué)習力場的精度均和第一性原理計算一致。相比之下,針對氨硼烷燃燒開發(fā)的ReaxFF力場預(yù)測結(jié)果非常差[78]。例如,文獻中晶體硼的單位原子體積為7.25 ?3/atom,機器學(xué)習預(yù)測值為7.01 ?3/atom,而ReaxFF力場的預(yù)測值高達12.58 ?3/atom,顯著高估了晶體硼的平衡結(jié)構(gòu)。以上這種明顯誤差主要來源于ReaxFF力場的濫用,考慮到上述工作的ReaxFF是針對氨硼烷燃燒體系,在其開發(fā)過程中不含任何硼單質(zhì)晶體結(jié)構(gòu),因此這套ReaxFF力場參數(shù)無法支持硼燃燒的分子動力學(xué)模擬。另外,這直接表明研究人員應(yīng)當避免任意拓展ReaxFF力場的應(yīng)用范圍。

Chang等[79]采用DP模型模擬了納米硼顆粒的融化行為。通過計算顆粒全局Lindemann指數(shù)和勢能變化,得到不同粒徑硼顆粒的熔點,歸一化后的熔點如圖10所示。對于3 nm硼顆粒,其融點高達1 992 K,而同樣粒徑的鋁鎂顆粒在該溫度下已完全融化并蒸發(fā)到氣相環(huán)境中[80-81]。這表明納米硼顆粒的燃燒機制仍然是以表面反應(yīng)為主。通過分析原子Lindemann系數(shù)演化,作者發(fā)現(xiàn)硼融化的過程包括初始加熱、表面預(yù)融、內(nèi)核預(yù)融、完全融化和繼續(xù)升溫等五個階段。該融化過程遵循Nanda等[82]提出的液核生成假設(shè),即表面先熔化隨著溫度逐漸向內(nèi)核區(qū)域擴散。上述融化機制表明,在達到融點之前,表面預(yù)融的硼顆粒會和周圍的顆粒發(fā)生燒結(jié)、團聚,阻礙氧化劑向顆粒內(nèi)部擴散,從而導(dǎo)致點火延遲?；跈C器學(xué)習勢的融點研究從原子角度解釋了硼顆粒在燃燒爆炸應(yīng)用時受限的物理化學(xué)機制。

圖10 鋁、鎂、硼的歸一化融點[79]

2.4 碳氫燃料

碳氫燃料具有高能量密度和低成本等特性,是火箭推進劑、噴射推進劑等航空燃油的主要組成部分,其反應(yīng)機制是燃燒反應(yīng)動力學(xué)中重要的研究對象。Zeng等[83]將深度勢模型引入燃燒領(lǐng)域研究,基于MN15/6-31G**級別的第一性原理數(shù)據(jù)開發(fā)機器學(xué)習勢,對甲烷氧化的基本反應(yīng)機理和動力學(xué)過程開展研究。分子模擬結(jié)果表明,甲醛(CH2O)和甲氧基(CH3O)是最關(guān)鍵的中間產(chǎn)物,主要反應(yīng)路徑為:CH4→CH3→CH3O→CH2OH→H2CO→HCO→CO。甲烷燃燒起始反應(yīng)是氧氣的抽氫反應(yīng),生成甲基(CH3)和HOO自由基。隨著反應(yīng)進行,不斷生成的OH、H和HOO自由基會與甲烷發(fā)生抽氫反應(yīng),生成大量的甲基。甲基進一步被氧化為甲醛,接著失去一個氫原子形成甲?；?HCO—),隨后可以通過碰撞解離,與氧分子反應(yīng)失去氫原子,從而形成一氧化碳(CO)。上述反應(yīng)路徑與文獻結(jié)果一致,以第一性原理的精度搭建了甲烷燃燒的反應(yīng)網(wǎng)絡(luò)。

高溫下碳氫燃料的不完全燃燒是碳煙前驅(qū)體多環(huán)芳香烴(Polycyclic aromatic hydrocarbons,PAH)的重要來源。碳煙的初生受燃料結(jié)構(gòu)、氧氣濃度、溫度等因素的影響,同時存在物理和化學(xué)機制的復(fù)雜過程。Wang等[84]通過Gaussian計算得到MN15/6-31G**級別的數(shù)據(jù)集從而搭建機器學(xué)習勢函數(shù)模型,并根據(jù)化學(xué)經(jīng)驗設(shè)計反應(yīng)物,對A1～A9多環(huán)芳烴化合物(PAH)生成路徑進行模擬研究。以單環(huán)和雙環(huán)PAH的生成過程為例,機器學(xué)習勢函數(shù)成功預(yù)測到文獻中丁二烯和乙炔反應(yīng)生成苯環(huán)、五元環(huán)碳碳鍵斷裂生成萘的過程。在生成三環(huán)芳香烴化合物的模擬中,觀測到的主要反應(yīng)為經(jīng)典的抽氫-乙炔加成(HACA)反應(yīng),即乙炔與1-萘結(jié)合生成苊并致使五元環(huán)閉合。此外,還發(fā)現(xiàn)一些新的反應(yīng)路徑,如苯基和C4H4反應(yīng)生成一條長碳鏈,然后通過閉環(huán)和抽氫反應(yīng)生成茚。這些新的反應(yīng)路徑對于構(gòu)建詳細的碳煙生成模型具有重要意義。

Chu等[85]基于DP方案,發(fā)展了PBE/DZVP精度級別的機器學(xué)習勢函數(shù),在此基礎(chǔ)上進一步建立了納米反應(yīng)器算法以加速反應(yīng)、探究密度對于反應(yīng)路徑的影響。納米反應(yīng)器用一個球形壓頭周期性地壓縮和放松體系,使體系的密度發(fā)生變化,從而改變分子碰撞頻率,加快反應(yīng)進程。此外,高溫手段使得分子動能增加,運動更為劇烈,同樣可以顯著加快反應(yīng)進程。圖11以乙炔體系多聚化反應(yīng)為例,對比了常規(guī)分子動力學(xué)計算結(jié)果和納米反應(yīng)器算法加速的計算結(jié)果。

圖11 3種分子動力學(xué)模擬過程中乙炔分子數(shù)量變化及物種演化[85]

其中,方案一表示納米反應(yīng)器中的反應(yīng)進程,方案二和方案三分別為1 500 K和3 000 K下傳統(tǒng)的NVT模擬。值得注意的是,納米反應(yīng)器算法可以顯著加快反應(yīng)進行,其反應(yīng)速率是同樣溫度MD結(jié)果的50倍。納米反應(yīng)器中的乙炔分子數(shù)量變化與1 500 K下縮放的結(jié)果趨勢相同,其反應(yīng)途徑相同,僅僅存在乙炔的聚合反應(yīng)。這表明,增加濃度或壓力只影響碰撞頻率,而不影響反應(yīng)路徑。高溫會影響體系的總能量,進而改變反應(yīng)進程,因此方案三中的反應(yīng)路徑與前兩種方案的結(jié)果并不相同,涉及到氫抽離和氫加成反應(yīng)。上述結(jié)果表明納米反應(yīng)器在不影響反應(yīng)路徑的前提下提高了碰撞頻率,為研究成核和低溫反應(yīng)等緩慢動力學(xué)過程的反應(yīng)機理提供了一種很有前景的工具。

作者進一步采用納米反應(yīng)器方法來加速碳煙的成核過程,考慮了文獻中報道的7種不同PAH分子的成核反應(yīng)進程。圖12展示了從納米反應(yīng)器中的不同前體收集的最終產(chǎn)物。對于小分子前驅(qū)體,最終生成的多環(huán)芳香烴以氣態(tài)形式排列。大多數(shù)自由基會形成交聯(lián)的二聚體或三聚體,并未觀察到大的團簇。質(zhì)量較大的前驅(qū)體分子最終會在納米反應(yīng)器的中心形成團簇。而對于大型前驅(qū)體,其最終產(chǎn)物呈球形堆積結(jié)構(gòu),不同層之間存在交聯(lián)。上述研究表明,碳煙的最終產(chǎn)物形態(tài)與前驅(qū)體的大小有很大關(guān)系。隨著前驅(qū)體尺寸的增加,產(chǎn)物顆粒從無序的團簇變成了有序的緊密堆積結(jié)構(gòu)。

圖12 不同前驅(qū)體在納米反應(yīng)器中生成的最終產(chǎn)物[85]

3 機器學(xué)習勢的計算優(yōu)勢

原子間相互作用勢的函數(shù)形式、理論假設(shè)和近似合理性決定了其精確度,進而影響分子動力學(xué)模擬結(jié)果的可靠性[86]。一般的力場假設(shè)的函數(shù)形式較為簡單,如施加兩體或三體截斷。目前含能材料分子動力學(xué)研究常采用的ReaxFF力場,是僅對能量項進行擬合,并未考慮原子受力項,其函數(shù)形式復(fù)雜,并行效率低。機器學(xué)習勢是利用其底層算法對高維函數(shù)的強大擬合能力,以一種基于數(shù)據(jù)驅(qū)動的方式近似第一性原理原子相互作用。圖13總結(jié)了目前含能材料和碳氫燃料中的機器學(xué)習勢模型精度。其中,Zeng等[83]開發(fā)的甲烷燃燒機器學(xué)習勢是基于Gaussian計算得到MN15/6-311G**精度的第一性原理數(shù)據(jù),其余機器學(xué)習勢的數(shù)據(jù)集均為CP2K計算的PBE/DZVP級別精度。

圖13 NNP勢和ReaxFF反應(yīng)力場精度對比

總體來說,各個體系中的機器學(xué)習勢表現(xiàn)均優(yōu)于ReaxFF力場,能量和受力的預(yù)測精度分別提高3.5倍和6.6倍(RDX、Al和碳煙體系)。最初的ReaxFF力場開發(fā)是針對碳氫燃燒并逐步拓展到CHON類含能材料研究,經(jīng)過近二十年的驗證和優(yōu)化,其參數(shù)精度較高、力場表現(xiàn)較好。從圖13中可以看出,經(jīng)過DP-GEN迭代的RDX和碳氫的勢函數(shù)模型精度都僅僅是略高于ReaxFF。以RDX為例,在能量預(yù)測方面,NNP力場相較于ReaxFF力場能量精度提高4倍,受力精度提高9倍。這是由于ReaxFF擬合過程僅考慮了能量項,并未考慮受力項,所以對于原子受力的偏差較大。值得注意的是,對于反應(yīng)力場訓(xùn)練集和目標體系并不匹配的力場而言,ReaxFF的表現(xiàn)不盡人意。以硼的機器學(xué)習勢為例,ReaxFF力場的能量和受力偏差分別是NNP模型的180倍和18倍。這樣顯著的偏差強調(diào)了力場訓(xùn)練集對于目標體系的適配性和移植性。此外,勢能面的采樣質(zhì)量對于機器學(xué)習勢的精度也至關(guān)重要。對于新型含能材料ICM-102,VRMD采樣方案可對勢能面上的高能區(qū)域進行采樣,從而提升數(shù)據(jù)集的整體質(zhì)量。采用VRMD采樣方案使得能量預(yù)測偏差降低58%,受力預(yù)測偏差降低37%,這種輔助采樣策略極大地提高了模型的整體精度。

ReaxFF力場在擬合過程中未考慮材料在極端條件下發(fā)生壓縮、膨脹等力學(xué)行為。在強烈的爆炸沖擊作用下,含能材料的溫度會瞬間急速上升,壓力升高至幾十GPa,體積減少30%[87]。機器學(xué)習勢可以在訓(xùn)練集中加入不同壓縮和拉伸條件下的高溫高壓數(shù)據(jù)集,這對于研究含能材料的沖擊壓縮效應(yīng)具有現(xiàn)實意義。圖14為不同計算方法得到的RDX晶體狀態(tài)方程。其中,機器學(xué)習勢預(yù)測結(jié)果與第一性原理計算結(jié)果一致。

圖14 RDX晶體的狀態(tài)方程曲線[44]

值得注意的是,機器學(xué)習勢具有一定的外推性能。在圖14中,機器學(xué)習勢的訓(xùn)練集僅考慮了0.9～1.0的壓縮比,即圖中的矩形部分。但是該模型能成功預(yù)測0.8～0.9區(qū)間的狀態(tài)方程曲線。相比之下,反應(yīng)力場由于并未考慮勢能面上的壓縮區(qū)域,無法復(fù)現(xiàn)RDX的狀態(tài)方程曲線,預(yù)測的平衡結(jié)構(gòu)體積偏大,在進行弛豫后將無法得到準確的晶體構(gòu)型。此外,反應(yīng)力場明顯高估了壓縮過程中的勢能,在爆炸或沖擊載荷的研究中可能會人為形成局部熱點,加速反應(yīng)的發(fā)生。因此,目前的反應(yīng)力場的精度并不能滿足爆炸沖擊下材料性能的研究需求。Chu等[44]提出的機器學(xué)習力場更能準確描述RDX晶體在燃燒爆炸等極端環(huán)境下的理化性質(zhì)和分解過程。

除了擁有更高的計算精度,機器學(xué)習勢計算效率也更為優(yōu)異。圖15以RDX晶體為例對比了不同方法的計算效率,體系規(guī)模從100個至10萬個原子。機器學(xué)習勢在保持DFT精度的基礎(chǔ)上,計算效率與體系規(guī)模呈線性關(guān)系,約為反應(yīng)力場的27倍,比DFT效率提升了4個數(shù)量級。這是由于機器學(xué)習勢模型計算是基于GPU硬件實現(xiàn)的。對于RDX晶體而言,機器學(xué)習勢的計算成本與體系規(guī)模呈O(N)的趨勢。而對于密度較低的氣相碳反應(yīng)體系,機器學(xué)習勢模型的計算成本與體系規(guī)模呈O(NlogN)的趨勢。機器學(xué)習勢的優(yōu)異表現(xiàn)使得用第一性精度探索具有數(shù)十萬甚至數(shù)千萬原子系統(tǒng)成為可能,可以從原子角度研究含能材料的復(fù)雜反應(yīng)機制,為建立詳細的含能材料反應(yīng)動力學(xué)模型提供了新的機遇。

圖15 密度泛函理論(DFT)、反應(yīng)力場(ReaxFF)和機器學(xué)習勢(NNP)的計算效率對比[44]

4 挑戰(zhàn)與展望

分子動力學(xué)模擬能夠在原子尺度上揭示含能材料的微觀化學(xué)反應(yīng)路徑,為研究復(fù)雜體系的反應(yīng)動力學(xué)機理提供了一條新的途徑。分子動力學(xué)模擬的關(guān)鍵是力場參數(shù),它能夠描述反應(yīng)的勢能面結(jié)構(gòu),決定了模擬的精度和效率。從頭算分子動力學(xué)準而不快,反應(yīng)力場計算快而不準,機器學(xué)習勢兼顧計算精度和計算效率,在含能材料燃燒與爆炸領(lǐng)域得到廣泛應(yīng)用。

機器學(xué)習算法不需要預(yù)設(shè)物理化學(xué)規(guī)律,可直接從計算數(shù)據(jù)出發(fā),準確推斷函數(shù)關(guān)系。得益于這一優(yōu)點,機器學(xué)習勢往往能夠重現(xiàn)含能材料反應(yīng)過程勢能面變化,進而研究其熱解或燃燒過程的詳細反應(yīng)機理。同時,數(shù)據(jù)驅(qū)動這一特性也給機器學(xué)習勢的應(yīng)用帶來了一系列的挑戰(zhàn)。機器學(xué)習勢的開發(fā)需要大量高質(zhì)量數(shù)據(jù),且很難先驗地判斷數(shù)據(jù)集的完備性。雖然,數(shù)據(jù)集的規(guī)?？赡芎艽?大于10萬量子化學(xué)計算數(shù)據(jù)),但它并不足以對目標體系進行充分采樣。例如,基于低溫分子模擬采樣搭建的數(shù)據(jù)集,其機器學(xué)習勢無法準確外推到高溫體系。因此,數(shù)據(jù)集需要盡可能地包括目標體系的整個勢能空間。其次,訓(xùn)練集一般是基于第一性原理計算開發(fā),而第一性原理計算本身也存在一定誤差,受到基組大小、泛函精度等因素的限制。通過提高第一性原理計算的理論方法可以減少計算誤差,然而也會極大地增加計算耗時,尤其是考慮到現(xiàn)有的勢函數(shù)開發(fā)往往需要十萬量級的第一性原理計算數(shù)據(jù),搭建高精度數(shù)據(jù)集的成本極其昂貴。未來含能材料的機器學(xué)習勢函數(shù)開發(fā)面臨著如下挑戰(zhàn):

一是如何對極端條件下復(fù)雜反應(yīng)勢能面進行充分采樣。此前的工作主要針對單質(zhì)體系,在實際應(yīng)用中含能材料的配方日益復(fù)雜,包含硝胺類材料、高分子化合物、金屬材料等多種組分。隨著元素種類的增加,反應(yīng)空間構(gòu)型的數(shù)量呈現(xiàn)指數(shù)級增長,需要包含更多的反應(yīng)路徑、中間體構(gòu)型。這對如何有效地構(gòu)建訓(xùn)練集提出了新的要求。目前存在主動學(xué)習、勢能面掃描等思路。Zhang等[54]提出了DP-GEN主動學(xué)習算法,通過比較不同模型的不確定性動態(tài)生成新的數(shù)據(jù)點,從而提高模型對全局勢能面的覆蓋度。另一方面,研究者可以根據(jù)對化學(xué)反應(yīng)的認知,設(shè)計勢能面掃描算法對化學(xué)反應(yīng)路徑上的結(jié)構(gòu)構(gòu)型進行采樣。Liu等[88]提出采用隨機勢能面行走方法(SSW),利用化學(xué)鍵振動對構(gòu)型施加微擾,從而高效搜索復(fù)雜結(jié)構(gòu)的全局勢能面。Zhang等[89]提出了一種準經(jīng)典軌跡法構(gòu)建局部勢能面的方法,可針對特定反應(yīng)開發(fā)高精度勢能面。增強采樣方法也是一種高效采樣方法,可以作為常規(guī)采樣方法的重要補充。本課題組最近結(jié)合虛擬現(xiàn)實技術(shù)開發(fā)了一種增強采樣工作流[52],在虛擬現(xiàn)實空間中直接對關(guān)鍵反應(yīng)路徑的勢能面進行采樣。和傳統(tǒng)的分子動力學(xué)模擬直接采樣相比,上述幾種方法都能夠提高勢能面采樣的效率,克服勢能面上存在的高勢壘,更好地描述復(fù)雜反應(yīng)過程。值得注意的是,含能材料應(yīng)用中常涉及極端的高溫高壓環(huán)境,在構(gòu)建訓(xùn)練數(shù)據(jù)集時需要把極端條件的構(gòu)型納入考慮。

二是如何提高機器學(xué)習勢訓(xùn)練集的理論精度。機器學(xué)習勢的精度依賴于訓(xùn)練集選取的量子化學(xué)計算方法,特別是特定的物理化學(xué)性質(zhì)或反應(yīng)體系對DFT的泛函精度要求較高(如反應(yīng)能壘、生成焓等),需要使用高精度的雜化甚至雙雜化泛函計算方法。然而高精度的DFT計算成本隨著系統(tǒng)規(guī)模的擴大而急劇增加,這限制了高精度DFT數(shù)據(jù)的數(shù)量,難以滿足機器學(xué)習勢開發(fā)的需求。Chen等[90]提出了DeePKS方法,通過機器學(xué)習方法開發(fā)DFT泛函,從而極大地縮減了高精度DFT計算的耗時,有望實現(xiàn)大規(guī)模的高精度DFT訓(xùn)練集構(gòu)建。Smith等[91]提出了一種程序框架實現(xiàn)了小樣本數(shù)據(jù)的機器學(xué)習勢函數(shù)開發(fā)。他們先使用普通的DFT數(shù)據(jù)訓(xùn)練模型,然后使用遷移學(xué)習方法在高精度的DFT數(shù)據(jù)上對模型進行再訓(xùn)練,從而獲得高精度的機器學(xué)習勢函數(shù)。在未來開發(fā)含能材料的機器學(xué)習勢函數(shù)時,結(jié)合上述方法有希望更好地構(gòu)建高精度DFT數(shù)據(jù)集。