LLM-105晶體形貌分子動力學模擬

李 蓉，甘 強，于 謙，程年壽，馮長根

(1.北京理工大學爆炸科學與技術(shù)國家重點實驗室，北京100081；2.中國工程物理研究院化工材料研究所，四川 綿陽 621999)

引 言

1-氧-2,6-二氨基-3,5-二硝基吡嗪(LLM-105)是一種高能鈍感含能材料，爆速高于TATB且易起爆[1]。在降感劑、傳爆藥、超高溫石油射孔等領(lǐng)域具有廣泛的應用前景[2]。LLM-105分子間與分子內(nèi)存在大量的氫鍵[3]，DMSO是其重結(jié)晶過程最優(yōu)的溶劑。但合成的LLM-105多存在“X”狀孿晶或呈針狀，安全性差、感度高，嚴重影響裝藥密度及成型性能。在實際結(jié)晶過程中，成核級數(shù)越大，越易爆發(fā)式成核，導致晶體缺陷。因此，控制LLM-105結(jié)晶過程、提高晶體品質(zhì)對其應用具有重要意義。

目前LLM-105形貌調(diào)控的主要方法是溶劑/非溶劑法。李海波等[4]采用DMSO/水體系制備出長徑比較大的六棱柱狀和小團粒狀晶體；Zhang等[5]用噴霧結(jié)晶法在DMSO溶劑中得到球狀LLM-105晶體；莊小博等[6]用DMSO溶劑誘導實驗發(fā)現(xiàn)，LLM-105在氫鍵作用下定向形成矩形棒狀晶體，且晶體可能沿著(1 0 0)方向擇優(yōu)生長；徐容等[7]用DMSO和65%HNO3水溶液反相滴加重結(jié)晶得到顆粒狀LLM-105，在DMSO溶劑中得到棒狀晶體；李媛等[8]分別用DMSO、N-甲基吡咯烷酮和H2SO4為溶劑進行重結(jié)晶，發(fā)現(xiàn)LLM-105在DMSO中的溶解度較大，得到平均粒徑5μm的棒狀晶體。當前研究主要集中于LLM-105晶體粒徑調(diào)控，而溶劑對LLM-105晶體形貌的影響規(guī)律報道較少。

近年來，量子力學和分子動力學方法已應用于含能材料模擬研究，平衡形態(tài)法、晶體生長法等已成功預測HMX、RDX、ANPyO等的晶體形貌[9]。1878年Gibbs[10]從熱力學出發(fā)提出晶體生長最小表面能原理。平衡形態(tài)法基于該原理提出比表面能是晶體平行于該晶面脫附時所斷裂的鍵能總和。晶體生長法[11]認為晶體的生長形態(tài)受到外部條件以及晶體內(nèi)部結(jié)構(gòu)特別是原子間鍵能的影響，是一非平衡過程。本研究采用分子動力學方法，通過平衡形態(tài)法和晶體生長法預測真空條件下LLM-105主要生長晶面。通過修正附著能分析DMSO環(huán)境下LLM-105主要控制晶面及其生長速率，通過徑向分布函數(shù)和晶面結(jié)構(gòu)分析不同晶面與DMSO之間的相互作用，考察溶劑DMSO在晶面上的擴散系數(shù)以分析LLM-105粒子成核能壘，以期為LLM-105晶體形貌控制提供理論參考。

1 分子動力學模擬

采用Materials Studio程序Visualizer模塊，根據(jù)LLM-105單晶數(shù)據(jù)[12]構(gòu)建晶胞模型。構(gòu)建LLM-105超晶胞(4×3×4)并置于周期性邊界條件的周期箱中，含192個LLM-105分子。采用COMPASS力場對LLM-105超晶胞進行能量最小化，優(yōu)化晶胞參數(shù)。用Morphology模塊計算LLM-105的晶體形態(tài)和晶體的主要生長晶面，對不同晶面分別構(gòu)建3D周期結(jié)構(gòu)，周期箱在c方向添加10nm的真空層。

