基于機(jī)器學(xué)習(xí)構(gòu)建的環(huán)三亞甲基三硝胺晶體勢*

王鵬舉 范俊宇 蘇艷 趙紀(jì)軍

(大連理工大學(xué), 三束材料改性教育部重點實驗室, 大連 116024)

環(huán)三亞甲基三硝胺(RDX)是一種高能低感度炸藥, 對其能量和性質(zhì)的準(zhǔn)確計算對于開展該炸藥的分子模擬至關(guān)重要. 本文基于機(jī)器學(xué)習(xí)算法, 采用高維神經(jīng)網(wǎng)絡(luò)模型, 對RDX 分子晶體結(jié)構(gòu)數(shù)據(jù)集進(jìn)行勢函數(shù)訓(xùn)練. 分別采用9 種不同的網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行測試訓(xùn)練, 并選取其中學(xué)習(xí)效果最好的勢函數(shù)對RDX 分子晶體結(jié)合能和晶格中原子受力進(jìn)行計算, 均能很好地重復(fù)出第一性原理的計算結(jié)果, 其測試集結(jié)合能的均方根誤差為59.2 meV/atom. 作為機(jī)器學(xué)習(xí)勢函數(shù)的應(yīng)用, 進(jìn)一步使用該勢函數(shù)對a 相RDX 晶體進(jìn)行分子動力學(xué)模擬,以驗證其適用性.

1 引 言

含能材料通常是由C, H, N 和O 四種元素在一定條件下形成的分子晶體, 可在短時間內(nèi)釋放出巨大能量, 被廣泛地應(yīng)用于常規(guī)武器裝備的發(fā)射藥、推進(jìn)劑、炸藥以及火工煙火藥劑等[1]. 典型的含能材料有三硝基甲苯(C7H5N3O6, TNT)、環(huán)三亞甲基三硝胺(C3H6N6O6, RDX)、環(huán)四亞甲基四硝胺(C4H8N8O8, HMX)等. 其中, RDX 俗稱黑索金,是一種常見的含能材料, 因其具有較好的熱穩(wěn)定性、高能量和高爆速等特點, 常用于各種軍用炸藥和推進(jìn)劑. 目前, 已知的RDX 晶體共有五種晶型,分別為a,β,γ,δ和ε相. 其中,a相是常溫常壓相;β相為亞穩(wěn)相[2?7], 在477 K 時可發(fā)生a到β相的相變[2];γ和δ相為高壓相, 分別在約4.0 和18.8 GPa下出現(xiàn)[8?14], 但δ相至今并沒有人得到其準(zhǔn)確結(jié)構(gòu);ε相為高溫高壓相, 出現(xiàn)的溫度和壓強(qiáng)分別為489 K和4.2 GPa[15].

含能材料的起爆過程是一個涉及高溫和高壓等極端條件的瞬態(tài)過程, 許多現(xiàn)象和細(xì)節(jié)仍然無法在實驗上直接獲得. 計算模擬研究因具有周期短、成本小和安全等特點被廣泛用于獲取含能材料的基本性質(zhì)、狀態(tài)方程和爆轟性能等. 相比于第一性原理計算所需求的大量時間和計算資源, 經(jīng)驗勢常被應(yīng)用于含能材料的理論研究. 為獲取含能材料晶體的基本性質(zhì), Sorescu 等[16]首先針對a相RDX晶體發(fā)展了一種分子間相互作用勢用于結(jié)構(gòu)和性質(zhì)的描述. 將該勢函數(shù)與晶格堆垛計算結(jié)合,得到的晶體結(jié)構(gòu)和晶格能都與實驗符合較好. 隨后, 他們又將這種相互作用勢用于30 多種其他硝銨分子晶體的預(yù)測, 發(fā)現(xiàn)當(dāng)原子使用包含電子關(guān)聯(lián)效應(yīng)的第一性原理計算的電荷時, 大部分體系結(jié)構(gòu)與實驗相一致, 但晶格能有較大差異, 并通過給這些原子電荷量乘以一個尺度因子, 來降低這種差異[17].

