含能材料感度評(píng)估的空間位阻指數(shù)算法和程序設(shè)計(jì)

崔 濤，劉知涵，謝 瑋，薛向貴

（1. 上海大學(xué)理學(xué)院化學(xué)系，上海 200444；2. 上海大學(xué)材料基因組工程研究院，上海 200444；3. 中國(guó)工程物理研究院化工材料研究所，四川 綿陽(yáng) 621999）

1 引言

含能材料的感度研究在近年愈發(fā)受到關(guān)注。以實(shí)驗(yàn)方式得到含能材料的撞擊感度需要合成含能材料大單晶因而比較困難，如何通過(guò)非實(shí)驗(yàn)的方式比較不同含能材料的撞擊感度一直是此領(lǐng)域的研究熱點(diǎn)。在過(guò)去對(duì)含能材料撞擊感度的研究中，Menikoff等［1］及Dick 等［2］報(bào)道了奧克托今（HMX）的β 相在撞擊載荷下的塑性變形與滑移運(yùn)動(dòng)直接相關(guān)，熱力學(xué)響應(yīng)和撞擊感度具有明顯的各向異性，晶體內(nèi)滑移面兩側(cè)的空間位阻對(duì)各向異性的感度起主要作用，且敏感晶面的空間位阻大。 Zyben［3］使用基于ReaxFF 反應(yīng)力場(chǎng)的分子動(dòng)力學(xué)模擬研究了太安（PETN）單晶在壓縮-剪切載荷下的物理化學(xué)響應(yīng)，計(jì)算得到的熱力學(xué)響應(yīng)和撞擊感度的各向異性與撞擊實(shí)驗(yàn)結(jié)果一致。Swell 等［4］報(bào)道了HMX 的塑性變形以沿著某些滑移系的位錯(cuò)運(yùn)動(dòng)為主。同年An 等［5］使用壓縮-剪切反應(yīng)動(dòng)力學(xué)（compress-and-shear reactive dynamics，CS-RD）模擬研究了奧克托今（RDX）晶體的撞擊感度的各向異性，所得計(jì)算結(jié)果與已有實(shí)驗(yàn)具有較高的一致性。但是，由于缺乏合適的力場(chǎng)等原因，苯并氧化呋咱（BTF）和梯恩梯（TNT）等仍未有公開的CS-RD 記錄報(bào)道，因此使用第一性原理和其他方法研究含能材料撞擊感度仍存在較大的發(fā)展空間。

除了分子動(dòng)力學(xué)模擬以外，An 等［5］還基于含能材料中分子重疊程度定義了空間位阻指數(shù)（Steric Hindrance Index，SHI）。使用SHI 定量評(píng)估和預(yù)測(cè)含能材料撞擊感度的優(yōu)勢(shì)在于只需晶體結(jié)構(gòu)作為輸入，不考慮原子馳豫和化學(xué)反應(yīng)，無(wú)須做代價(jià)昂貴的分子動(dòng)力學(xué)模擬，因而高效且獨(dú)立于力場(chǎng)。盡管An等［5］定義了SHI，但受制于自動(dòng)計(jì)算算法及其程序?qū)崿F(xiàn)的缺乏，僅手動(dòng)計(jì)算了RDX 晶體的（100）/｛-110｝<110>和（110）/｛010｝<100>兩個(gè)滑移系的SHI。本文在An 的研究基礎(chǔ)上，開發(fā)出具有材料普適性的SHI 自動(dòng)計(jì)算算法，并使用Python3 進(jìn)行程序?qū)崿F(xiàn)。用戶在程序中僅需輸入正交晶系的含能材料的cif 文件，即可計(jì)算任意一個(gè)滑移系的SHI，相比于涉及力場(chǎng)的CS-RD 計(jì)算方法極大的簡(jiǎn)化了計(jì)算過(guò)程。

2 設(shè)計(jì)目標(biāo)與算法

2.1 空間位阻指數(shù)的定義