根據(jù)文獻[13]，采用分層建模方法構(gòu)建DMSO/LLM-105“雙層結(jié)構(gòu)”初始模型。用Amorphous Cell模塊建立溶劑層，選擇900個溶劑分子，溶劑密度為1.507g/cm3。采用同樣的方法進行結(jié)構(gòu)優(yōu)化，根據(jù)不同晶面的大小確定溶劑層面積。構(gòu)建DMSO/LLM-105“雙層結(jié)構(gòu)”初始模型，其中上層為DMSO溶劑分子，下層分別為LLM-105(0 2 0)、(0 1 1)、(1 1 0)、(1 0 -1)晶面。雙層結(jié)構(gòu)模型經(jīng)優(yōu)化后，選取NPT系綜進行分子動力學模擬，溫度設(shè)為295K，模擬時間步長為1fs，總模擬步數(shù)100萬步，共保存1000幀軌跡。模擬平衡以溫度和能量為參照，通常溫度和能量的值均須達到穩(wěn)定體系才算平衡，一般在5%～10%以內(nèi)的波動即認為體系已經(jīng)達到平衡[14]?？梢缘玫皆榆壽E結(jié)構(gòu)圖和平衡結(jié)構(gòu)圖，對原子軌跡圖前300ps用于熱力學平衡，后700ps用于計算分析。溫度控制采用Anderson方法，壓力控制用Parrinello積分法，模擬過程中庫侖力和范德華長程非鍵作用力分別采用Ewald和Atom based方法計算，截斷半徑為0.95nm。

2 結(jié)果與討論

2.1 LLM-105晶體形貌預測

(1)平衡形態(tài)法

晶體在恒溫和恒容的條件下，基于表面能分析晶體形貌即為平衡形態(tài)法，已成功運用于真空條件下含能材料晶體形貌的預測[15-16]。某一晶面的線性生長速率與其表面能成比例，定義為Gibbs-Wulff 晶體生長定律，如式(1)[17]所示。

(1)

式中：ri為平衡形態(tài)的晶體中心引向第i個晶面的距離；σi為第i個晶面的比表面自由能。

所預測LLM-105的晶體形貌如圖1所示，主要生長晶面參數(shù)如表1所示。如果某一晶面的表面能小，則其形成臨界晶核的能壘低，具有相對較高的成核速率。

從圖1可見，LLM-105晶體類似于形狀不規(guī)則的金字塔，長徑比為1.766，體積為21317.920×10-3nm3，表面積為4216.635×10-2nm2，相對表面積為1.134，相似球形度為0.882。平衡形態(tài)法預測LLM-105晶體主要生長晶面為(0 2 0)、(0 1 1)、(1 1 0)及(1 0 -1)晶面。從表1可見，各晶面表面能按照由大到小的順序為(1 1 0)>(1 0 -1)>(0 1 1)>(0 2 0)，其中(0 2 0)晶面表面能最低，生長最慢。

表1 平衡形態(tài)法預測得到的LLM-105主要生長晶面的參數(shù)Table 1 Parameters of the main crystal faces of LLM-105 obtained by equilibrium morphology

注:n為晶面多重度

根據(jù)Bravais法則[18]，實際晶面往往平行于面網(wǎng)密度大的晶面，對生長質(zhì)點吸引力小，生長速度慢，晶面的法線方向生長速率R與晶面間距dhkl成反比，即R∝1/dhkl。用平衡形態(tài)法計算晶面能的同時，還得到了真空條件下LLM-105不同晶面面積比例及晶面間距dhkl。從表1可見，LLM-105各晶面間距dhkl按照從大到小排列為(0 2 0)>(0 1 1)>(1 1 0)>(1 0 -1)，其中(0 2 0)生長較慢，與表面能分析結(jié)果一致。

(2)晶體生長法

在真空條件下晶面的線性生長速度R與其附著能Eatt成正比，即R∝Eatt。其中附著能是范德華力與靜電相互作用能之和。晶體生長法預測得到LLM-105在真空中的晶型如圖2所示。表2分別給出了(0 2 0)、(0 1 1)、(1 1 0)及(1 0 -1)4個主要生長晶面的參數(shù)。