相較于靜態(tài)計算模擬, 分子動力學(xué)模擬不僅能計算不同溫度、壓力下含能材料的基本性質(zhì), 還能觀測到一定條件下的分解過程. 常見的用于分子動力學(xué)模擬的經(jīng)驗勢有Lennard-Jones 勢[18]、CHARMM 勢[19,20]、Compass 勢[21]等. Liu 等[19]使 用CHARMM 力場對液態(tài)硝基甲烷進(jìn)行分子動力學(xué)模擬. 結(jié)合Murnaghan 狀態(tài)方程, 其體彈性模量和一階導(dǎo)數(shù)計算結(jié)果與實驗符合很好[22], 而通過Hogoniot 曲線計算得到的壓力卻高于實驗值. 因傳統(tǒng)經(jīng)驗勢依賴于體系中已存在的化學(xué)鍵, 無法對舊化學(xué)鍵的斷裂和新化學(xué)鍵的生成進(jìn)行較好描述,Strachan 等[23]使用ReaxFF 反應(yīng)力場[24]研究RDX在沖擊波下的分解過程, 發(fā)現(xiàn)低沖擊波速(< 6 km/s)下, RDX 分子中的N—N 鍵斷裂為起始分解的主導(dǎo)過程, 隨著波速增加和時間推移, N2和OH 成為主要產(chǎn)物. 此后, 人們不斷優(yōu)化ReaxFF 力場參數(shù),并研究撞擊等條件[25,26]對含能材料的影響. 鑒于現(xiàn)有經(jīng)驗勢的參數(shù)對體系具有較強(qiáng)的依賴性, 其可移植性較差, 需要更加泛用的勢函數(shù)來對含能材料進(jìn)行準(zhǔn)確描述.

近年來, 機(jī)器學(xué)習(xí)方法在物理和材料等領(lǐng)域得到了廣泛的應(yīng)用[27]. 通過大量的實驗或計算得到的數(shù)據(jù), 機(jī)器學(xué)習(xí)方法可以構(gòu)建從結(jié)構(gòu)、能量到性質(zhì)的直接映射關(guān)系, 從而直接預(yù)測未知體系的性質(zhì), 以節(jié)省第一性原理計算所需要的大量計算時間.目前常見的機(jī)器學(xué)習(xí)算法有人工神經(jīng)網(wǎng)絡(luò)(ANN)[28]、高維神經(jīng)網(wǎng)絡(luò)(HDNN)[29]、庫侖矩陣[30]、高斯近似勢(GAP)[31]和核嶺回歸(KRR)[32]等方法. 針對分子體系, 通常使用這些方法對有機(jī)分子[30,33,34]、芳香族化合物分子晶體[35]以及分子團(tuán)簇[36]等進(jìn)行了大量研究, 并構(gòu)建了可以準(zhǔn)確地描述體系能量和性質(zhì)的勢函數(shù). 此外, Elton 等[37]使用多種機(jī)器學(xué)習(xí)方法對不同含能材料分子的爆壓和爆速等性質(zhì)進(jìn)行預(yù)測, 其結(jié)果均與實驗結(jié)果相近. 基于此, 我們將機(jī)器學(xué)習(xí)方法用于含能材料RDX 晶體的研究, 采用Behler 和Parrinello[29]提出的HDNN 方法, 對第一性原理計算得到的15199 個RDX 晶體結(jié)構(gòu)進(jìn)行勢函數(shù)擬合, 并對不同神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)所擬合得到的勢函數(shù)進(jìn)行評估. 通過對大量RDX 晶體數(shù)據(jù)集的學(xué)習(xí)和機(jī)器學(xué)習(xí)模型的參數(shù)擬合, 以期獲得可以準(zhǔn)確描述晶體中分子內(nèi)和分子間相互作用的勢函數(shù), 為不同RDX 晶型能量和性質(zhì)的準(zhǔn)確計算提供新的研究方法.

2 研究方法

2.1 數(shù)據(jù)集獲取

首先, 使用Monte Carlo 和模擬退火算法[38,39]結(jié)合Mors6-SCFF 力場[40]對RDX 晶體進(jìn)行結(jié)構(gòu)搜索, 以期在已知的a,β,γ和ε相之外得到更多的穩(wěn)定晶型, 模擬晶胞中均包含8 個RDX 分子.之 后,使 用VASP (Viennaab-initiosimulation package)[41]軟件包對搜索得到的晶型進(jìn)行結(jié)構(gòu)優(yōu)化. 選取廣義梯度近似(GGA)下的PBE (Perdew-Burke-Ernzerhof)[42]泛函和綴加平面波(PAW)[43],能量截斷為1000 eV,k點設(shè)置為2 × 2 × 2. 優(yōu)化后共得到248 種不同的RDX 晶型結(jié)構(gòu).

