1，1-二氨基-2，2-二硝基乙烯晶体热膨胀各向异性的模拟研究

钱 文，张朝阳，宗和厚，熊 鹰，张伟斌，舒远杰

（中国工程物理研究院化工材料研究所 四川绵阳621900）

1 引 言

1，1-二氨基-2，2-二硝基乙烯，俗称FOX-7，是1998年由瑞典国防研究院（Swedish Defence Research Agency，FOI）的一个研究小组首先报道合成的[1]，近几年来成为高能量密度化合物的研究热点之一．FOX-7分子（如图1）属于推－拉型烯烃，分子内存在大的共轭体系，硝基氧原子与氨基氢原子之间还存在氢键，含两个硝基因此高能而又有两个胺基降感，是钝感高能炸药（IHE）的理想组分．

图1 FOX-7的分子构型Fig．1 Molecular structure of FOX-7

通常的结晶条件下均得到α晶型FOX-7晶体，属于单斜晶系 P21／n空间群[2，15]．FOX-7晶体广泛应用于低感度炸药配方，炸药晶体热膨胀性质具有重要意义：材料内部晶体热膨胀的各向异性，将引发显著的内应力[3]，如果内应力的作用超过了内聚能，将最终导致内部裂纹的出现；而材料内部晶体受热的各向异性膨胀，引起材料热形变[4]，产生不可逆的尺寸或密度变化，将影响精密测量的准确度．为了确保含FOX-7炸药部件的精确测量和库存可靠性，需要对FOX-7晶体特性进行深入研究．

针对FOX-7晶体的性质各国科学家在实验方面开展了一系列工作．Bemm 等[2]和 Gilardi等[15]先后用X射线衍射研究了FOX-7的晶体结构，发现了其层状结构；Suhithi等[5]用X射线衍射和拉曼光谱研究了在高压下FOX-7晶体的状态与结构变化，发现压力在1．1GPa以下时晶胞b轴的压缩程度远大于a轴和c轴的压缩程度；Kempa和Herrmann等[6]用 DSC、TG、TMA 等热分析方法和X射线衍射研究了FOX-7晶体的相变行为，观察到晶体的α、β、γ三相；Burnham 和 Weese等[7，8]用 DSC、TG、HPC、ARC 及 TG-DTA-FTIR-MS联用技术详细研究了FOX-7的固体相变、热分解和活化能等热性质；Evers等[9]采用X射线单晶和粉末衍射详细研究了α、β两种晶型的FOX-7在200K-393K温度范围的晶体特性、热膨胀和相变．理论研究主要集中在对FOX-7晶体的各种性能计算上，第一性原理计算和分子模拟的结果不仅与实验相互印证，更成为深入研究本质机理的有效途径．赵纪军等[10]用密度泛函理论方法计算了高压下FOX-7的晶胞参数和能带结构，所得结果与实验值符合；Hu，Anguang等[11]采用DFT方法分析了FOX-7晶体热分解反应，发现了在高压下分子间氢转移导致的化学分解；Dan C．Sorescu等[12]采用密度泛函理论方法计算了FOX-7的结构性质，结果与X射线衍射实验数据相符，并通过自建分子间作用势和分子动力学模拟研究了其各向异性性质；DeCarlos E．Taylor等[13]以SAPT理论（symmetry adapted perturbation theory）为基础建立分子力场（势函数），对FOX-7晶体进行等温等压下的动力学模拟，结果与实验结果相符并且其随温度和压力变化体现显著的各向异性；边亮等[14]采用巨正则系综蒙特卡洛、密度泛函理论方法结合分子动力学模拟，详细研究了FOX-7晶体的密度、表面能、形貌和爆轰性能．以上研究为深入阐释FOX-7晶体的热膨胀机理及进行相关预测提供了新的思路，本文中运用分子动力学模拟（MD）结合量子化学计算（QM）方法研究了FOX-7晶体的热膨胀的各向异性，将能量变化及分子堆积模式与热膨胀各向异性建立了关联．

2 模拟计算方法与细节

