HMX/FOX-7共晶炸药分子动力学模拟

林鹤，张琳，朱顺官，李洪珍，彭新华

(1.南京理工大学 化工学院，江苏 南京210094;2.中国工程物理研究院 化工材料研究所，四川 绵阳621900)

0 引言

目前，随着高能炸药广泛应用，对炸药的要求越来越高。环四甲撑四硝胺(HMX)是一种爆速高、能量大的传统单质猛炸药。但由于HMX 的机械感度高，应用受到限制。应用溶剂侵蚀法制备表面光滑球形化HMX 晶粒［1］，超声波法原位合成钝感混合炸药HMX/TATB［2］和NTO(3-硝基-1，2，4-三唑-5-酮)在HMX 包覆［3］等方法均未能很大程度降低HMX 的机械感度。

超分子是分子通过分子间作用力形成的分子聚集体。共晶是2 种中性分子(2 种以上)通过主客体相互作用力组装而成的超分子复合物。主客体相互作用力主要有:氢键、范德华力，π-π 堆积和卤键［4-5］。在药物领域，药物活性组分(API)与共晶试剂(CCF)形成共晶，改善API 的物理化学性质，如:溶解度［6］、溶出速率和稳定性［7］。对于炸药领域，共晶在不破坏原有炸药分子化学结构情况下通过分子间作用力增大炸药晶体的密度［8］，提高炸药的爆速。此外，对于机械感度高的炸药，通过与钝感炸药形成共晶，从而降低机械感度。文献［9 -10］成功制备了摩尔比为1∶1的CL-20 和TNT 共晶炸药，还制备了17 种含有TNT 的共晶炸药，TNT 的晶体密度从1.70 g/cm3增加到1.73 g/cm3.文献［11 -12］利用分子动力学(MD)模拟HMX/TATB 共晶炸药并进行合成表征，HMX 的机械感度得到降低。同时共晶药物在起爆药领域也有进展，如文献［13 -14］制备研究了高氯酸铵盐(SY)起爆药，研究表明SY 的热安定性好、火焰感度低、燃速高和对静电不敏感。

1，1-二氨基-2，2-二硝基乙烯(FOX-7)是一种新型密度大(理论密度1.885 g/cm3［15］)钝感单质炸药。HMX 含有NO2基团，FOX-7 含有氨基，HMX 与FOX-7 可能以共晶形式存在。HMX 与FOX-7 形成共晶可以降低HMX 的机械感度。本文以共晶形成的机理搭建FOX-7 和HMX 共晶模型，并运用Materials Studio 软件［16］中Reflex［17］模块和Discover 模块，对FOX-7 和HMX 共晶进行了动力学模拟。通过考察共晶模型的结合能，径向分布函数，分子的键长、键角、二面角和XRD 衍射图判断共晶是否生成，为HMX/FOX-7 共晶的形成提供理论依据。

1 模拟方法

1.1 力场的选择

Compass［18］力场是在早期第2 代CFF9x 和PCFF 力场的基础上开发的一个ab initio 力场，能准确计算物质的热物理性质，广泛应用于有机分子、无机分子、聚合物的分子模拟。Compass 力场对HMX体系的适用性已得到证实［19］，HMX/FOX-7 共晶炸药采用Materials Studio 程序包中的Compass 力场进行分子力学(MM)和MD 的模拟。

1.2 共晶模型建立与MD 模拟

根据共晶形成的氢键规则［20］，FOX-7 和HMX可能通过FOX-7 的氨基与HMX 的硝基或FOX-7 的硝基与HMX 的亚甲基之间的分子间氢键组装成摩尔比为1∶1共晶。HMX/FOX-7 共晶炸药中FOX-7作为氢键的给体，HMX 作为氢键的受体。因此MD模拟过程中HMX 与FOX-7 的模型按照摩尔比为1∶1建立。基于β-HMX 的中子衍射单晶数据［21］，建立β-HMX 的晶胞和超晶胞(4 ×3 ×3)模型。FOX-7分子采取随机取代β-HMX 超晶胞和主要生长面2种模型。主要生长面模型采取切割分面法切割β-HMX 超晶胞，用36 个FOX-7 取代36 个HMX 分子，真空层厚度为30Å 构建。FOX-7 随机取代HMX超晶胞如图1所示。FOX-7 分子β-HMX 晶体的主要生长面用Morphology 模块中的grow morphology方法确定主要生长面。