為獲取更多的數(shù)據(jù)進(jìn)行訓(xùn)練, 繼續(xù)對這248 種結(jié)構(gòu)進(jìn)行加壓優(yōu)化、改變晶格參數(shù)、二面角大小、斜切以及對部分原子位移等操作產(chǎn)生新的結(jié)構(gòu), 并使用上述相同的設(shè)置進(jìn)行能量和原子受力計算. 剔除部分不合理的結(jié)構(gòu)后, 最終得到15199 個用于訓(xùn)練的結(jié)構(gòu).

2.2 機(jī)器學(xué)習(xí)模型構(gòu)建

使用HDNN[29]模型在aenet (atomic energy network)[44,45]軟件包進(jìn)行機(jī)器學(xué)習(xí)計算. 這種方法對體系中每一個原子的能量進(jìn)行單獨計算, 所有原子能量的加和為體系的總能量. 因此, 首先通過每個原子和周圍不同種類原子的位置關(guān)系來構(gòu)建基函數(shù)[29], 公式如下所示:

其中i,j和k代表不同原子;Rij,Rik和Rjk代表原子間距離;θijk為i,j,k三個原子連線的夾角.Rs為偏移參數(shù), 取值為0. (1a)式中的η和(1b)式中的ζ,λ,η為可調(diào)節(jié)參數(shù), 取值如表1 和表2 所列.fc為平滑截斷函數(shù), 表達(dá)式為[29]:

Rc為截斷距離, 在本文中設(shè)置為6.5 ?. 因RDX晶體中有H, C, N 和O 4 種原子, 對于(1a)式中的每一種η的取值, 有4 個不同的基函數(shù), 而對于(1b)式中每一組參數(shù)取值, 會有10 個不同的基函數(shù). 因此, 對于每一種原子, 共得到72 個基函數(shù)結(jié)果作為輸入層神經(jīng)元. 隱藏層的層數(shù)和每層神經(jīng)元的個數(shù)是神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)的超參數(shù), 需要手動調(diào)整來獲得最佳學(xué)習(xí)效果. 因此, 我們分別采用雙隱層

(10-10, 20-20, 30-30), 三隱層(10-10-5, 20-20-10,30-30-10)和四隱層(10-10-10-5, 20-20-10-10, 20-20-20-10)共9 種隱藏層結(jié)構(gòu)進(jìn)行神經(jīng)網(wǎng)絡(luò)構(gòu)建.輸出層只包含一個神經(jīng)元, 代表該原子的能量, 所有原子能量加和為體系總能量, 其表達(dá)式為[46]

表1 4×3 個G2 型徑向基函數(shù)((1a)式)中η 取值Table 1. η of 4 × 3G2 type radial basis functions(Eq. (1a)).

表2 10 × 6 個G4 型角向基函數(shù)((1b)式)中η,λ, ζ 取值Table 2. η, λ, and ζ of 10 × 6 G4 type angular basis functions (Eq. (1b)).

原子在一個方向上的受力可以寫為總能量對該方向位移偏導(dǎo)的相反數(shù), 表達(dá)式為[46]

其總體神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)如圖1 所示.

對于隱藏層神經(jīng)元, 激活函數(shù)選擇的是線性彎曲的雙曲正切函數(shù), 其表達(dá)式為[47]

對于每一種隱藏層網(wǎng)絡(luò)結(jié)構(gòu)所對應(yīng)的神經(jīng)網(wǎng)絡(luò), 分別隨機(jī)選取90%數(shù)據(jù)集(即13680 個)作為訓(xùn)練集進(jìn)行最多2000 代誤差反向傳播學(xué)習(xí)[48], 對神經(jīng)網(wǎng)絡(luò)中的連接權(quán)重和偏置進(jìn)行更新優(yōu)化, 以期使誤差達(dá)到最小, 其余10%的數(shù)據(jù)集(即1519 個)作為測試集來評估訓(xùn)練得到的勢函數(shù)效果.