如圖1 所示，含能材料受到某一方向的撞擊后，在沿撞擊平面（圖1 示例中為（110）晶面）一定夾角θ（圖1 示例中θ=45°）產(chǎn)生剪切滑移面（也叫剪切帶）。通過(guò)沿撞擊方向等軸縮放晶體來(lái)壓縮分子。將每個(gè)分子中的各個(gè)原子視作球體，將壓縮后分子中的各個(gè)原子完整幾何投影至剪切滑移面（圖1 示例中為｛010｝面）的垂直平面上（圖1 示例中為底面）。對(duì)指定的滑移系建立新的空間直角坐標(biāo)系并旋轉(zhuǎn)晶胞以處理任意撞擊方向和滑移系，將旋轉(zhuǎn)后的晶胞內(nèi)分子根據(jù)質(zhì)心的x 坐標(biāo)進(jìn)行分層，如圖2。如果相鄰平面上的分子在投影后重疊，則在投影平面上的位置給出一個(gè)陰影點(diǎn)，即重疊數(shù)加一。將重疊數(shù)加權(quán)的重疊區(qū)域面積相加，相對(duì)于投影橫截面的比率即被定義為空間位阻指數(shù)?？臻g位阻指數(shù)可以繪制二維等值線圖，其中重疊區(qū)域表示在剪切過(guò)程中可能發(fā)生空間位阻的區(qū)域。

圖1 分子在剪切形變過(guò)程中的空間位阻示意圖Fig.1 Schematic diagram of steric hindrance of molecules during shear deformation

圖2 相鄰層分子在投影區(qū)域的重疊部分Fig.2 Overlapping part of adjacent layer molecules on the projection plane

2.2 壓縮

為模擬真實(shí)情況下含能材料受到撞擊時(shí)所產(chǎn)生的形變，對(duì)晶胞沿著撞擊面進(jìn)行壓縮。假設(shè)分子為剛性即在撞擊下不變形，只是壓縮分子質(zhì)心。由于分子晶體中分子間的范德華作用較弱，當(dāng)受到撞擊時(shí)，主要是分子間距離變化，由于分子本身相對(duì)剛性較高其變形可以近似地忽略。具體步驟如下：

設(shè)壓縮的壓縮比率為r（r>0），晶胞的晶格常數(shù)為a，b，c，撞擊面的晶面指數(shù)為（jkl）。參考實(shí)際情況j，k，l≥0。注意到由于晶面指數(shù)的定義，j，k，l 不能同時(shí)為零。任取晶胞中一原子A，進(jìn)行壓縮前坐標(biāo)為（x，y，z），進(jìn)行壓縮后坐標(biāo)為（x'，y'，z'），由撞擊面晶面指數(shù)j、k、l 及壓縮比r 可得原子各軸壓縮比分別為ra、rb、rc，則有：

2.3 旋轉(zhuǎn)

由于每一個(gè)滑移系的滑移面和滑移方向上的空間位阻存在差別，為了保證算法的一致性需要建立新的空間直角坐標(biāo)系，使得滑移方向?yàn)樾伦鴺?biāo)系的X 軸正方向，滑移面為新坐標(biāo)系的XOZ 平面，如圖3 所示。

在對(duì)晶胞進(jìn)行壓縮后需要對(duì)晶胞進(jìn)行旋轉(zhuǎn)。該步輸入的信息為滑移面｛pqr｝以及滑移方向。同樣地p，q，r 和u，v，w 均不能同時(shí)為零。滑移面的法線向量為｛pqr｝T，而滑移方向在滑移面上，所以有：

壓縮后坐標(biāo)系的三條坐標(biāo)軸為X、Y、Z 軸，旋轉(zhuǎn)后坐標(biāo)系的三條坐標(biāo)軸為X'、Y'、Z'軸。則壓縮坐標(biāo)系原子A 坐標(biāo)仍記為為（x，y，z），旋轉(zhuǎn)后的坐標(biāo)記為（x'，y'，z'），旋轉(zhuǎn)之前有：

旋轉(zhuǎn)后，根據(jù)原滑移系與旋轉(zhuǎn)后的滑移系｛010｝<100>各軸之間的夾角，計(jì)算出在旋轉(zhuǎn)后的坐標(biāo)系下晶胞內(nèi)各原子在各軸最大及最小坐標(biāo)：x'min，x'max，y'min，y'max，z'min，z'max。

