LLM-105晶体形貌分子动力学模拟

李 蓉，甘 强，于 谦，程年寿，冯长根

(1.北京理工大学爆炸科学与技术国家重点实验室，北京100081；2.中国工程物理研究院化工材料研究所，四川 绵阳 621999)

引 言

1-氧-2,6-二氨基-3,5-二硝基吡嗪(LLM-105)是一种高能钝感含能材料，爆速高于TATB且易起爆[1]。在降感剂、传爆药、超高温石油射孔等领域具有广泛的应用前景[2]。LLM-105分子间与分子内存在大量的氢键[3]，DMSO是其重结晶过程最优的溶剂。但合成的LLM-105多存在“X”状孪晶或呈针状，安全性差、感度高，严重影响装药密度及成型性能。在实际结晶过程中，成核级数越大，越易爆发式成核，导致晶体缺陷。因此，控制LLM-105结晶过程、提高晶体品质对其应用具有重要意义。

目前LLM-105形貌调控的主要方法是溶剂/非溶剂法。李海波等[4]采用DMSO/水体系制备出长径比较大的六棱柱状和小团粒状晶体；Zhang等[5]用喷雾结晶法在DMSO溶剂中得到球状LLM-105晶体；庄小博等[6]用DMSO溶剂诱导实验发现，LLM-105在氢键作用下定向形成矩形棒状晶体，且晶体可能沿着(1 0 0)方向择优生长；徐容等[7]用DMSO和65%HNO3水溶液反相滴加重结晶得到颗粒状LLM-105，在DMSO溶剂中得到棒状晶体；李媛等[8]分别用DMSO、N-甲基吡咯烷酮和H2SO4为溶剂进行重结晶，发现LLM-105在DMSO中的溶解度较大，得到平均粒径5μm的棒状晶体。当前研究主要集中于LLM-105晶体粒径调控，而溶剂对LLM-105晶体形貌的影响规律报道较少。

近年来，量子力学和分子动力学方法已应用于含能材料模拟研究，平衡形态法、晶体生长法等已成功预测HMX、RDX、ANPyO等的晶体形貌[9]。1878年Gibbs[10]从热力学出发提出晶体生长最小表面能原理。平衡形态法基于该原理提出比表面能是晶体平行于该晶面脱附时所断裂的键能总和。晶体生长法[11]认为晶体的生长形态受到外部条件以及晶体内部结构特别是原子间键能的影响，是一非平衡过程。本研究采用分子动力学方法，通过平衡形态法和晶体生长法预测真空条件下LLM-105主要生长晶面。通过修正附着能分析DMSO环境下LLM-105主要控制晶面及其生长速率，通过径向分布函数和晶面结构分析不同晶面与DMSO之间的相互作用，考察溶剂DMSO在晶面上的扩散系数以分析LLM-105粒子成核能垒，以期为LLM-105晶体形貌控制提供理论参考。

1 分子动力学模拟

采用Materials Studio程序Visualizer模块，根据LLM-105单晶数据[12]构建晶胞模型。构建LLM-105超晶胞(4×3×4)并置于周期性边界条件的周期箱中，含192个LLM-105分子。采用COMPASS力场对LLM-105超晶胞进行能量最小化，优化晶胞参数。用Morphology模块计算LLM-105的晶体形态和晶体的主要生长晶面，对不同晶面分别构建3D周期结构，周期箱在c方向添加10nm的真空层。

根据文献[13]，采用分层建模方法构建DMSO/LLM-105“双层结构”初始模型。用Amorphous Cell模块建立溶剂层，选择900个溶剂分子，溶剂密度为1.507g/cm3。采用同样的方法进行结构优化，根据不同晶面的大小确定溶剂层面积。构建DMSO/LLM-105“双层结构”初始模型，其中上层为DMSO溶剂分子，下层分别为LLM-105(0 2 0)、(0 1 1)、(1 1 0)、(1 0 -1)晶面。双层结构模型经优化后，选取NPT系综进行分子动力学模拟，温度设为295K，模拟时间步长为1fs，总模拟步数100万步，共保存1000帧轨迹。模拟平衡以温度和能量为参照，通常温度和能量的值均须达到稳定体系才算平衡，一般在5%～10%以内的波动即认为体系已经达到平衡[14]。可以得到原子轨迹结构图和平衡结构图，对原子轨迹图前300ps用于热力学平衡，后700ps用于计算分析。温度控制采用Anderson方法，压力控制用Parrinello积分法，模拟过程中库仑力和范德华长程非键作用力分别采用Ewald和Atom based方法计算，截断半径为0.95nm。

