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

（1. 上海大學理學院化學系，上海 200444；2. 上海大學材料基因組工程研究院，上海 200444；3. 中國工程物理研究院化工材料研究所，四川 綿陽 621999）

1 引言

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

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

2 設計目標與算法

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

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

圖1 分子在剪切形變過程中的空間位阻示意圖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 壓縮

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

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

2.3 旋轉(zhuǎn)

由于每一個滑移系的滑移面和滑移方向上的空間位阻存在差別，為了保證算法的一致性需要建立新的空間直角坐標系，使得滑移方向為新坐標系的X 軸正方向，滑移面為新坐標系的XOZ 平面，如圖3 所示。

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

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

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

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

2.4 分層

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

2.5 投影

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

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

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

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

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

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

2.6 滑移系的自動搜尋

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

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

3.1 投影算法的評判

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

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

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

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

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

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

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

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

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

為 與An 等［5］報 道 的RDX 的CS-RD 對 比，同 樣 選定壓縮比（r）為0.1 和0.2。壓縮比增加則系內(nèi)晶體分子碰撞后溫度與NO2含量都會相應增加。對滑移系的SHI 計算結(jié)果分別見表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］報 道 的PETN 的CS-RD 對 比，選定壓縮比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］報 道 的β-HMX 的CS-RD 對 比，選定壓縮比ratio=0.17。

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

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

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

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

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

3.6 各正交晶系含能材料的SHI 與落錘實驗測得撞擊感度比較

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

表3 r=0.2 時24 個滑移系的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 時8 個滑移系的SHI 與CS-RD 結(jié)果［3］比較Table 4 Comparison of SHI and CS-RD results for 8 slip systems with r=0.1

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

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

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

表7 RDX、PETN、BTF、TNT 在22 個低指數(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 對含能材料10%的壓縮剪切后22 個滑移系的SHI 對照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

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

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

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

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

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

4 結(jié)論

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

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

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