圖3 原坐標(biāo)系與新坐標(biāo)系示意圖Fig.3 Schematic diagram of original coordinate system and new coordinate system

2.4 分層

在實(shí)際的剪切形變過(guò)程中，在滑移方向上相鄰的分子之間的碰撞是影響含能材料感度的重要因素。旋轉(zhuǎn)坐標(biāo)軸后，按照垂直于剪切滑移面的方向以分子質(zhì)心的x 坐標(biāo)對(duì)分子進(jìn)行分層，并且假設(shè)歸屬于同一層內(nèi)的分子無(wú)法碰撞。分層時(shí)如兩相鄰分子質(zhì)心x 坐標(biāo)不大于質(zhì)心x 坐標(biāo)相差最大的兩分子質(zhì)心x 坐標(biāo)差值的5%，則將兩分子歸屬同一平面內(nèi)。

2.5 投影

在剪切形變過(guò)程中，晶胞內(nèi)的分子沿著滑移方向，即X 軸正方向運(yùn)動(dòng)并發(fā)生碰撞，為了定量計(jì)算不同分子層間碰撞的劇烈程度，將分子層投影到Y(jié)OZ平面上。投影區(qū)域設(shè)為矩形［y'min，y'max］*［z'min，z'max］。將投影區(qū)域細(xì)分為數(shù)量ny*nz個(gè)最小單元格組成的矩形。則Y 軸上步長(zhǎng)為（y'max-y'min）/ny，Z軸上步長(zhǎng)為（z'max-z'min）/nz。投 影 區(qū) 域 被 劃 分 為ny*nz個(gè) 邊 長(zhǎng) 為（y'max-y'min）/ny，（z'max-z'min）/nz的小矩形。記所有的小矩形構(gòu)成的集合為R=｛rij｝，i=1，2…，ny，j=1，2，…，nz。

分子層是由若干個(gè)分子組成的，而每個(gè)分子包含m 個(gè)原子，每個(gè)原子在空間中是一個(gè)以該原子坐標(biāo)為球心，原子半徑為半徑的球。每個(gè)原子投影到Y(jié)OZ 平面上之后實(shí)質(zhì)上是一個(gè)圓。由于分子層的不同原子在投影區(qū)域內(nèi)所形成的圓是相交的且有可能多個(gè)圓相交導(dǎo)致計(jì)算分子層所形成的面積較為復(fù)雜，因此使用ny*nz矩形的面積之和對(duì)分子層的投影面積進(jìn)行近似。

任取某一分子層Lk，對(duì)于任意分子Mj屬于Lk，任意原子A 屬于Mj，則A 屬于Lk。設(shè)Lk在YOZ 平面上的投影所對(duì)應(yīng)的小矩形的集合為Sk。選取如下方法求Sk：

rij是劃分的最小方格的中心點(diǎn)，對(duì)于在分子層Lk內(nèi)的原子A，若劃分的最小格子中心點(diǎn)在A 所投影的YOZ 平 面 上 的 圓 內(nèi)，則 將rij計(jì) 入Sk.。如 圖4 所 示，r_ij和r_（i+1）j不屬于Sk，r_i，（j+1），r_（i+1）（j+1）屬于Sk。Sk集合對(duì)應(yīng)的中心點(diǎn)數(shù)即對(duì)應(yīng)的最小方格數(shù)與投影區(qū)域全部最小方格數(shù)量之比即為SHI。

圖4 投影方法示意圖Fig.4 Schematic diagram of projection method

設(shè)晶胞經(jīng)過(guò)壓縮與旋轉(zhuǎn)并按照分子質(zhì)心的x 坐標(biāo)排序后形成了t 個(gè)分子層，對(duì)每個(gè)分子層進(jìn)行投影后，得到每一個(gè)分子層中被包含在原子投影中的最小方格中心格點(diǎn)的集合Si。隨后計(jì)算相鄰分子層重疊的中心點(diǎn)，即為計(jì)算投影的重疊面積。設(shè)集合Gi是集合Si與集合Si+1，的交集，其中i=1，2，…，t-1。另外需要指出的是，晶胞在晶體內(nèi)呈周期性排列，所以需要計(jì)算St與S1的并集，記為Gt。