圖1 高維神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意圖, C代表原子坐標(biāo), G 代表基底函數(shù), H 代表隱藏層神經(jīng)元, E 代表原子能量, 下角標(biāo)1, 2, ···, n 為原子序號, ET 為體系總能量Fig. 1. Structure of high-dimensional neural network. C, G,H, and E represent coordinates of atom , basis functions,hidden layer neurons, and energy of atom , respectively.Subscripts, 1, 2, ···, n are the serial numbers of atoms, and ET is the total energy of the system.

進(jìn)一步, 使用機(jī)器學(xué)習(xí)所得到的勢函數(shù)對2 ×2 × 2 晶胞的a相RDX 晶體進(jìn)行分子動力學(xué)模擬計算. 模擬計算應(yīng)用LAMMPS (large-scale atomic/molecular massively parallel simulator)[49]軟件包,使用麥克斯韋-玻爾茲曼分布進(jìn)行速度初始化, 時間步長設(shè)置為0.2 fs, 在300 K 溫度及13194.4 ?3體積(對應(yīng)密度為1.79 g/cm–3)條件下進(jìn)行正則系綜(NVT)模擬, 模擬時長為20 ps.

3 結(jié)果與討論

3.1 勢函數(shù)訓(xùn)練

圖2 為a相RDX 晶體結(jié)構(gòu)和4 種常見RDX分子構(gòu)型示意圖.a相RDX 是在常溫常壓下穩(wěn)定存在的晶型, 最易在實驗上觀測到, 其單胞中含有8 個RDX 分子, 空間群為Pbca[2,3,7]. 同時, 根據(jù)硝基與三嗪環(huán)軸線夾角的不同, RDX 分子有AAA,AAE, AEE 和EEE 四種常見構(gòu)型, 如圖2(b)所示, 其中, 在a相中RDX 分子為AAE 構(gòu)型[10].

圖2 (a) a 相RDX 晶體結(jié)構(gòu)示意圖; (b) 4 種常見RDX 分子構(gòu)型示意圖[10], 白球代表氫原子, 灰球代表碳原子, 藍(lán)球代表氮原子, 紅球代表氧原子Fig. 2. (a) Structure of a-RDX crystal; (b) structures of four usual types of RDX molecules.[10] The white, grey, blue, and red balls represent hydrogen, carbon, nitrogen, and oxygen atoms, respectively.

圖3 (a) 9 種隱藏層結(jié)構(gòu)在400 代之后測試集最低RMSE 隨迭代步數(shù)變化示意圖; (b) 30-30-10 隱藏層網(wǎng)絡(luò)結(jié)構(gòu)在1500 代之后訓(xùn)練集和測試集RMSE 隨迭代步數(shù)變化示意圖Fig. 3. (a) Diagram of test set lowest RMSEs variation along with training iteratons of nine types hidden layer neural structures after 400 iterations; (b) diagram of training and test sets RMSEs variation along with training iteratons of 30-30-10 hidden layer neural structures after 1500 iterations.

目前所確知的RDX 晶型只有a,β,γ和ε相4 種, 而機(jī)器學(xué)習(xí)尤其是神經(jīng)網(wǎng)絡(luò)模型需要大量的數(shù)據(jù)來訓(xùn)練. 因此, 首先進(jìn)行RDX 晶體結(jié)構(gòu)搜索,隨后對搜索得到的結(jié)構(gòu)使用改變晶格參數(shù)和斜切等方法產(chǎn)生大量新結(jié)構(gòu)并計算其能量和原子受力,共獲取15199 個結(jié)構(gòu)用于訓(xùn)練勢函數(shù)以及對訓(xùn)練效果進(jìn)行評估. 采用Behler 等[29]提出的HDNN方法, 如圖1 所示, 即對每一個原子采用(1a)式和(1b)式兩種共72 個基函數(shù), 通過周圍原子環(huán)境構(gòu)建輸入層神經(jīng)元, 其參數(shù)選取如表1 和表2 所列. 隱藏層分別采用9 種不同的網(wǎng)絡(luò)結(jié)構(gòu). 輸出神經(jīng)元為該原子的能量, 所有原子能量加和為體系總能量.