2．1 分子模型的构建与分析

FOX-7晶体结构基于α-FOX-7单晶X射线衍射数据[15]构建，晶胞参数为a＝0．6940nm，b＝0．6637nm，c＝1．1341nm，β＝90．61°，α＝γ＝90°（295K），呈波浪型π堆积．为在有限的计算资源下，尽可能获得能够代表实际体系的分子模型，需要以单胞结构为基础构建超晶胞作为分子模拟的初始构型．分别构建了3＊3＊2超晶胞（图3）和5＊5＊3超晶胞（图4），与单胞一起作比较分析．

图3 FOX-7的3＊3＊2超晶胞结构Fig．3 3＊3＊2supercell of FOX-7

图4 FOX-7的5＊5＊3超晶胞结构Fig．4 5＊5＊3supercell of FOX-7

2．2 FOX-7晶体热膨胀的分子动力学模拟

MD模拟计算结果的优劣很大程度上也决定于力场的合理选择．FOX-7属于由碳、氢、氧、氮四种原子组成的硝胺类分子晶体，PCFF（The polymer consistent force field）[16－18]及 COMPASS（Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies）[19]对其均应有一定的适用性．通过303K下对单胞与超晶胞的模拟计算（温度外其他参数设置如2．3），对其适用性做了比较．在NPT（等温等压系综）下分子动力学模拟细节如：各分子起始速度按Maxwell分布取样，Velocity Verlet算法[20]进行求解，van der waals力和coulomb作用分别用atom-based和ewald方法[21，22]计算，采用 Anderson法[23]控温，为实现对晶体各向异性形变的模拟，采用Parrinello法[24]控压，外压力取常压10－4GPa，总共进行了200ps的MD模拟以保证体系达到平衡，每500fs取样一次，共得到400帧的轨迹文件，选取平衡后的轨迹文件用于统计分析得到晶胞参数．303K的 MD模拟计算结果及与室温（283K-303K）时CCDC标准晶胞参数值[15]的比较见表1．

表1 303K下FOX-7晶胞参数模拟计算值及其与标准值的偏差Table 1 Simulation values of cell parameters for FOX-7at 303Kand their deviation

首先，将原胞作为初始构型直接进行分子动力学模拟，所得到的晶胞参数模拟值与实验值的偏差很大，超晶胞由于更接近真实体系，所得的数据偏差较小，可以进行进一步比较；其次，总体而言PCFF与COMPASS力场的计算值与实验值偏差都较小，相对地在选择PCFF力场计算时得到的晶胞参数偏差较大，而选择COMPASS力场计算时，晶胞参数计算结果与实验值的相对偏差均在3%以内，其中5＊5＊3超晶胞的棱边夹角α、β、γ相对偏差很小，均在0．5%以内，证明COMPASS力场和5＊5＊3超晶胞模型具有较好的适用性．

此外，晶格能被作为验证力场适用性的能量参数．晶格能可通过升华热的文献值计算，公式如

其中R是标准气体常数，T为热力学温度．FOX-7晶体升华热文献值[25]为108．784kJ／mol，计算所得晶格能为 －113．739kJ／mol；晶格能亦可通过在选定力场下晶体总能及气态分子总能和的差值计算，近似公式为

对FOX-7晶体Z＝4，COMPASS力场下所得晶格能为－109．526kJ／mol，可见两个晶格能值相差小于10kJ／mol，因此可以认为COMPASS力场对FOX-7晶体具有较好的适用性．

选取COMPASS力场，在NPT系综下，以FOX-7的5＊5＊3超晶胞模型为初始结构，其他参数设置不变，分别完成在303K、318K、333K、348K、363K、378K下对FOX-7晶体的分子动力学模拟，得到平衡构象，实现对一定温

度范围内晶体各向异性热膨胀的模拟．

2．3 FOX-7晶体随各向棱长膨胀的能量变化计算