2.6 滑移系的自動(dòng)搜尋

考慮某一撞擊面（pk，qk，rk）中pk，qk，rk不同時(shí)為0。希望通過(guò)計(jì)算機(jī)自動(dòng)搜尋np個(gè)滑移面，并在每個(gè)滑移面上搜尋nd個(gè)滑移方向作為備選滑移系。事實(shí)上，對(duì)于一個(gè)晶胞，低指數(shù)面方向上原子排布更密集，通過(guò)枚舉法達(dá)到自動(dòng)搜尋滑移系的目的。

3 計(jì)算結(jié)果與分析

3.1 投影算法的評(píng)判

構(gòu)造如下的簡(jiǎn)單晶胞，晶格常數(shù)為a=b=c=10 ?，晶胞內(nèi)有兩個(gè)分子M1，M2。每個(gè)分子由四個(gè)原子組成，設(shè)為Ai，i=1，2，…，8，其中，A1，A2，A3，A4屬于M1，A5，A6，A7，A8屬于M2，原子半徑均為1。8 個(gè)原子在晶胞內(nèi)的坐標(biāo)分別為A1（1，1，1），A2（1，9，1），A3（1，1，9），A4（1，9，9），A5（9，1，1），A6（9，9，1），A7（9，1，9），A8（9，9，9）。不 對(duì) 晶 胞 進(jìn) 行 壓 縮，同 時(shí) 滑 移 面 為｛010｝，滑移方向?yàn)?100>，即不旋轉(zhuǎn)晶胞，直接進(jìn)行分層與投影。該晶胞分成兩層，投影區(qū)域?yàn)椋?，10］*［0，10］，每個(gè)分子層在投影區(qū)域內(nèi)的投影如圖5所示。

圖5 構(gòu)造簡(jiǎn)單晶胞分子層投影示意圖Fig.5 Schematic diagram of molecular layer projection in a single unit cell

該晶胞的兩分子層在投影區(qū)域內(nèi)的重疊部分面積為：

在不同的劃分精度ny，nz，投影算法的空間位阻指數(shù)計(jì)算結(jié)果見(jiàn)表1。

在ny=nz=100 時(shí)，SHI 為0.1264，誤 差 為0.6%。ny=nz=500 時(shí)，該方法的誤差為0.08%，在ny，nz不太大時(shí)，近似計(jì)算方式就可以得到較良好的模擬結(jié)果。

3.2 RDX 晶體的SHI 與CS?RD 模擬結(jié)果比較

選取α-RDX 單晶的晶格常數(shù)為a=13.182 ?，b=11.574 ?，c=10.709 ?［6］。

垂直于下列五個(gè)低指數(shù)平面的撞擊面：（100），（210），（111），（110）和（120）。參照An 等［5］的實(shí)驗(yàn)，對(duì)5 個(gè)撞擊面選擇計(jì)算主要滑移系：（A）（210）/｛120｝<-210>；（B）（100）/｛-110｝<110>；（C）（111）/｛021｝<100>；（D）（120）/｛010｝<100>；（E）（110）/｛010｝<100>，同時(shí)計(jì)算非主要滑移系的SHI 作為對(duì)照。對(duì)RDX 晶體 的CS-RD 實(shí) 驗(yàn) 表 明，（A），（B）為 敏 感 滑 移 系，而（C），（D），（E）為非敏感滑移系。由于非主要滑移系在剪切形變的過(guò)程中較少產(chǎn)生，因而非主要滑移系屬于敏感或非敏感晶系未被納入討論。

從表2 對(duì)比可得敏感滑移系1，7 的空間位阻指數(shù)顯著大于非敏感滑移系9，16，20 的空間位阻指數(shù)。

為 與An 等［5］報(bào) 道 的RDX 的CS-RD 對(duì) 比，同 樣 選定壓縮比（r）為0.1 和0.2。壓縮比增加則系內(nèi)晶體分子碰撞后溫度與NO2含量都會(huì)相應(yīng)增加。對(duì)滑移系的SHI 計(jì)算結(jié)果分別見(jiàn)表2 和表3。