為了找到最適于描述RDX 晶體勢函數(shù)的網(wǎng)絡(luò)結(jié)構(gòu), 對于9 種不同的隱藏層結(jié)構(gòu)所對應(yīng)的神經(jīng)網(wǎng)絡(luò), 進(jìn)行最多2000 代的誤差反向傳播學(xué)習(xí)[48],并使用測試集的均方根誤差(RMSE)對訓(xùn)練效果進(jìn)行評估. 圖3(a)為400 代后9 種網(wǎng)絡(luò)結(jié)構(gòu)的測集的最小RMSE 隨迭代步數(shù)的變化圖. 由圖3(a)可以看出, 30-30-10 隱層結(jié)構(gòu)在1268 代后具有最低的RMSE, 說明30-30-10 隱層網(wǎng)絡(luò)結(jié)構(gòu)訓(xùn)練得到的勢函數(shù)對能量的計算相比于第一性原理誤差最小, 更加適合于對RDX 晶體的描述.

圖3(b)是1500 代后30-30-10 隱藏層網(wǎng)絡(luò)結(jié)構(gòu)的訓(xùn)練集和測試集RMSE 隨迭代步數(shù)變化示意圖. 由圖3(b)可以看出, 這一階段訓(xùn)練集的RMSE隨迭代步數(shù)近似線性下降, 而測試集的RMSE 隨迭代步數(shù)增加無明顯變化趨勢, 其最低點出現(xiàn)在1847 代, RMSE 為59.2 meV/atom, 即圖中兩條藍(lán)色虛線交點處. 在1847 代后, 測試集的RMSE 上升,而訓(xùn)練集的RMSE 繼續(xù)下降, 說明此時對勢函數(shù)的訓(xùn)練出現(xiàn)了過擬合. 因此, 選取第1847 代訓(xùn)練后所得到的勢函數(shù)進(jìn)行后續(xù)的計算, 此時訓(xùn)練集的平均絕對誤差(MAE)為29.2 meV/atom, RMSE 為47.1 meV/atom, 測試集MAE 為35.1 meV/atom,RMSE 為59.2 meV/atom.

圖4 (a) 訓(xùn)練集(黑色叉)和測試集(紅色十字)所有結(jié)構(gòu)第一性原理計算形成能和機(jī)器學(xué)習(xí)計算形成能對應(yīng)關(guān)系圖; (b) 訓(xùn)練集(黑色叉)和測試集(紅色十字)所有結(jié)構(gòu)第一性原理計算原子受力和機(jī)器學(xué)習(xí)計算原子受力對應(yīng)關(guān)系示意圖Fig. 4. (a) Correlation of machine learning binding energies with the corresponding ab initio reference energies for all structures in the training (black skew crosses) and testing (red crosses) sets; (b) correlation of machine learning atomic forces with the corresponding ab initio reference forces for all structures in the training (black skew crosses) and testing (red crosses) sets.

使用此勢函數(shù)對包括訓(xùn)練集和測試集所有的15199 個結(jié)構(gòu)進(jìn)行能量和晶胞中原子受力的計算,其結(jié)果與第一性原理計算結(jié)果對比如圖4 所示,MAE 與RMSE 列在表3 中. 圖4(a)為第一性原理計算結(jié)合能和機(jī)器學(xué)習(xí)勢函數(shù)計算形成能的對應(yīng)關(guān)系圖, 圖中藍(lán)色虛線代表第一性原理計算能量和機(jī)器學(xué)習(xí)計算能量相等, 數(shù)據(jù)點離藍(lán)色虛線越近說明計算誤差越小. 由圖4(a)可以看出, 除個別結(jié)構(gòu)外, 絕大多數(shù)數(shù)據(jù)點都在藍(lán)色虛線附近, 說明該勢函數(shù)可以很好地重復(fù)出第一性原理計算所得到的結(jié)合能, 而且訓(xùn)練集和測試集在虛線兩側(cè)分布均勻, 比例大致相當(dāng), 說明數(shù)據(jù)集的選取具有代表性.圖4(b)為第一性原理計算晶胞中原子受力和機(jī)器學(xué)習(xí)計算原子受力的對應(yīng)關(guān)系示意圖. 因每個原子有x,y,z三個方向的受力, 每個結(jié)構(gòu)晶胞中有168 個原子, 計算所得的原子受力過多, 無法在一張圖中做出, 這里對每一個結(jié)構(gòu)隨機(jī)選取了一個原子在一個方向上的受力, 做出此示意圖. 由圖4(b)可以看出, 圖中所有數(shù)據(jù)點都在藍(lán)色虛線附近,而且?guī)讉€受力極大(> 200 eV/?)的原子, 機(jī)器學(xué)習(xí)勢函數(shù)的計算結(jié)果也和第一性原理計算結(jié)果相近. 如表3 所列, 與第一性原理計算相比較,訓(xùn)練所得到的勢函數(shù)不僅能較好得計算晶體能量信息, 同時也能得到可信的原子受力信息, 說明該勢函數(shù)適用于對不同晶型的RDX 晶體能量與受力描述.