分子晶体的热膨胀性质与晶格能相关．对于FOX-7晶体，其晶体内相互作用越强，则热膨胀系数越小，因此可以通过能量计算来进一步阐释热膨胀各向异性．经典的处理方法难以描述晶体内部各种复杂的相互作用，需要应用量子理论．晶体结合能为晶体总能与自由原子的总能之差，计算晶体总能是关键所在，通过密度泛函理论方法对FOX-7晶体进行计算可得到其能量变化情况．

图5 FOX-7晶体分子动力学模拟的温度／能量－时间平衡曲线Fig．5 Temperature／energy-time equilibrium curves for MD simulation of FOX-7crystal

采用基于密度泛函理论赝势平面波方法的CASTEP程序[26，27]，原子核和电子之间的相互作用用超软赝势[28]来描述，截断能取750eV；虽然LDA（local density approximation）[29，30]是一个相对简单且简化的近似，但它仍能满足某些交换相干能在原理上就原本应该满足具有的求和规则，在计算晶体结构方面获得了出乎意料的成功，而GGA（Generalized Gradient Approximation）[31]考虑了粒子数密度的一阶梯度，但实际计算时可能出现矫枉过正的情况，反而不如LDA的结果好，因此交换关联项采取了LDA近似．通过对沿各棱边方向随0．1%膨胀率（100%-100．5%，该范围接近于303K-378K热膨胀的尺度）的各晶胞模型进行几何构型全优化和总能计算，求算膨胀时各向能量变化率，实现从能量角度阐述FOX-7晶体热膨胀各向异性．

3 结果与讨论

3．1 体系平衡判定

体系达到平衡的标志是能量和温度均趋于平衡，以在303K时为例，平衡时的温度－时间曲线及能量－时间曲线分别如图5所示，后100ps时体系处于平衡，选取相应轨迹文件作为统计分析的原始数据．

3．2 平衡构象分析

图6 303K下FOX-7超晶胞MD模拟后的平衡构象Fig．6 Equilibrium structure of FOX-7supercell after MD simulation at 303K

比较各温度下的平衡构象发现分子在晶体中的取向及堆积方式未见变化．分子内键角的微小扭转将会一定程度上影响分子间相互作用，但相对于温度的影响，这种偏差微乎其微，下面的数据也表明模拟值与实验值的偏差较小，证明了模拟方法较为准确可靠．

3．3 MD计算结果分析

与303K的MD模拟类似，分别对FOX-7晶体在318K、333K、348K、363K、378K下 MD模拟得到的平衡轨迹文件进行统计分析，最后得到晶体在一定温度范围内（303K-378K未发生相变）随温度升高的晶胞参数变化，如图7和表2所示．

表2 不同温度下FOX-7晶体的晶胞参数计算值Table 2 Simulation values of cell parameters for FOX-7 crystal at different temperatures

图7 FOX-7的晶胞棱长与温度的线性关系（303K-378K）Fig．7 Linear relation between cell length and temperature for FOX-7（303K-378K）

图6为303K时MD模拟平衡后FOX-7超晶胞结构，可以看出FOX-7分子的硝基和碳碳双键平面间的扭转角有微小变化但分子堆积方式不变，

可见，模拟得到的晶胞，随温度上升，其棱边夹角α、β、γ偏差和变化不大，棱长a、b、c与温度成良好的线性关系．

3．4 各向热膨胀系数的计算

为研究晶体热膨胀的各向异性，根据微分线膨胀系数计算公式

求算各方向的线膨胀系数，计算结果如表3所示：

可见，随温度升高FOX-7晶体沿三个方向均膨胀；线膨胀系数值与文献值[9，12，13]相近，虽然本文计算值与文献[13]相比较大，但无论文献[9]的实验值还是文献[12]的理论计算值都能与本文结果相互印证；沿b轴方向的线膨胀系数值αb明显较大，与文献[9]中实验值 12．28×10－5K－1及文献[12，13]理论计算值符合很好；而αc与αa的值接近但不相等，能与文献[9]实验值和文献[13]计算值相互符合．由此可见，本文的FOX-7晶体热膨胀值和各向异性特点是可信且较佳的结果．