初始模型通过MM 能量最小后，在Compass力场下，选取NVT 系综，采用Anderson［22］恒温器，温度设为298 K，用Velocity Verlet 积分法求解牛顿运动方程，原子初始速度按Boltzmann 随机分布方法确定，步长为1 fs，模拟时间为100 ps.模拟后得到最后平衡模型(FOX-7 随机取代HMX 超晶胞平衡模型见图1)。模拟过程中Van der Waals 和 coulomb 分 别 按Atom based［23］和Ewald［24］方法计算。

1.3 X 射线粉末衍射

将平衡后的模型运用Materials Studio 软件中Reflex 模块进行X-射线衍射(XRD)图模拟，选择Cu靶，2θ 为5°～75°范围内扫描，步长为0.02°.

图1 HMX/FOX-7 共晶模型Fig.1 Models of HMX/FOX-7 co-crystal in 1∶1 molar ratio

表1 HMX 的主要生长面Tab.1 Main stable surface of HMX

2 模拟与讨论

2.1 HMX/FOX-7 共晶体系的平衡

以HMX/FOX-7 共晶炸药(1，0，0)晶面为例，判别体系是否平衡必须同时满足温度平衡和能量平衡。从图2可知体系平衡后的温度上下波动不超过10%，从图3可知体系平衡后的能量波动较小，HMX/FOX-7 共晶模拟体系已达到平衡。其他模型具有相同结论。由图3可知，HMX 和FOX-7 存在着很强的分子间作用力，有利于HMX 与FOX-7 形成共晶。

图2 温度随时间变化Fig.2 Temperature vs simulation time

图3 能量随时间变化Fig.3 Energy vs simulation time

2.2 HMX 主要生长面的预测

HMX 经MM 能量最小化后，采用Morphology 模块中的growth morphology 预测模型，根据Compass力场和Compass 力场自带电荷，对HMX 晶体主要生长晶面进行预测。预测过程中，能量采用Forcite 模块计算，Van Der Waals 和Electrostatic 分别按照Atom based 和Ewald 方法计算。由表1知:HMX 有5个主要生长面，且其中(011)和(11)晶面占了90.99%.晶面附着能主要有范德华力和静电力构成，且总附着能Ea(0，1，1)＞Ea(0，2，0)＞Ea(1，1，)＞Ea(1，0)＞Ea(1，0，0).

2.3 HMX/FOX-7 共晶结合能

HMX/FOX-7 共晶可以通过HMX 的硝基和FOX-7 的氨基或HMX 的亚甲基和FOX-7 的硝基通过分子间形成氢键组装。二者之间的结合能可定义为ΔE=Ec- (EHMX+ EFOX-7)，Eb= - ΔE，式中:Ec为HMX/Fox-7 共晶平衡后的单点能(kcal/mol);EHMX为HMX/Fox-7 共晶平衡后去除FOX-7 之后的单点能(kcal/mol);EFOX-7为HMX/Fox-7 共晶平衡后去除HMX 之后的单点能(kcal/mol);结合能Eb定义为二者之间作用能ΔE 的负值(kcal/mol).

由表2可知:6 种共晶模型均有很强的结合能，结合能大小顺序为:Eb(020)＞Eb(100)＞Eb(附机)＞Eb(11)＞Eb(10)＞Eb(011);HMX/FOX-7在(020)晶面上结合能最大，体系能量最低，更加稳定，有利于HMX/FOX-7 共晶的形成。