2 结果与讨论

2.1 LLM-105晶体形貌预测

(1)平衡形态法

晶体在恒温和恒容的条件下，基于表面能分析晶体形貌即为平衡形态法，已成功运用于真空条件下含能材料晶体形貌的预测[15-16]。某一晶面的线性生长速率与其表面能成比例，定义为Gibbs-Wulff 晶体生长定律，如式(1)[17]所示。

(1)

式中：ri为平衡形态的晶体中心引向第i个晶面的距离；σi为第i个晶面的比表面自由能。

所预测LLM-105的晶体形貌如图1所示，主要生长晶面参数如表1所示。如果某一晶面的表面能小，则其形成临界晶核的能垒低，具有相对较高的成核速率。

从图1可见，LLM-105晶体类似于形状不规则的金字塔，长径比为1.766，体积为21317.920×10-3nm3，表面积为4216.635×10-2nm2，相对表面积为1.134，相似球形度为0.882。平衡形态法预测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 平衡形态法预测得到的LLM-105主要生长晶面的参数Table 1 Parameters of the main crystal faces of LLM-105 obtained by equilibrium morphology

注:n为晶面多重度

根据Bravais法则[18]，实际晶面往往平行于面网密度大的晶面，对生长质点吸引力小，生长速度慢，晶面的法线方向生长速率R与晶面间距dhkl成反比，即R∝1/dhkl。用平衡形态法计算晶面能的同时，还得到了真空条件下LLM-105不同晶面面积比例及晶面间距dhkl。从表1可见，LLM-105各晶面间距dhkl按照从大到小排列为(0 2 0)>(0 1 1)>(1 1 0)>(1 0 -1)，其中(0 2 0)生长较慢，与表面能分析结果一致。

(2)晶体生长法

在真空条件下晶面的线性生长速度R与其附着能Eatt成正比，即R∝Eatt。其中附着能是范德华力与静电相互作用能之和。晶体生长法预测得到LLM-105在真空中的晶型如图2所示。表2分别给出了(0 2 0)、(0 1 1)、(1 1 0)及(1 0 -1)4个主要生长晶面的参数。

从图2可见，与平衡形态法相比，晶体生长法预测的晶体形貌更趋于长柱状，长径比为2.094。从表2可见，晶面附着能绝对值按照由大到小排序为(1 1 0)>(1 0 -1)>(0 2 0)>(0 1 1)，其中(0 1 1)晶面占总晶面面积最大，与中心的距离最小，附着能最小，因此(0 1 1)晶面生长速率最小，是生长最慢的晶面。这一结果与平衡形态法预测结果有差异。两种方法分别从表面能和附着能进行研究，但是平衡形态法忽略了分子与基团间的键极性及生长条件的影响，仅能粗略预测晶体形貌，晶体生长法结果更为准确。考虑溶剂效应的影响，对已得到的附着能进行修正，得到修正附着能模型(AE模型)。溶剂层与晶面的相互作用能Eb如式(2)[16]所示

Eb=Etot-Esurf-Esolv

(2)

式中：Etot为整体结构能量；Esurf为晶面层的能量；Esolv为溶剂层的能量。

表2 晶体生长法预测得到LLM-105主要生长晶面的参数Table 2 Parameters of the main crystal faces of LLM-105 obtained by growth morphology method

注：S为单个晶面面积占总晶面表面积的百分比。

(3)

(4)

式中：Aacc为晶体单元胞晶面溶剂可达面积；Amodel为双层结构模型总晶面面积。

分子动力学模拟得到295K时DMSO/LLM-105双层模型能量最低的平衡结构，如图3所示。表3为LLM-105晶面的修正附着能。经过溶剂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分子内部有强相互作用，在真空条件下生长速率很快，极难控制晶体形貌，容易爆发式成核。加入溶剂DMSO后晶面附着能显著减小，有助于LLM-105晶体形貌控制。

2.2 晶面结构分析

采用径向分布函数(RDF)分析DMSO与LLM-105之间的相互作用，截断半径取1.6nm，每个计算都包括分子内和分子外作用。对4种晶面双层模型分别选取4对可能形成较强作用的原子对，分别是LLM-105中的N原子与DMSO中的H、O、S原子，以及LLM-105中的O原子与DMSO中的S原子。设置最大氢与受体距离为0.25nm，在供体、氢与受体之间键角最小为90°，RDF结果如图5所示。图6为LLM-105各个晶面与DMSO溶剂界面结构，红色虚线表示在0.0891nm原子间距离范围内的原子，蓝色虚线表示氢键相互作用原子。