從表3 可得敏感滑移系1，5 的空間位阻指數(shù)顯著大于非敏感滑移系11，17，21 的空間位阻指數(shù)。

3.3 PETN 晶體的SHI 與CS?RD 模擬結(jié)果比較

所研究PETN 單晶的晶格常數(shù)為a=13.29 ?，b=13.49 ?，c=6.83 ?［7］。

為 與Zyben 等［3］報(bào) 道 的PETN 的CS-RD 對(duì) 比，選定壓縮比ratio=0.1。

從表4 可得敏感滑移系1、3 的空間位阻指數(shù)顯著大于中等敏感以及非敏感滑移系4、7、8 的空間位阻指數(shù)。

3.4 β?HMX 的SHI 與CS?RD 模擬結(jié)果比較

所研究β-HMX 單晶的晶格常數(shù)為a=6.5209 ?，b=10.7610 ?，c=7.3063 ?［8］。

為 與ZHOU 等［9］報(bào) 道 的β-HMX 的CS-RD 對(duì) 比，選定壓縮比ratio=0.17。

對(duì)非正交滑移系的計(jì)算，需要用atomsk 軟件將非正交晶系轉(zhuǎn)換為正交晶系。以單斜晶系中β-HMX 為例，結(jié)果見(jiàn)表5。

非正交晶系SHI 計(jì)算結(jié)果和CS-RD 計(jì)算偏差較大，因此需要詳細(xì)分析對(duì)于非正交轉(zhuǎn)為正交晶系與初始晶系為正交晶系之間的區(qū)別，此外亦可能是在SHI算法中未考慮化學(xué)因素對(duì)碰撞所產(chǎn)生的影響。

3.5 SHI 與撞擊后8 ps 體系溫度以及撞擊后10 ps 時(shí)NO2/RDX 含量的比較

An 等［5］基于RDX 的CS-RD 計(jì)算了每個(gè)滑移系在8 ps 時(shí)的反應(yīng)溫度與10 ps 時(shí)NO2/RDX 含量。這兩個(gè)數(shù)據(jù)可以體現(xiàn)反應(yīng)的劇烈程度，即8 ps 體系溫度及10 ps 時(shí)的NO2/RDX 含量越高滑移系的敏感程度越高。CS-RD 計(jì)算結(jié)果與SHI 比較見(jiàn)表6（r=0.1），其相應(yīng)的圖示見(jiàn)圖6。

如表6 和圖6 顯示，5 個(gè)主要滑移系的SHI 與8 ps時(shí)的溫度，10 ps 時(shí)的NO2/RDX 大體上趨勢(shì)相同，可有效地區(qū)分敏感與非敏感滑移系，敏感滑移系的SHI 顯著大于非敏感化學(xué)系SHI。在r=0.1 時(shí)RDX 的5 個(gè)主要滑移系的SHI 與CS-RD 反應(yīng)體系在8 ps 時(shí)的溫度相關(guān)系數(shù)為0.9160，即SHI 與CS-RD 反應(yīng)體系在8 ps 時(shí)溫度呈強(qiáng)相關(guān)性。SHI 與CS-RD 反應(yīng)體系10 ps 時(shí)NO2/RDX 含量相關(guān)系數(shù)為0.8134，即SHI 與CS-RD 反應(yīng)體系10 ps NO2/RDX 含量呈強(qiáng)相關(guān)性。撞擊后8 ps體系溫度及10 ps NO2/RDX 含量表明敏感滑移系（A）、（B）較非敏感滑移系（C）、（D）、（E）反應(yīng)更劇烈。

3.6 各正交晶系含能材料的SHI 與落錘實(shí)驗(yàn)測(cè)得撞擊感度比較