3．5 热膨胀各向异性的能量分析

在接近于303K-378K温度范围的热膨胀尺度100%-100．5%，采用DFT方法计算得到沿各棱边方向随0．1%膨胀率的晶胞总能、膨胀时各向能量变化及对应的各晶胞棱长，能量值与棱长的关系如图8所示．

随各棱边的不断膨胀，晶胞能量也逐渐增大，膨胀时的能量差值具有物理意义，如表4所示．

表3 FOX-7晶体热膨胀系数值及其与文献值比较Table 3 Comparisons of CTE values for FOX-7crystal between our result and the literature

图8 晶胞总能计算值与棱长的关系（棱长膨胀率0．1%）Fig．8 Correlation between the calculated total energy and the cell length

表4 FOX-7晶体随各向膨胀的总能变化Table 4 Total energy changes for the expansion along each direction for FOX-7crystal

3．6 晶体热膨胀各向异性与分子堆积的关联

首先，FOX-7晶体各向总能量变化与不同膨胀率的棱长成线性关系，即E∝L，而各棱长变化与温度成良好的线性关系，即L∝T，可以推知E∝T，即各向能量变化与温度成线性关系，从而将各向能量变化与温度进行关联．

由沿三个棱边方向的能量（单位：kJ／mol）与棱长（单位：nm）的线性关系，以及棱长与热力学温度（单位：K）的线性关系可以联立得到：

将各组数据代入经整理可得：

再将具体的温度和能量值代回方程（3．2）（3．3）（3．4）便可得到FOX-7晶体热膨胀系数与沿各棱边方向的能量变化的关系式，如303K下

从而建立FOX-7晶体能量变化与热膨胀系数的关联．

4 结 论

和c两方向相互作用不同，因此受热膨胀过程沿各棱边方向的能量变化各不相同，导致b方向热膨胀最为显著，而a和c方向热膨胀系数不同．

[1]Latypov N V，Bergman J，Langlet A，etal．Synthesis and reactions of 1，1-diamino2，2-dinitroeth-ylene[J]．Tetrahedron，1998，54：11525．

[2]Bemm U，Östmark H．1，1-Diamino-2，2-dinitroethylene：a novel energetic material with infinite layers in two dimensions[J]．ActaCrystallogr．，1998，C54：1997．

[3]Hsueh C H，Becher PF．Thermal stresses due to thermal expansion anisotropy in materials with preferred orientation[J]．JournalofMaterialsScience Letters，1991，19（10）：1165．

[4]Boas W，Honeycombe R W K．Theanisotropy of thermal expansion as a cause of deformation in metals and alloys[A]．ProceedingsoftheRoyalSocietyofLondon．SeriesA，MathematicalandPhysicalSciences[C]，（Feb．25），1947，188（1015）：427．

[5]Suhithi M，Chak P，Maija M，etal．Equation of state and structural changes in diaminodinitroethylene from experimental studies and ab initio quantum calculations[A]．12thDetonation（Int．）Symposium[C]．，SanDiego，California，2002．

[6]Kempa P B，Herrmann M，MolinaMetzger F J，et al．Phase transitions of FOX-7studied by X-ray diffraction and thermal analysis[A]．35thInt．Annu．Conf．ICT[C]，2004．

[7]Burnham A K，Weese R K，Wang R，etal．Thermalproperties of FOX-7[A]，36thInternationalAnnualConferenceofICT＆32ndInternationalPyrotechnicsSeminar[C]，Karlsruhe，Germany，June 28，2005through July 1，2005．

[8]Burnham A K，Weese R K，Wang R，etal．Solidsolid phase transition kinetics of FOX-7[A]，2005 NATASAnnualConference[C]．UniversalCity，CA，UnitedStates，September 18-21，2005．

[9]Evers J，Klapötke T M，Mayer P，etal．α-andβ－FOX-7，polymorphs of a high energy density material，studied by X-ray single crystal and powder investigations in the temperature range from 200to 423K[J]．Inorganicchemistry，2006，45（13）：4996．