2.4 HMX/FOX-7 共晶体系的径向分布函数

径向分布函数(RDF)为系统的区域密度与平均密度的比。参考分子的附近(r 值小)区域密度不同于系统的平均密度，但当参考分子距离远时区域密度应与平均密度相同，即当r 值大时RDF 接近1.g(r)通常可理解为给定某个粒子的坐标，其他粒子在空间的几何分布(离给定粒子多远)。通常情况下，分子间作用力包括氢键和范德华力。氢键长度在1.1～3.1Å 之间，强范德华力相互作用键长范围为3.1～5.0Å，弱范德华力相互作用键长大于5.0Å.

如图4(a)～图4(b)所示，HMX/FOX-7 6 种共晶模型在r =1.8Å 附近均有出峰，这表明HMX 的氧原子和FOX-7 的氢原子形成了分子间氢键，且(10)面共晶模型、(11)面共晶模型、(011)面共晶模型、(100)面共晶模型、随机取代共晶模型、(020)面共晶模型的峰值依次减少，说明形成氢键的作用力依次减少;在r=4.5Å 处，6 种共晶模型都有出峰，且峰值均大于形成氢键的峰值，说明6 种共晶模型均存在强范德华力作用。从图4(c)～图4(d)可知，HMX/FOX-7 6 种共晶模型在r =2.3Å 左右均有出峰，且随机取代模型在r =1.8Å 时有一个很强的峰出现，可见HMX 的亚甲基的氢与FOX-7的氧形成了氢键。氢键的强弱顺序为:随机取代共晶模型＞(10)共晶模型＞(11)共晶模型＞(011)共晶模型＞(100)共晶模型＞(020)共晶模型。同样存在强的范德华力作用。由此可以预测(10)面共晶模型存在着较强分子间作用力，较其他模型易形成HMX/FOX-7 共晶。

2.5 HMX/FOX-7 共晶的键长、键角和二面角的变化

MD 模拟过程中，考察HMX/FOX-7 共晶体系中键长，键角和二面角的变化。如图5所示，在MD 模拟过程中，FOX-7 中氨基中N…H 键，HMX 中硝基和HMX 构型都发生扭转和畸变。FOX-7 中氨基的键长由0.75Å 伸长到1.1Å.HMX 中硝基的角度从118Å 变化到138Å.特别是HMX 中部二面角波动，从-175°变化到175°，波动范围极大，进一步证实了HMX/FOX-7 共晶的形成。这些键长，键角，二面角发生扭转是由于HMX 和FOX-7 之间的氢键和范德华力的作用，使得原先的HMX 和FOX-结构造成破坏。

表2 HMX 各晶面与FOX-7 的结合能Tab.2 Binding energies between FOX-7 and the surfaces of HMX

2.6 XRD 图

经MD 平衡后，6 种共晶模型用Reflex 模块对其XRD 图进行了模拟。图1(b)中MD 平衡后的随机取代模型里有许多分子间和分子内氢键存在，是HMX/FOX-7 形成共晶的动力因素之一。

如图6所示，纯组分HMX 和FOX-7 的特征峰强在HMX/FOX-7 共晶模型衍射图中有所下降或消失，且随机共晶模型、(011)共晶模型、(11)共晶模型、(10)共晶模型在2θ 为6°附近均有新的衍射峰出现，(100)共晶模型、(020)共晶模型在2θ 为16°附近也有新的衍射峰出现。这些特征衍射峰消失或新的衍射峰出现不是HMX 和FOX-7 衍射峰简单的叠加。衍射峰的位置与物相有关，衍射峰的强度与物质的含量及晶体的大小等有关。任何一个物相都有一套晶面间距与衍射强度特真值，2 种不同物相的结构稍有差异其衍射图中的晶面间距与衍射强度。根据这些特征峰消失，衍射峰强度下降和新的衍射峰出现，证明了HMX/FOX-7 共晶模型是一种有别于纯组峰HMX 和FOX-7 的晶体，具有新的晶胞参数。主要是由于HMX 和FOX-7 形成共晶时，分子间氢键及范德华力对FOX-7 和HMX 的结构破坏，形成了一种新的结构。在上面分析中可知HMX/FOX-7 共晶MD 模拟过程中，FOX-7 中氨基键长，HMX 中硝基键角和HMX 中部二面角均有扭转、畸变。