以表2 的RDX 的22 個(gè)低指數(shù)晶面滑移系為例，計(jì)算PETN，BTF［10］，TNT［11］對(duì)應(yīng)的SHI 及22 個(gè)滑移系的平均SHI 并與落錘實(shí)驗(yàn)測(cè)得的H50（用2.5 kg 的落錘儀對(duì)含能材料進(jìn)行試驗(yàn)，在達(dá)到50%爆炸幾率時(shí)的落高）比較，結(jié)果見(jiàn)表7 及圖7。

表3 r=0.2 時(shí)24 個(gè)滑移系的SHI 與CS-RD 結(jié)果［5］比較Table 3 Comparison of SHI and CS-RD results for 22 slip systems at r=0.2

表4 r=0.1 時(shí)8 個(gè)滑移系的SHI 與CS-RD 結(jié)果［3］比較Table 4 Comparison of SHI and CS-RD results for 8 slip systems with r=0.1

對(duì)不同含能材料的撞擊感度進(jìn)行排序的方法之一是根據(jù)剪切壁壘最小原理［1-5］，按照能量判據(jù)選擇在撞擊作用下最容易發(fā)生滑移的晶面，然后比較主要的滑移面中最低SHI。然而目前對(duì)BTF 和TNT 沒(méi)有合適的力場(chǎng)模擬，沒(méi)有公開的對(duì)TNT 和BTF 關(guān)于CS-RD 的報(bào)道，因此無(wú)法準(zhǔn)確獲知TNT 和BTF 主要滑移面。對(duì)于當(dāng)前無(wú)法獲取主要滑移面的體系，可假設(shè)低指數(shù)密排面為可能的滑移面，低指數(shù)密排方向?yàn)榭赡艿幕品较?。?duì)22 個(gè)低指數(shù)滑移系計(jì)算SHI，然后利用SHI 平均值判斷撞擊感度大小。例如，圖7 為PETN、BTF、RDX、TNT4 種“含能材料0.1”壓縮剪切后22 個(gè)滑移系的SHI；它們的平均SHI 由大至小依次為0.8707，0.7940，0.4228，0.0924（表7）。相應(yīng)的，表8 中給出實(shí)驗(yàn)H50由小至大排列為PETN、BTF、RDX、TNT。四種含能材料的22 個(gè)滑移系的平均SHI 與其落錘實(shí)驗(yàn)對(duì)應(yīng)的H50的相關(guān)系數(shù)為-0.9061，呈強(qiáng)負(fù)相關(guān)性，即正交晶系含能材料的撞擊感度越高，對(duì)應(yīng)的平均SHI越大。由此推斷，這種策略可以適用于對(duì)廣泛的含能分子晶體進(jìn)行相同基準(zhǔn)下的系統(tǒng)比較和搜索優(yōu)化設(shè)計(jì)，對(duì)于推測(cè)正交晶系含能材料的撞擊感度具有較高的參考價(jià)值。

表5 r=0.17 時(shí)8 個(gè)滑移系的SHI 與CS-RD 結(jié)果［9］比較Table 5 Comparison of SHI and CS-RD results for 8 slip systems at r=0.17

表6 5 個(gè)主要滑移系的SHI、8 ps 時(shí)的溫度以及10 ps 時(shí)的NO2/RDX 含量［5］Table 6 SHI of the five main slip systems，temperature at 8 ps，and NO2/RDX content at 10 ps

圖6 RDX 的5 個(gè)主要滑移系的SHI 與碰撞后8 ps 時(shí)的溫度及10 ps 時(shí)NO2/RDX 含量的關(guān)系Fig.6 SHI of the five main slip systems of RDX compared with the temperature at 8 ps and NO2/RDX content at 10 ps after shock

3.7 原子半徑的影響