[10]Zhao J J，Liu H．High-pressure behavior of crystalline FOX-7by density functional theory calculations[J]．ComputationalMaterialsScience，2008，42：698．

[11]Hu A G，Larade B，Hakima A R，etal．A first principles density functional study of crystalline FOX-7chemical decomposition process under external pressure[J]．Propellants，Explosives，Pyrotechnics，2006，31（5）：355．

[12]Sorescu D C，Boatz J A，Thompson D L．Classical and quantum-mechanical studies of crystalline FOX-7（1，1-diamino-2，2-dinitroethylene）[J]．Journal ofPhysicalChemistryA，2001，105（20）：5010．

[13]Taylor D E，Rob F R，Betsy M，etal．A molecular dynamics study of 1，1-diamino-2，2-dinitroethylene（FOX-7）crystal using a symmetry adapted perturbation theory-based intermolecular force field[J]．PhysicalChemistryChemicalPhysics，2011，13（37）：16629．

[14]Bian L，Shu Y J，Li H R．Computational investigation on the surface electronic density，morphology and detonation property of 1，1-diamino-2，2-dinitroethylene（FOX-7）crystal[J]．ChineseJ．Struct．Chem．，2012，31（12）：1736．

[15]Gilardi R．Cambridge crystallographic data centre，CCDC127539，1999．

[16]Sun H，Mumby S J，Maple J R，etal．An ab initio CFF93all-atom forcefield for polycarbonates[J]．J．Am．Chem．Soc．，1994，116：2978．

[17]Sun H．Ab initio calculations and forcefield development for computer simulation of polysilanes[J]．Macromolecules，1995，28：701．

[18]Sun H，Rigby D．Polysiloxanes：ab initio forcefield and structural，conformational and thermophysical properties[J]．Spectrochim．Acta，PartA，1997，53：1301．

[19]Sun H．COMPASS：AnabInitioforce-field optimized for condensed-phase applications[J]．J．Phys．Chem．B．，1998，102（38）：7338．

[20]Verlet L．Computer experiments on classical fluids．I．thermodynamical properties of Lennard-Jones molecules[J]，Phys．Rev．，1967，159：98．

[21]Ewald P P．Die Berechnung optischer und elektrostatischer Gitterpotentiale[J]．Ann．Phys．Leipzig，1921，64：253．

[22]Karasawa N，Goddard W A．Acceleration of convergence for lattice sums[J]．J．Phys．Chem．，1989，93：7320．

[23]Andersen H C．Molecular dynamics simulations at constant pressure and／or temperature[J]．J．Phys．Chem．，1980，72：2384．

[24]Parrinello M，Rahman A．Polymorphic transitions in single crystals：a new molecular dynamics method[J]．J．Appl．Phys．，1981，52：7182．

[25]Politzer P，Concha M C，Grice M E，etal．Computational investigation of the structures and relative stabilities of amino／nitro derivatives of ethylene[J]．J．Mol．Struct．（THEOCHEM），1998，452：75．

[26]Segall M D，Lindan P J D，Probert M J，etal．First-principles simulation：ideas，illustrations，and the CASTEP code[J]．J．Phys．：CondensedMater，2002，14：2717．

[27]Clark S J，Segall M D，Pickard C J，etal．First principles methods using CASTEP[J]．Zeitschrift fuerKristallographie，2005，220（5-6）：567．

[28]Vanderbilt D．Soft self-consistent pseudopotentials in a generalized eigenvalue formalism[J]．Phys．Rev．B，1990，41：7892．

[29]Ceperley D M，Alder B J．Ground state of the electron gas by a stochastic method[J]．Phys．Rev．Lett．，1980，45：566．

[30]Perdew J P，Zunger A．Self-interaction correction to density-functional approximations for many-electron systems[J]．Phys．Rev．B，1981，23：5048．

[31]Perdew J P，Burke K，Ernzerhof M．Generalized gradient approximation made simple[J]，Phys．Rev．Lett．，1996，77：3865．