分子间相互作用随距离变化通常分为氢键(0.26～0.31nm)和范德华力作用(0.31～0.50nm)，大于0.50nm范德华作用微弱。从图5可见，DMSO与LLM-105中4对原子径向分布函数中，LLM-105各个晶面上S—O和S—N原子对，在距离0.50nm之内没有峰值出现，说明没有强氢键或范德华力作用。对于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晶体形态学上是重要晶面，与晶体生长法预测结果一致。上述4个晶面均有大量活泼的O原子和H原子，因此可以选择合适的表面活性剂，通过氢键吸附降低晶体表面能，调节LLM-105的晶体形貌。

2.3 扩散系数分析

研究溶剂分子在含能材料中的扩散快慢有助于分析含能材料晶体的生长过程。根据Einstein方程[20]计算位移的时间平均获得扩散系数，计算式[21]如下：

(5)

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

(6)

式中：D为扩散系数；d为空间维度；t为扩散时间；r(0)和r(t)分别为起始坐标和t时刻的分子坐标。

在分子动力学模拟中三维立体晶胞模型d=3。在Discover模块中求解所有原子的运动学方程，再进行统计分布得到均方差位移MSD。溶剂DMSO分子在LLM-105中的扩散系数结果见表4。结合表3得到的修正附着能数据，得出了扩散系数与附着能的关系如图7所示。

表4 溶剂DMSO在LLM-105中的扩散系数Table 4 Diffusion coefficient of solvent DMSO in LLM-105

注：(MSD/t)为斜率

由表4可见，(0 1 1)晶面扩散系数最大，说明该晶面上溶剂最容易扩散，成核能垒高不易聚集成核。结合表3中(0 1 1)晶面修正附着能最小，说明与DMSO溶剂分子相容性较好，在LLM-105间隙中发生空间位移。从图7可见，修正附着能随着溶剂扩散系数的增大而减小，推测溶剂DMSO在不同晶面上扩散系数不同导致不同晶面修正附着能不同。(1 0 -1)晶面修正附着能大、成核能垒低，使LLM-105晶体容易团聚成核，导致分子扩散系数小。

LLM-105在DMSO中溶解度较大，表面能快速减小，当表面能减小到一定程度时，LLM-105颗粒会快速聚集在一起，溶剂诱导作用增强[22]。且DMSO属于极性溶剂，更加促进了诱导扩散，这一过程与含能材料在增塑剂中的扩散过程相似[23]。

3 结 论

(1)考察了真空和溶剂DMSO条件下LLM-105主要生长晶面情况，结果表明，真空条件下采用平衡形态法和晶体生长法，分别预测得到(0 2 0)和(0 1 1)晶面为LLM-105晶体生长控制晶面。分析溶剂DMSO作用后的修正附着能表明，(0 1 1)晶面为生长控制晶面，生长速率最小。

(2)溶剂中4个晶面附着能远小于真空条件下，说明选择合适的极性溶剂有利于LLM-105晶体形貌调控。

(3)径向分布函数和晶面结构分析表明，DMSO中的H原子与(0 1 1)晶面LLM-105中的N原子对之间具有较强氢键和范德华力相互作用，极性晶面(0 1 1)显露完整的硝基，有利于形成分子间相互作用。

(4)DMSO在不同晶面上扩散系数考察表明，修正附着能随着溶剂扩散系数的增大而减小，溶剂DMSO在(0 1 1)晶面上扩散系数最大，因此成核能垒高，不容易聚集成核。

参考文献:

[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] 周杰文, 张创军, 王友兵,等. 耐热炸药的现状及研究进展[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重结晶与性能研究[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] 徐容, 廖龙渝, 王述存, 等. 重结晶方法对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的制备及粒度控制技术研究[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] 肖继军,朱卫华,朱伟, 等.高能材料分子动力学[M].北京:科学出版社, 2013.

[15] 杨利, 任晓婷, 严英俊, 等. 六硝基茋的晶体结构及形貌模拟[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] 何前军, 黄志良, 张联盟. 氢键影响羟基磷灰石结晶形态的理论模型研究[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] 段晓惠, 卫春雪, 裴重华,等. 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缺陷晶体及其复合体系结构与性能的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.