在空間位阻指數(shù)的計(jì)算過(guò)程中，原子的半徑會(huì)對(duì)SHI 計(jì)算結(jié)果絕對(duì)值造成影響。原子半徑分為軌道半徑、范德華半徑、共價(jià)半徑等。本文所研究的含能材料的晶胞一般為分子晶體，起作用的半徑為范德華半徑。Pauling［20］于20 世 紀(jì)30 年 代 提 出 了 范 德 華 半 徑及一組原子的范德華半徑，但在隨后的近100 年里，化學(xué)家們根據(jù)不同的理論基礎(chǔ)給出了不同的計(jì)算結(jié)果。Bondi［21］于1964 年 根 據(jù) 晶 體 結(jié) 構(gòu) 數(shù) 據(jù)、原 子 的 碰 撞 界面等數(shù)據(jù)得到了一組范德華半徑。 1994 年，Allinger［22］根據(jù)分子力學(xué)得出了孤立狀態(tài)下的原子的范德華半徑，他的計(jì)算結(jié)果較由晶體結(jié)構(gòu)數(shù)據(jù)得出的計(jì)算結(jié)果偏大。胡盛志等［23］在2003 年以晶體中原子的平均體積為出發(fā)點(diǎn)，提出了另一組計(jì)算結(jié)果。此外將原子的共價(jià)半徑也納入考量，以RDX 晶體所包含的H、C、N、O 原子為例，四組范德華半徑及共價(jià)半徑見(jiàn)表9。

表7 RDX、PETN、BTF、TNT 在22 個(gè)低指數(shù)晶面下的SHI 及SHI 均值Table 7 SHI of RDX，PETN，BTF，TNT in 22 low-index crystal plane slip systems and the average SHI values

圖7 對(duì)含能材料10%的壓縮剪切后22 個(gè)滑移系的SHI 對(duì)照Fig.7 SHI comparison of 22 slip systems based on 10% compression shear of energetic materials

表8 PETN、BTF、RDX、TNT 的H50Table 8 H50 of PETN，BTF，RDX and TNT

為研究不同的原子半徑對(duì)空間位阻指數(shù)的影響，計(jì)算了各原子半徑下壓縮比為0.1 時(shí)的SHI，結(jié)果見(jiàn)表10。

由表10 可得原子半徑標(biāo)度的選擇并不會(huì)定性改變對(duì)敏感與非敏感滑移系的區(qū)分。非敏感滑移系（C）、（D）、（E）的空間位阻指數(shù)仍小于敏感滑移系（A）、（B）的空間位阻指數(shù)。

研究中原子半徑采用Allinger 的計(jì)算結(jié)果。

表9 H、C、N、O 的原子半徑［23］Table 9 Atomic radius for H，C，N，O ?

表10 不同原子半徑的SHI 計(jì)算結(jié)果Table 10 SHI with different atomic radius

4 結(jié)論

設(shè)計(jì)了含能材料空間位阻指數(shù)計(jì)算算法，并使用Python 語(yǔ)言開發(fā)了相應(yīng)的計(jì)算機(jī)程序。應(yīng)用于幾種典型的正交晶系含能材料PETN、BTF、RDX 和PETN 等算出其SHI，并與壓縮-剪切反應(yīng)動(dòng)力學(xué)的模擬結(jié)果比較表明：（1）空間位阻指數(shù)與反應(yīng)分子動(dòng)力學(xué)模擬所得結(jié)果如溫度、二氧化氮含量正相關(guān)。（2）對(duì)于正交晶系含能材料多個(gè)滑移系的SHI 的計(jì)算結(jié)果可知，SHI 大小與其撞擊感度有明顯的正相關(guān)。

文中提出的評(píng)價(jià)含能材料撞擊感度的算法適用于正交晶系含能材料快速查找出敏感滑移系與非敏感滑移系，以及評(píng)價(jià)正交晶系含能材料中撞擊感度的高低。算法中未考慮不同原子間碰撞化學(xué)因素的影響以及碰撞時(shí)的剪切應(yīng)力能壘，原則上將來(lái)可擴(kuò)展至非正交晶系含能材料。

基于空間位阻方法的感度評(píng)估不需要復(fù)雜耗時(shí)的分子動(dòng)力學(xué)模擬計(jì)算，因此能被用于高效評(píng)估含能材料感度，尤其適合于高通量計(jì)算的含能晶體結(jié)構(gòu)搜索和篩選，并可作為描述因子構(gòu)建關(guān)聯(lián)結(jié)構(gòu)-感度的機(jī)器學(xué)習(xí)模型。本研究對(duì)于篩選高能量密度低撞擊感度的含能材料具有一定的實(shí)用價(jià)值。