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．