從圖2可見，與平衡形態(tài)法相比，晶體生長法預測的晶體形貌更趨于長柱狀，長徑比為2.094。從表2可見，晶面附著能絕對值按照由大到小排序為(1 1 0)>(1 0 -1)>(0 2 0)>(0 1 1)，其中(0 1 1)晶面占總晶面面積最大，與中心的距離最小，附著能最小，因此(0 1 1)晶面生長速率最小，是生長最慢的晶面。這一結(jié)果與平衡形態(tài)法預測結(jié)果有差異。兩種方法分別從表面能和附著能進行研究，但是平衡形態(tài)法忽略了分子與基團間的鍵極性及生長條件的影響，僅能粗略預測晶體形貌，晶體生長法結(jié)果更為準確?？紤]溶劑效應的影響，對已得到的附著能進行修正，得到修正附著能模型(AE模型)。溶劑層與晶面的相互作用能Eb如式(2)[16]所示

Eb=Etot-Esurf-Esolv

(2)

式中：Etot為整體結(jié)構(gòu)能量；Esurf為晶面層的能量；Esolv為溶劑層的能量。

表2 晶體生長法預測得到LLM-105主要生長晶面的參數(shù)Table 2 Parameters of the main crystal faces of LLM-105 obtained by growth morphology method

注：S為單個晶面面積占總晶面表面積的百分比。

(3)

(4)

式中：Aacc為晶體單元胞晶面溶劑可達面積；Amodel為雙層結(jié)構(gòu)模型總晶面面積。

分子動力學模擬得到295K時DMSO/LLM-105雙層模型能量最低的平衡結(jié)構(gòu)，如圖3所示。表3為LLM-105晶面的修正附著能。經(jīng)過溶劑DMSO作用后預測的晶體形貌如圖4所示。

表3 溶劑中DMSO/LLM-105晶面的修正附著能Table 3 Modified attachment energies of LLM-105 crystal faces in solvent

注：R′為溶劑作用之后的生長速率。

對比真空條件與溶劑條件下LLM-105附著能(表2和表3)，可見溶劑中4個晶面的附著能遠遠小于真空條件下，這是由于LLM-105分子內(nèi)部有強相互作用，在真空條件下生長速率很快，極難控制晶體形貌，容易爆發(fā)式成核。加入溶劑DMSO后晶面附著能顯著減小，有助于LLM-105晶體形貌控制。

2.2 晶面結(jié)構(gòu)分析

采用徑向分布函數(shù)(RDF)分析DMSO與LLM-105之間的相互作用，截斷半徑取1.6nm，每個計算都包括分子內(nèi)和分子外作用。對4種晶面雙層模型分別選取4對可能形成較強作用的原子對，分別是LLM-105中的N原子與DMSO中的H、O、S原子，以及LLM-105中的O原子與DMSO中的S原子。設(shè)置最大氫與受體距離為0.25nm，在供體、氫與受體之間鍵角最小為90°，RDF結(jié)果如圖5所示。圖6為LLM-105各個晶面與DMSO溶劑界面結(jié)構(gòu)，紅色虛線表示在0.0891nm原子間距離范圍內(nèi)的原子，藍色虛線表示氫鍵相互作用原子。

分子間相互作用隨距離變化通常分為氫鍵(0.26～0.31nm)和范德華力作用(0.31～0.50nm)，大于0.50nm范德華作用微弱。從圖5可見，DMSO與LLM-105中4對原子徑向分布函數(shù)中，LLM-105各個晶面上S—O和S—N原子對，在距離0.50nm之內(nèi)沒有峰值出現(xiàn)，說明沒有強氫鍵或范德華力作用。對于H—N原子對僅觀察到LLM-105(0 1 1)晶面中峰值明顯，具有強氫鍵作用和范德華力相互作用。(1 0 -1)晶面上O—N原子對在距離0.3nm左右有一個峰值，但概率密度并不大，表明有較弱的氫鍵作用。