為驗證勢函數(shù)對已知結(jié)構(gòu)計算的可靠性, 對于a,β,γ和ε四種晶型, 使用VASP[42]軟件包,在0 K溫度下分別在標(biāo)準(zhǔn)大氣壓到6 GPa 以1 GPa 為間隔, 共7 種壓強(qiáng)條件下進(jìn)行結(jié)構(gòu)優(yōu)化并對其能量和原子受力進(jìn)行計算, 進(jìn)而使用機(jī)器學(xué)習(xí)所得的勢函數(shù)對其結(jié)合能和原子受力進(jìn)行計算, 并與之比較, 其結(jié)合能對比結(jié)果如圖5 所示. 可以看出, 能量誤差最大的結(jié)構(gòu)為γ相在1 和2 GPa 條件下的結(jié)構(gòu), 此時機(jī)器學(xué)習(xí)勢函數(shù)所計里算結(jié)合能較第一性原理所計算結(jié)合能低20.1 meV/atom,其余26 種結(jié)構(gòu)的絕對誤差均小于20 meV/atom.對于這28 種結(jié)構(gòu), 結(jié)合能的MAE 為8.8 meV/atom,RMSE 為10.0 meV/atom, 原子受力的MAE 為0.88 eV/?, RMSE 為1.11 eV/?, 均遠(yuǎn)小于訓(xùn)練集和測試集的結(jié)合能和原子受力的MAE 和RMSE(如表3 所列), 說明對于已知穩(wěn)定結(jié)構(gòu), 機(jī)器學(xué)習(xí)勢函數(shù)的計算效果優(yōu)于對于隨機(jī)不穩(wěn)定結(jié)構(gòu)的計算, 因此可將該勢函數(shù)用于對更多未知的穩(wěn)定結(jié)構(gòu)的預(yù)測.

表3 訓(xùn)練集與測試集機(jī)器學(xué)習(xí)計算形成能和原子受力與第一性原理計算比較MAE 和RMSETable 3. MAE and RMSE of machine learning binding energies and atomic forces corresponding ab initio reference energies and forces in the training and test sets.

圖5 0 K 下4 種已知RDX 晶型在標(biāo)準(zhǔn)大氣壓到6 GPa 壓強(qiáng)下機(jī)器學(xué)習(xí)計算結(jié)合能的誤差Fig. 5. Errors of machine learning binding energies of four known RDX crystals from 1 atm to 6 GPa at 0 K.

圖6 (a) 4 種RDX 晶型結(jié)構(gòu)示意圖; (b) 7—10 GPa 下4 種晶型第一性原理計算(黑色方塊)和機(jī)器學(xué)習(xí)勢函數(shù)計算(紅色圓圈)的焓值Fig. 6. (a) Structures of four RDX crystals; (b) enthalpies from 7 to 10 GPa calculated by ab initio (black block ) and machine learning potential (red circle).

因存在實驗上觀測到而理論上未確定具體結(jié)構(gòu)的高壓RDX 晶型, 我們在搜索到的248 種構(gòu)型中選取了4 種能量較低的晶型, 其結(jié)構(gòu)如圖6(a)所示. 對這4 種晶型分別在7—10 GPa 四種壓強(qiáng)條件下優(yōu)化, 并使用VASP[42]軟件包和訓(xùn)練得到的勢函數(shù)進(jìn)行焓值計算, 其結(jié)果如圖6(b)所示. 由圖6(b)可以看出, 在4 種壓強(qiáng)條件下, 機(jī)器學(xué)習(xí)勢函數(shù)均能給出與第一性原理計算相近的能量排序,焓值的RMSE 為13.7 meV/atom, 與已知結(jié)構(gòu)能量的RMSE 相近, 進(jìn)一步證明了該勢函數(shù)對RDX分子晶體能量描述的準(zhǔn)確性.