图4 HMX/FOX-7 共晶结构的径向分布函数Fig.4 The radial distribution function g(r)of HMX/FOX-7 co-crystal

图5 FOX-7 中共晶键长、键角和二面角分布Fig.5 Distribution of specific bond length，bond angle and dihedral angle in HMX/FOX-7 co-crystal

3 结论

1)共晶模型中FOX-7 和HMX 存在较强的结合能，结合能大小顺序为:Eb(0 2 0)＞Eb(1 0 0)＞Eb(随机)＞Eb(1 1)＞Eb(1 0)＞Eb(0 1 1).

2)共晶模型的径向分布函数表明FOX-7 和HMX 分子间有氢键和范德华力作用，其中(1 0 2)共晶模型分子间作用力最大。

用MD 模拟6 种HMX/FOX-7 共晶模型，从理论上预测HMX/FOX-7 共晶形成。基于HMX 的晶体生长、结合能、径向分布函数的研究，FOX-7 取代HMX(111)晶面的共晶模型较其他5 种共晶模型形成。具体模拟结果如下:

3)MD 模拟过程中，FOX-7 氨基键长，HMX 硝基键角和HMX 中部二面角均发生了扭转和畸变，其中HMX 中部二面角的畸变最为严重。

图6 HMX/FOX-7 共晶模型XRD 粉末衍射Fig.6 X-ray power diffraction pattern of HMX/FOX-7 co-crystal models

4)6种共晶模型的XRD 衍射峰位置和强度均有别于纯组分HMX 或FOX-7，表明晶胞参数发生变化，一种新物相的形成。

References)

［1］ 徐容，李洪珍，黄明，等.球形化HMX 制备及性能研究［J］.含能材料，2010，18(5):505 -509.XU Rong，LI Hong-zhen，HUANG Ming，et al.Preparation and properties of rounded HMX［J］.Chinese Journal of Energetic Materials，2010，18(5):505 -509.(in Chinese)

［2］ 黄凤臣，邓明哲.原位合成钝感混合炸药HMX/TATB［J］.火炸药学报，2011，34(1):45 -48.HUANG Feng-chen，DENG Ming-zhe.Preparation of insensitive composite explosive HMX/TATB using ultrasonic method［J］.Chinese Journal of Explosives ＆ Propellants，2011，34(1):45 -48.(in Chinese)

［3］ Jung J W，Kim K J.Effect of supersaturation on the morphology of coated surface in coating by solution crystallization［J］.Industrial＆ Engineering Chemistry Research，2011，50(6):3475 -3482.

［4］ Clarke S M，Friscic T，Jones W，et al.Observation of a two-dimensional halogen-bonded cocrystal at sub-monolayer coverage using synchrotron X-ray diffraction［J］.Chemical Communications，2011，47(9):2526 -2528.

［5］ Alhalaweh A，George S，Bostrom D，et al.1∶1 and 2∶1 urea-succinic acid cocrystals:structural diversity，solution chemistry，and thermodynamic stability［J］.Crystal Growth ＆ Design，2010，11(10):4847 -4855.

［6］ Smith A J，Kavuru P，Wojtas L，et al.Cocrystal of quercetin with improved solubility and oral bioavailability［J］.Molecular Pharmaceutics，2011，8(5):1867 -1876.

［7］ Cherukuvada S，Babu N J，Nangia A.Nitrofurantoin-p-aminobenzoic acid cocrystal:hydration stability and dissolution rate studies［J］.Journal of Pharmaceutical Sciences，2011，100(8):3233 -3244.