晶面上垂直顯露的硝基有利于與周圍溶劑分子形成氫鍵[19]。由圖6可見，(0 2 0)晶面顯露的是氫原子，其他晶面都有強極性的硝基基團。其中(0 1 1)晶面顯露完整的硝基，(1 1 0)晶面主要顯露的是N—O基團，(1 0 -1)晶面顯露原子排列整齊但不完整。按晶面極性來劃分，(0 1 1)屬于極性晶面，(1 1 0)與(1 0 -1)極性較弱，(0 2 0)晶面屬于非極性晶面。上述分析可見，(0 1 1)晶面上的硝基與DMSO形成氫鍵，使溶劑脫附變得困難，阻礙了LLM-105分子在晶面上沉積，晶體生長速率變緩，因而在LLM-105晶體形態(tài)學上是重要晶面，與晶體生長法預測結(jié)果一致。上述4個晶面均有大量活潑的O原子和H原子，因此可以選擇合適的表面活性劑，通過氫鍵吸附降低晶體表面能，調(diào)節(jié)LLM-105的晶體形貌。

2.3 擴散系數(shù)分析

研究溶劑分子在含能材料中的擴散快慢有助于分析含能材料晶體的生長過程。根據(jù)Einstein方程[20]計算位移的時間平均獲得擴散系數(shù)，計算式[21]如下：

(5)

MSD=〈r(t)-r(0)〉2

(6)

式中：D為擴散系數(shù)；d為空間維度；t為擴散時間；r(0)和r(t)分別為起始坐標和t時刻的分子坐標。

在分子動力學模擬中三維立體晶胞模型d=3。在Discover模塊中求解所有原子的運動學方程，再進行統(tǒng)計分布得到均方差位移MSD。溶劑DMSO分子在LLM-105中的擴散系數(shù)結(jié)果見表4。結(jié)合表3得到的修正附著能數(shù)據(jù)，得出了擴散系數(shù)與附著能的關(guān)系如圖7所示。

表4 溶劑DMSO在LLM-105中的擴散系數(shù)Table 4 Diffusion coefficient of solvent DMSO in LLM-105

注：(MSD/t)為斜率

由表4可見，(0 1 1)晶面擴散系數(shù)最大，說明該晶面上溶劑最容易擴散，成核能壘高不易聚集成核。結(jié)合表3中(0 1 1)晶面修正附著能最小，說明與DMSO溶劑分子相容性較好，在LLM-105間隙中發(fā)生空間位移。從圖7可見，修正附著能隨著溶劑擴散系數(shù)的增大而減小，推測溶劑DMSO在不同晶面上擴散系數(shù)不同導致不同晶面修正附著能不同。(1 0 -1)晶面修正附著能大、成核能壘低，使LLM-105晶體容易團聚成核，導致分子擴散系數(shù)小。

LLM-105在DMSO中溶解度較大，表面能快速減小，當表面能減小到一定程度時，LLM-105顆粒會快速聚集在一起，溶劑誘導作用增強[22]。且DMSO屬于極性溶劑，更加促進了誘導擴散，這一過程與含能材料在增塑劑中的擴散過程相似[23]。

3 結(jié) 論

(1)考察了真空和溶劑DMSO條件下LLM-105主要生長晶面情況，結(jié)果表明，真空條件下采用平衡形態(tài)法和晶體生長法，分別預測得到(0 2 0)和(0 1 1)晶面為LLM-105晶體生長控制晶面。分析溶劑DMSO作用后的修正附著能表明，(0 1 1)晶面為生長控制晶面，生長速率最小。

(2)溶劑中4個晶面附著能遠小于真空條件下，說明選擇合適的極性溶劑有利于LLM-105晶體形貌調(diào)控。

(3)徑向分布函數(shù)和晶面結(jié)構(gòu)分析表明，DMSO中的H原子與(0 1 1)晶面LLM-105中的N原子對之間具有較強氫鍵和范德華力相互作用，極性晶面(0 1 1)顯露完整的硝基，有利于形成分子間相互作用。

