基于機器學習構建的環三亞甲基三硝胺晶體勢*

王鵬舉 范俊宇 蘇艷 趙紀軍

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

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

1 引 言

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

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

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

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

2 研究方法

2.1 數據集獲取

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

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

2.2 機器學習模型構建

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

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

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

(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 種隱藏層結構進行神經網絡構建.輸出層只包含一個神經元, 代表該原子的能量, 所有原子能量加和為體系總能量, 其表達式為[46]

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

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

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

其總體神經網絡結構如圖1 所示.

對于隱藏層神經元, 激活函數選擇的是線性彎曲的雙曲正切函數, 其表達式為[47]

對于每一種隱藏層網絡結構所對應的神經網絡, 分別隨機選取90%數據集(即13680 個)作為訓練集進行最多2000 代誤差反向傳播學習[48], 對神經網絡中的連接權重和偏置進行更新優化, 以期使誤差達到最小, 其余10%的數據集(即1519 個)作為測試集來評估訓練得到的勢函數效果.

圖1 高維神經網絡結構示意圖, C代表原子坐標, G 代表基底函數, H 代表隱藏層神經元, E 代表原子能量, 下角標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.

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

3 結果與討論

3.1 勢函數訓練

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

圖2 (a) a 相RDX 晶體結構示意圖; (b) 4 種常見RDX 分子構型示意圖[10], 白球代表氫原子, 灰球代表碳原子, 藍球代表氮原子, 紅球代表氧原子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 種隱藏層結構在400 代之后測試集最低RMSE 隨迭代步數變化示意圖; (b) 30-30-10 隱藏層網絡結構在1500 代之后訓練集和測試集RMSE 隨迭代步數變化示意圖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 種, 而機器學習尤其是神經網絡模型需要大量的數據來訓練. 因此, 首先進行RDX 晶體結構搜索,隨后對搜索得到的結構使用改變晶格參數和斜切等方法產生大量新結構并計算其能量和原子受力,共獲取15199 個結構用于訓練勢函數以及對訓練效果進行評估. 采用Behler 等[29]提出的HDNN方法, 如圖1 所示, 即對每一個原子采用(1a)式和(1b)式兩種共72 個基函數, 通過周圍原子環境構建輸入層神經元, 其參數選取如表1 和表2 所列. 隱藏層分別采用9 種不同的網絡結構. 輸出神經元為該原子的能量, 所有原子能量加和為體系總能量.

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

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

圖4 (a) 訓練集(黑色叉)和測試集(紅色十字)所有結構第一性原理計算形成能和機器學習計算形成能對應關系圖; (b) 訓練集(黑色叉)和測試集(紅色十字)所有結構第一性原理計算原子受力和機器學習計算原子受力對應關系示意圖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.

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

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

表3 訓練集與測試集機器學習計算形成能和原子受力與第一性原理計算比較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 晶型在標準大氣壓到6 GPa 壓強下機器學習計算結合能的誤差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 晶型結構示意圖; (b) 7—10 GPa 下4 種晶型第一性原理計算(黑色方塊)和機器學習勢函數計算(紅色圓圈)的焓值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).

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

3.2 分子動力學模擬

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

圖7 (a)a 相RDX 2×2×2晶胞在NVT系綜下分子動力學模擬溫度隨時間變化圖;(b)壓強隨時間變化圖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 結論

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

進一步,使用訓練得到的勢函數對a相RDX晶體的2×2×2晶胞在NVT系綜下進行分子動力學模擬,發現其溫度和壓強在2.5 ps后波動穩定,證明該勢函數對分子內和分子間相互作用均能較好地描述.本文針對含能材料RDX 晶體進行機器學習勢函數的構建,為今后分子晶體方向的機器學習研究提供幫助和啟示.

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