［8］ Klapotke T M，Sabate C M.Bistetrazoles:nitrogen-rich，highperforming，insensitive energetic compounds［J］.Chemistry of Materials，2008，20(11):3629 -3637.

［9］ Atzger A J，Bolton O.Improved stability and smart-material functionality realized in an energetic cocrystal［J］.Angewandte Chemie International Edition，2011，50(38):8960 -8963.

［10］ Landenberger K B，Matzger A J.Cocrystal engineering of a prototype energetic materials supramolecular chemistry of 2，4，6-trinitrotoluene［J］.Crystal Growth ＆ Design，2010，10(12):5341-5347.

［11］ 卫春雪，段晓惠，李金山，等.环四甲撑四硝胺/1，3，5-三氨基-2，4，6-三硝基苯共晶炸药的分子模拟［J］.化学学报，2009，67(24):2822 -2826.WEI Chun-xue，DUAN Xiao-hui，LI Jin-shan，et al.Molecular simulation on co-crystal structure of HMX/TATB［J］.Acta Chimica Sinica，2009，67(24):2822 -2826.(in Chinese)

［12］ Shen J P，Duan X H，Luo Q P，et al.Preparation and characterization of a novel cocrystal explosive［J］.Crystal Growth ＆ Design，2011，(11):1759 -1765.

［13］ Ma P，Zhang L，Zhu S G，et al.Synthesis crystal structure and DFT calculation of an energetic perchlorate amine salt［J］.Journal of Crystal Growth，2011，335(1):70 -74.

［14］ 张垒，欧仙荣，朱顺官，等.胺盐共晶物为主装药的柔性导爆索研究［J］.爆破器材，2010，39(5):7 -9.ZHANG Lei，OU Xian-rong，ZHU Shun-guan，et al.Research on the initiating performance of MDF charged with amine perchlorate cocrystal compound［J］.Explosive Materials，2010，39(5):7 -9.(in Chinese)

［15］ Badgujar D M，Talawar M B，Asthana S N，et al.Advances in science and technology of modern energetic materials:An overview［J］.Journal of Hazardous Materials，2008，151(2 -3):289 -305.

［16］ Accelrys Inc.Materials Studio 3.0［CP］.San Diego，CA:Accelrys Inc.，2004.

［17］ Engel G E，Wilke S，Konig O，et al.Powder solve-a complete package for crystal structure solution from powder diffraction patterns［J］.Applied Crystallography，1999，32:1169 -1179.

［18］ Sun H.COMPASS:An ab initio force-field optimized for condensed-phase applications overview with details on alkane and benzene compounds［J］.Journal of Physical Chemistry B，1998，102:7338 -7364.

［19］ 朱伟，肖继军，赵峰，等.HMX/TATB 复合材料弹性性能的MD 模拟［J］.化学学报，2007，65(13):1223 -1228.ZHU Wei，XIAO Ji-yun，ZHAO Feng，et al.Molecular dynamics simulation of elastic properties of HMX/TATB composite［J］.Acta Chimica Sinica，2007，65(13):1223-1228.(in Chinese)

［20］ Panunto T W，Johnson R，Etter M C，et al.Hydrogen-bond formation in nitroanilines:the first step in designing acentric materials［J］.Journal of American Chemistry Society，1987，109:7786 -7797.

［21］ Choi C S，Boutin H P.A study of the crystal structure of β-Cyclotetramethylene tetranitramine by neutron diffraction［J］.Acta Crystallography，1970，B26:1240.

［22］ Andersen H C.Molecular dynamics simulations at constant pressure and/or temperature［J］.Journal of Chemical Physics，1980，72:2374 -2383.

［23］ Tosi M P.Cohesion of ionic solids in the Born model［J］.Solid State Physics，1964，16:1 -120.

［24］ Ewald P P.Evaluation of optical and electrostatic lattice potentials［J］.Annals of Physic，1921，64:253 -287.