(4)DMSO在不同晶面上擴散系數(shù)考察表明，修正附著能隨著溶劑擴散系數(shù)的增大而減小，溶劑DMSO在(0 1 1)晶面上擴散系數(shù)最大，因此成核能壘高，不容易聚集成核。

參考文獻:

[1] Pagoria P F. Synthesis, scale-up, and characterization of 2,6-diamino-3,5-dinitropyrazine-1-oxide (LLM-105) [J]. Propellants, Explosives, Pyrotechnics, 1998, 23:156-160.

[2] 周杰文, 張創(chuàng)軍, 王友兵,等. 耐熱炸藥的現(xiàn)狀及研究進展[J]. 兵器裝備工程學報, 2016, 37(3):111-115.

ZHOU Jie-wen, ZHANG Chuang-jun, WANG You-bing, et al. Status and research progress of heat resistant explosives[J]. Jounal of Ordance Equipment Engineering, 2016, 37(3):111-115.

[3] He W D, Zhou G, Wong N B, et al. Intramolecular H-bonds in LLM-105 and its derivatives: a DFT study[J]. Journal of Molecular Structure: Theochem, 2005, 723(1/3):217-222.

[4] 李海波, 程碧波, 劉世俊, 等. LLM-105重結(jié)晶與性能研究[J]. 含能材料, 2008, 16(6):686-688.

LI Hai-bo，CHENG Bi-bo, LIU Shi-jun, et al. Recrystallization and properties of LLM-105 [J]. Chinese Journal of Energetic Materials, 2008, 16(6):686-688.

[5] Zhang J, Wu P, Yang Z, et al. Preparation and properties of submicrometer-sized LLM-105 via spray-crystallization method[J].Chinese Journal of Energetic Materials, 2015, 39(5): 653-657.

[6] 莊小博, 黃兵, 高冰, 等. 納米LLM-105自組裝制備矩形微米棒[J]. 含能材料, 2016, 24(5):433-438.

ZHUANG Xiao-bo, HUANG Bing, GAO Bing, et al. Preparation of rectangular micro-rods nano-LLM-105 self-assembly[J]. Chinese Journal of Energetic Materials, 2016, 24(5): 433-438.

[7] 徐容, 廖龍渝, 王述存, 等. 重結(jié)晶方法對2,6-二氨基-3,5-二硝基吡嗪-1-氧化物晶體特性及性能影響[J]. 兵工學報, 2015, 36(11):2099-2103.

XU Rong, LIAO Long-yu, WANG Shu-cun, et al. Influence of recryztallization method on crystal properties and performance of LLM-105[J]. Acta Amamentarii, 2015, 36(11):2099-2013.

[8] 李媛, 黃鳳臣, 王友兵, 等. 細顆粒LLM-105的制備及粒度控制技術(shù)研究[J]. 火工品, 2017(1):34-37.

LI Yuan, HUANG Feng-chen, WANG You-bing, et al. Preparation and research on particle size control of LLM-105[J]. Initiators & Pyrotechnics, 2017(1): 34-37.

[9] 湯嶄, 楊利, 喬小晶, 等. HMX晶體形貌的計算模擬[J]. 火炸藥學報, 2009, 32(4):10-13.

TAN Zhan, YANG Li, QIAO Xiao-jing, et al. Calculated simulation of the crystal morphology of HMX[J]. Chinese Journal of Explosives & Propellants(Huozhayao Xuebao), 2009, 32(4):10-13.

[10] Gibbs J W. On the equilibrium of heterogeneous substances [J]. Transactions of the Connecticut Academy of Arts and Sciences, 1878, 3: 441-458.

[11] Hartman P, Chan H K. Application of the periodic bond chain (PBC) theory and attachment energy consideration to derive the crystal morphology of hexamethylmelamine[J]. Pharmaceutical Research, 1993, 10(7):1052-8.