3.2 分子動力學(xué)模擬

為了考察該勢函數(shù)的泛用性, 使用機(jī)器學(xué)習(xí)所得勢函數(shù)對a相RDX 晶體的2 × 2 × 2 晶胞進(jìn)行分子動力學(xué)模擬. 模擬的時間步長設(shè)置為0.2 fs,總時長為20 ps, 在300 K 溫度及不同體積條件下進(jìn)行NVT 系綜模擬. 在13194.4 ?3體積(對應(yīng)密度為1.79 g/cm3)下,其溫度和壓強(qiáng)隨時間變化如圖7 所示.由圖7(a)可知,在1 ps之前,體系溫度較高, 而在1.0—2.5 ps之間, 溫度的波動較大, 2.5 ps之后溫度在300 K附近平穩(wěn)波動.由圖7(b)可知,在2.5 ps之前,壓強(qiáng)的波動幅度較大,在2.5 ps之后壓強(qiáng)平穩(wěn)波動,其平均壓強(qiáng)為4.63 GPa.由此說明,體系在2.5 ps之前為體系弛豫過程,結(jié)構(gòu)發(fā)生一定程度上的變形和重構(gòu),在2.5 ps之后,結(jié)構(gòu)穩(wěn)定,原子均在分子的平衡位置附近做規(guī)律運動.上述結(jié)果表明,該勢函數(shù)可以較好地模擬分子內(nèi)和分子間相互作用, 初步了證明其在分子模擬中的適用性.

圖7 (a)a 相RDX 2×2×2晶胞在NVT系綜下分子動力學(xué)模擬溫度隨時間變化圖;(b)壓強(qiáng)隨時間變化圖Fig.7.Variations in time of the tem perature (a)and pressure(b)in the NVT ensem ble for 2×2×2 a-RDX crystal.

4 結(jié)論

本文采用HDNN[29]方法,對15199個RDX分子晶體結(jié)構(gòu)數(shù)據(jù)集進(jìn)行機(jī)器學(xué)習(xí)訓(xùn)練并得到相應(yīng)的勢函數(shù).通過對9種不同隱藏層結(jié)構(gòu)的神經(jīng)網(wǎng)絡(luò)進(jìn)行最多2000代誤差反向傳播訓(xùn)練[48],并對其訓(xùn)練效果進(jìn)行評估,選取測試集結(jié)合能RMSE最低的30-30-10隱藏層結(jié)構(gòu)神經(jīng)網(wǎng)絡(luò)模型的訓(xùn)練所得勢函數(shù).相較第一性原理計算結(jié)果,機(jī)器學(xué)習(xí)擬合的勢函數(shù)所計算的結(jié)合能和原子受力的誤差如表3所列,均能很好地重復(fù)出第一性原理的計算結(jié)果.使用該勢函數(shù)對4種已知晶型的RDX晶體在不同壓強(qiáng)條件下計算結(jié)合能并與第一性原理計算結(jié)果相對比, 其結(jié)合能RMSE僅為10.0 meV/atom,原子受力RMSE為1.11 eV/?,說明該勢函數(shù)對于已知穩(wěn)定結(jié)構(gòu)的計算更加合理.

進(jìn)一步,使用訓(xùn)練得到的勢函數(shù)對a相RDX晶體的2×2×2晶胞在NVT系綜下進(jìn)行分子動力學(xué)模擬,發(fā)現(xiàn)其溫度和壓強(qiáng)在2.5 ps后波動穩(wěn)定,證明該勢函數(shù)對分子內(nèi)和分子間相互作用均能較好地描述.本文針對含能材料RDX 晶體進(jìn)行機(jī)器學(xué)習(xí)勢函數(shù)的構(gòu)建,為今后分子晶體方向的機(jī)器學(xué)習(xí)研究提供幫助和啟示.

感謝北京計算科學(xué)研究中心管鵬飛研究員、北京應(yīng)用物理與計算數(shù)學(xué)研究所宋華杰研究員的幫助,以及大連理工大學(xué)超級計算中心提供的計算支持.