[12] Liu Jinjian, Liu Zuliang, Cheng Jian, et al. Synthesis, crystal structure and catalytic properties of 2,6 -diamino-3,5-dinitropyridine -1-oxide coblat(Ⅲ)[J]. Chinese Journal of Inorganic Chemistry, 2013, 29(2):289-294.

[13] 石文艷, 王風云, 夏明珠,等. 2,6-二氨基-3,5-二硝基吡啶-1-氧化物晶體形貌的MD模擬[J]. 含能材料, 2016, 24(1):19-26.

SHI Wen-yan,WANG Feng-yun, XIA Ming-zhu,et al.Molecular dynamics simulation on the crystal morpholoqy of 2,6-diamino-3,5-dinitropyridine-1-oxide[J].Chinese Journal of Energetic Materials,2016,24(1):19-26.

[14] 肖繼軍,朱衛(wèi)華,朱偉, 等.高能材料分子動力學[M].北京:科學出版社, 2013.

[15] 楊利, 任曉婷, 嚴英俊, 等. 六硝基茋的晶體結(jié)構(gòu)及形貌模擬[J]. 火炸藥學報, 2009, 32(6):1-5.

YANG Li, REN Xiao-ting, YAN Ying-jun, et al. Crystal structure and morphology simulation of HNS[J]. Chinese Journal of Explosives and Propellants(Huozhayao Xuebao), 2009, 32(6):1-5.

[16] 任曉婷, 楊利, 張國英, 等. TATB晶體形貌的計算模擬[J]. 火炸藥學報, 2010, 33(6):43-46.

REN Xiao-ting, YANG Li, ZHANG Guo-ying, et al. Computional simulation of the crystal morphology of TATB[J]. Chinese Journal of Explosives and Propellants(Huozhayao Xuebao), 2010, 33(6):43-46.

[17] 何前軍, 黃志良, 張聯(lián)盟. 氫鍵影響羥基磷灰石結(jié)晶形態(tài)的理論模型研究[J]. 人工晶體學報, 2006, 35(1):147-151.

HE Qian-jun, HUANG Zhi-liang, ZHANG Lian-meng. Study on the model of crystal morphology of synthesized hydroxyapatite affected by hydrogen bond[J]. Journal of Synthetic Crystals, 2006, 35(1):147-151.

[18] Donnay J D H, Harker D. A new law of crystal morphology extending the law of bravais[J]. American Mineralogist, 1938, 22: 446-467.

[19] 段曉惠, 衛(wèi)春雪, 裴重華,等. HMX晶體形貌預測[J]. 含能材料, 2009, 17(6):655-659.

DUAN Xiao-hui, WEI Chun-xue，PEI Chong-hua, et al. Prediction of crystal morphology of HMX[J]. Chinese Journal of Energetic Materials, 2009, 17(6):655-659.

[20] Haesslin H W. Dimethylsiloxane-ethylene oxide block copolymers, 2. Preliminary results on dilute solution properties[J]. Macromolecular Chemistry & Physics, 1985, 186(2):357-366.

[21] 曹強. CL-20/TNT共晶、RDX缺陷晶體及其復合體系結(jié)構(gòu)與性能的MD模擬研究[D]. 南京: 南京理工大學, 2015.

CAO Qiang. Molecular dynamics simulation study on the structure and performance of CL-20/TNT cocrystal/RDX defective crystal and their composite materials[D]. Nanjing: Nanjing University of Science and Technology, 2015.

[22] 莊小博. 典型不敏感炸藥納米顆粒的自組裝研究[D]. 綿陽: 西南科技大學, 2016.

ZHUANG Xiao-bo. Self-assembly of a typical insensitive explosive nanoparticles[D]. Mianyang: Southwest University of Science and Technology, 2016.

[23] 李紅霞, 強洪夫, 王廣, 等. 基于MD方法的增塑劑擴散行為的模擬研究[J]. 含能材料, 2009, 17(1): 36-41.

LI Hong-xia, QIANG Hong-fu, WANG Guang, et al. Molecular dynamics simulation of plasticizer diffusion[J]. Chinese Journal of Energetic Materials, 2009, 17(1):36-41.