Anisotropy of Thermal-expansion for β-Octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine:Quantum Chemistry Calculation and Molecular Dynamics Simulation

Wen Qian,Chao-yang Zhang,Yuan-jie Shu,Ying Xiong,He-hou Zong,Wei-bin Zhang

Institute of Chemical Materials,China Academy of Engineering Physics,Mianyang 621900,China

Anisotropy of Thermal-expansion for β-Octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine:Quantum Chemistry Calculation and Molecular Dynamics Simulation

Wen Qian,Chao-yang Zhang,Yuan-jie Shu?,Ying Xiong,He-hou Zong,Wei-bin Zhang

Institute of Chemical Materials,China Academy of Engineering Physics,Mianyang 621900,China

Molecular dynamics simulations on octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine(HMX) at 303-383 K and atmospheric pressure are carried out under NPT ensemble and COMPASS force f i eld,the equilibrium structures at elevated temperatures were obtained and showed that the stacking style of molecules don’t change.The coefficient of thermal expansion (CTE)values were calculated by linear fitting method.The results show that the CTE values are close to the experimental results and show anisotropy.The total energies of HMX cells with separately increasing expansion rates(100%-105%)along each crystallographic axis was calculated by periodic density functional theory method,the results of the energy change rates are anisotropic,and the correlation equations of energy change-CTE values are established.Thus the hypostasis of the anisotropy of HMX crystal’s thermal expansion,the determinate molecular packing style,is elucidated.

Octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine,Molecular dynamics simulation,Thermal expansion,Anisotropy,Density functional theory

I.INTRODUCTION

Octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine, commonly known as HMX or octogen,is a kind of powerful and relatively insensitive nitroamine explosive [1].HMX has high density,high energy,as well as good thermal stability,so that it’s widely applied in ammunitions or warheads.HMX crystal is polymorphous but only β-HMX is stable at room temperature and is actually used[2,3].HMX is hard to be independently molded,thus HMX crystals are often mixed with additives such as binders,plasticizers,and solidif i ed agents to form HMX-based PBX(polymer bonded explosive)[1].For weapon utility,HMX-based PBX needs to be resistant to thermal impact;however,the anisotropic thermal expansion of explosive crystals will bring notable internal stress,and micro cracks will appear inside the explosive part when the effect of internal stress is bigger than that of cohesive energy[4], and these internal micro voids could induce hot spot and increase the sensitivity of the energetic materials [5];moreover,the anisotropic thermal expansion of the crystals will induce thermal deformation of PBX,make dimension or density of the explosive part irreversibly change,thus influence the measurement precision[6]. On all accounts,a deep research on the anisotropy of thermal expansion of HMX is of great importance.

The properties of crystalline HMX has been studied experimentally in many ways.The molecular structures of α,β and γ-HMX were gotten from IR spectrum, while the crystal densities and dielectric constants were obtained by other measurements[7].CTE values of the β and δ-HMX were measured by thermo-mechanical analysis(TMA),and the phase transition kinetics was measured by differential scanning calorimetry(DSC) [10].The phase transition and thermal expansion of HMX crystal were studied by X-ray diffraction(XRD) technique[8,9],then more precise CTE values were obtained by ref i nements of the cell parameters data from XRD[11,12,13].Phase transition process of HMX crystal was also observed by high-temperature atomic force microscope(AFM)[14].In spite of this,in order to investigate the hypostasis of HMX crystal properties, the experimental measurements are not enough and further theoretical studies are needed.

With the progress of physics,chemistry,and information technology,more and more researchers used molecular simulation and computational methods to theoretically study the properties of HMX crystal.Isothermalisobaric molecular dynamics simulations(NPT-MD)of α,β,and γ-HMX were performed under a modif i ed RDX intermolecular potential,and the space group symmetries and structural parameters for the three phases were predicted[15].Thermal expansion and sublimation enthalpies for all three pure polymorphs of HMX have been investigated and compared to experimental results by molecular dynamics(MD)sim-ulation using a fl exible molecule force field[16].Condensed phase density functional theory(DFT)was used to study various uniaxial compressions of β-HMX and extract energy versus extension information[17].NPTMD have been applied to study the structural,vibrational and thermodynamic properties of β-HMX, the initial decomposition and pressure-induced phase transition were discussed[18].First-principle planewave method with an ultrasoft pseudo potential and generalized gradient approximation was used to calculate the lattice parameters,thermodynamics and electronic properties of β-HMX[19].Self-consistent charge density-functional tight-binding method(DFTB)and MD simulation were used to study the initial chemical processes of condensed-phase HMX under shock wave loading[20].ab initio MD simulations combined with multi-scale shock technique have been performed to investigate the initial chemical processes of HMX under shock wave loading,a new mechanism for initiation decomposition was suggested[21].Compressive shear reactive dyanmics(CS-RD)model with ReaxFF-lg force field has been used to predict the anisotropy of shock sensitivity and chemical process initiation in β-HMX [22].MD simulations of nanoparticles of HMX crystal polymorphs were used to study the interaction potential energy and heat of sublimation[23].

In this work,MD simulations combined with fi rstprinciple calculations are used to study the anisotropy of HMX crystal’s thermal expansion,the methods and results are helpful in elucidating the anisotropic thermal expansion of HMX crystal at atmospheric pressure.

II.SIMULATION AND CALCULATION METHOD

A.Structure modeling

Modeling process was carried out in Materials Studio software package[24].The molecular structure is shown in Fig.1(a),it can be seen that HMX has a simple molecular structure that consists of an eight-membered ring formed by alternating carbon and nitrogen atoms,with a nitro group attached to each nitrogen atom which brings in high energy.

HMX crystal structure in Fig.1(b)is based on the neutron diffraction data of β-HMX crystal[25], β-HMX crystal belongs to monoclinic system and P21/c space group,with the cell parameters as:a=6.540?A, b=11.050?A,c=8.700?A,β=124.3?,and α=γ=90?.

To approach reality and obtain more precise results,a super cell based on the HMX primitive cell needs to be constructed as initial configuration of MD simulation. In this work,a 4×2×3 super cell containing 1344 atoms is constructed as shown in Fig.1(c).

B.MD simulation

FIG.1 Structures of HMX.(a)Molecular structure,(b) crystal structure,(c)super cell through a-axis,and(d)super cell through c-axis.

Proper selection of the force f i eld is the key to precise MD simulation results.PCFF(polymer consistent force field)[26-28]is one of the most appropriate force fields for MM(molecular mechanics)and MD simulations.COMPASS(condensed-phase optimized molecular potentials for atomistic simulation studies)[29]is an ab initio force field especially applicable to condensedphase compound.HMX is a kind of nitro-amino compound.PCFF and COMPASS are both thought to be applicable to HMX molecular crystal.

Primitively the comparison of these two force fields is done by comparison of the simulation results for cell parameters at certain condition.The simulation details are as follows:MD simulations for HMX crystal are carried out under NPT(with constant particle number, pressure and temperature)ensemble;and velocity verlet arithmetic[30]is used;the initial velocity is sampled by Maxwell distribution;van der Waals force is calculated by atom-based method and coulomb interaction is calculated by Ewald method[31,32];the target temperature and pressure are set to 303 K and 0.1 MPa, the temperature is controlled by Anderson method[33] and the pressure is controlled by Parrinello method[34]; the MD simulation lasts 200 ps in total,a time step 1 fs is used;and equilibrium of the system is determined both by energy and temperature.The equilibrium trajectory documents are selected and the cell parameter values are gained from them.The simulation results of cell parameters for HMX crystal using di ff erent systems and force fields are shown in Table I.

Firstly,in the MD simulation on primitive cell,the system is not so easy to equilibrate,and the result of cell parameters is not so precise.On the other hand,the MD simulation on HMX super cell can be well equili-brated and the result is much closer to the experimental value[12].Hence 4×2×3 super cell is used in the coming MD simulations.

TABLEI Simulation values of cell parameters for HMX crystal at 303 K,the relative deviation is in parentheses.

Secondly,for HMX super cell,compared with experimental values,relative deviations of cell parameters are much larger when choosing PCFF,while relative deviations of the unit cell edge lengths are less than or equal to 5%and that of inter-axial angles are less than 0.5% when choosing COMPASS,which proves that COMPASS force f i eld is of better applicability.

For the force f i eld validation,it is also standard to use an energetic component,so that calculation of sublimation energy is carried out to verify the system as well as the force f i eld,and the simplif i ed equations are as follows

where Ecrystaland Emoleculeis the total energy of the crystal and that of a molecule respectively.The number of molecules contained in each cell n is 2 for β-HMX crystal.

The result from the simulation using COMPASS force f i eld shows that Elattice=-193.794 kJ/mol and?H=188.839 kJ/mol.Because the experimental value of the sublimation energy is 184.766 kJ/mol[35],so that the sublimation energy calculated from the simulation result matches with the experimental value (?＜10 kJ/mol).Hence the COMPASS force f i eld is most fit for the MD simulation.

MD simulations are carried out on 4×2×3 super cell for HMX crystal under NPT ensemble using COMPASS force f i eld,periodic boundary condition is used in all simulations.To simulate thermal expansion process under atmospheric condition,the target pressure is set to 0.1 MPa and the temperatures are elevated from 303 K to 378 K by 20 K each step;Anderson method is used to control the temperature and Parrinello method is used to control the pressure.Initial velocity is sampled by Maxwell distribution,and the time step is 1 fs in each case;each MD simulation lasts 100 ps to ensure the system properly equilibrated both in energy and temperature,another 100 ps time of simulation is followed to get abundant data for subsequent statistical analysis.The equilibrium trajectory documents are selected to do analysis about the thermal expansion properties.

C.DFT calculation

As a kind of molecular crystal,the stronger the interactions inside HMX crystal are,the larger the binding forces will be,and the CTE values are thought to be smaller,so that it’s efficient to elucidate the anisotropy of HMX crystal’s thermal expansion by energy calculations.The classic calculation methods are hard to describe all kinds of complex interactions inside the crystal,and quantum chemistry is introduced to solve the problem.Binding energy is the difference between total energies of the crystal and that of the free particles,thus calculation of the total energy of the crystal is the key to examine the energy change.In this work,periodic DFT calculations on HMX crystal are carried out to get the energy change along each crystallographic axis,the details are as follows.

The primitive cell for HMX crystal is used as input structure.CASTEP[36]program,based on the DFT method using plane-wave basis set and ultra-soft pseudo-potentials(USP)[37],is used to do energy calculations on the crystal structures.The local density approximation function(LDA),which was f i rst established by Ceperley and Alder[38]and later parameterized by Perdew and Zunder[39](CA-PZ),is used to describe the exchange correlation terms.Geometry optimizations using LDA/CA-PZ are carried out on HMX cells with three unit cell edge length expanding from 100%to 105%respectively,and the total energy for each cell is calculated one by one.Then the correlation between total energy of crystal cell and the axial length can be revealed.

III.RESULTS AND DISCUSSION

A.Equilibrium structure after MD simulation

FIG.2 Equilibrium structure of HMX super cell after MD simulation at 303 K through(a)a-axis and(b)c-axis.

From the equilibrium structure graphs(Fig.2)of HMX super cell after MD simulation at 303 K along three different directions we can see that both the tropism of the molecules and the way how molecules stack in the crystal stay the same,which also prove the validity of COMPASS and the MD simulation method.

B.Results of the MD simulations

The equilibrium trajectory documents from MD simulations on HMX super cell at elevated temperature from 303 K to 383 K are analyzed,and the cell parameters obtained are listed in Fig.3 and Table II.

It can be seen that the inter-axial angles(α,β,γ) don’t change much and the relative deviations are all below 0.5%,which prove the validity of the simulation results.The unit cell edge lengths(a,b,c)are linearly related to temperature(T)as shown in Fig.3, the linear CTE values are calculated to characterize the anisotropic thermal expansion,and the equation for linear CTE is as follows,

the results of simulated linear CTE values and the comparison between our results and experimental values from Ref.[11-13]are shown in Table III.

It can be seen that HMX crystal expands along three axial directions as the temperature increases,and the CTE values are notably different along each crystallographic axial direction(δb＞δa＞δc),showing the anisotropy of the thermal expansion.

Due to the differences between the referenced CTE results from different measurement techniques,we think that the calculation results can act as a criterion of the diverse results,so that we mainly discuss about the similarities.Compared with experimental values using TMA technique[10](about 3.7×10-5K-1near room temperature and 8×10-5-10×10-5K-1at higher temperature),values from single-crystal XRD experiment are more precise and can show anisotropic CTEs in different directions which are more meaningful in comparison with our simulation results.As shown in Table IV, all the results show obvious anisotropy in the thermal expansion process,and for most of them,δb＞δa＞δc, which is the same as that in our results,while the CTE values have slight difference.The linear CTE values for b-axis is very close to Saw’s experimental value(LLNL) [11],which was obtained by ref i nements of the cell parameters from XRD data at elevated temperature(heating from 203 K to 323 K),but for a-axis and c-axis,our results seem a little larger than Xue’s[12],which were also gained through XRD experimentally.Besides,for the latest experimental values from Deschamps(Naval Research Laboratory)[13],δcis closer to our result.

FIG.3 The linear relationship between the unit cell edge length and temperature for HMX crystal.

TABLE II Relative deviation of simulation results of interaxial angles for HMX crystal.

TABLE III Comparison between simulation CTE results and experimental values.

FIG.4 Linear relationship between the total energy and the axial length of HMX crystal.

C.Results of the energy calculations

Calculations of the total energies of the cells with axial length respectively expanding from 100%to 105%are done using LDA/CA-PZ method,thus the anisotropy of thermal expansion of HMX crystal can be investigated in the point view of energy,and the results are shown in Fig.4.It can be seen that the total energy increase with the expansion of each crystallographic axial length,and the energy is found to be linearly related to the axial length.Because of the approximation algorithm used in the energy calculation,the absolute values of the total energy make no sense,while the difference between the initial energy and the energy after expansion(shown in Table IV)is more physically meaningful.

TheenergychangeratesalongeachaxialdirectionofHMXcrystalaredEa/dLa=674.3890, dEb/dLb=595.9972,dEc/dLc=733.4319,respectively. When each crystallographic axis use the same expansion rate,the energy changes are different;the energy change rate during expansion is anisotropic (dEb/dLb＜dEa/dLa＜dEc/dLc),which is related to the anisotropy of CTE(δb＞δa＞δc).As a matter of fact, the bigger energy change rate during expansion is,the smaller the CTE value will be.

D.Discussions about energy change and anisotropic thermal-expansion

It can be concluded that in HMX crystal’s thermal expansion process,the energy change is linearly related to the expanding axial length(E∝L),and the length change is linearly related to the temperature(L∝T),accordingly it can be inferred that energy change along each crystallographic axis is linearly related to the temperature(E∝T).

TABLE IV Total energies changes with increasing expansion rates.

The linear relationship between total energy of HMX crystal and axial length(in Fig.4)and that between axial length and temperature(in Fig.3)can be summed together as follows:

and the linear relationship between total energy and temperature can be easily extrapolated as

It can be seen that the energy change rates along each axial direction of HMX crystal during heating process are dEa/dT=0.0300,dEb/dT=0.0948,and dEc/dT=0.0246,respectively.Obviously the energies increase with temperature slightly,which can explain the slight decrease of CTE values at increasing temperature;what’s more important,it can be validated that the anisotropy of thermal expansion is related to the difference of energy change needed along each axis during heating(dEb/dT＞dEa/dT＞dEc/dT).

And the correlation of energy change and CTE can be obtained by substituting the temperature and energy value back to Eqs.(4)-(6),for example the equations at 303 K are

As a matter of fact,energy calculation is the tool to characterize intermolecular interactions in the crystal.HMX molecules stack in a determinate way in HMX molecular crystal,causing different situations of intermolecular interactions along different axial directions; as a result,the energy change needed in thermal expansion is also different along different axial directions.

IV.CONCLUSION

MD simulations on HMX super cell under NPT ensemble using COMPASS force f i eld with periodic boundary condition are successfully carried out to get the CTE values of HMX crystal at 303-383 K and atmospheric pressure,and DFT calculations on HMX crystal are done to gain the energy changes during expansion of each crystallographic axial length from 100%to 105%respectively.The simulation and calculation results show that:(i)The thermal expansion is anisotropic(δb＞δa＞δc)and CTE values are close to the experimental results in literature.(ii)HMX molecules stack in a determinate way in the crystal,both the tropism of the molecules and the way how molecules stack in the crystal don’t change during thermal expansion process.(iii)When each axis use the same expansion rate,the energy change rate along each axial direction is anisotropic(dEb/dLb＜dEa/dLa＜dEc/dLc). (iv)Accordingly it can be concluded that the energy change along each crystallographic axis is linearly related to the temperature,and the anisotropy of thermal expansion of HMX crystal is related to the difference of energy change needed in thermal expansion along each axis(dEb/dT＞dEa/dT＞dEc/dT).(v)Thus the hypostasis of the anisotropy of HMX crystal’s thermal expansion can be elucidated as:because of the determinate stack style of molecules in HMX molecular crystal,the intermolecular interactions are anisotropic along three axial directions,inducing energy change needed in thermal expansion to be anisotropic.

V.ACKNOWLEDGMENTS

This work was supported by the Director’s Foundation from Institute of Chemical Materials,China Academy of Engineering Physics(No.626010948).

[1]J.Akhavan,The Chemistry of Explosives,Cambridge: The Royal Society of Chemistry(2004).

[2]W.C.McCrone,Anal.Chem.22,1225(1950).

[3]H.H.Cady,A.C.Larson,and D.T.Cromer,Acta Crystallogr.16,617(1963).

[4]C.H.Hsueh and P.F.Becher,J.Mater.Sci.Lett.19, 1165(1991).

[5]S.D.McGrane,A.Grieco,K.J.Ramos,D.E.Hooks, and D.S.Moore,J.Appl.Phys.105,7(2009).

[6]W.Boas and R.W.K.Honeycombe,Proc.R.Soc.London,Ser.A 188,427(1947).

[7]M.Bederd,H.Huber,J.L.Myers,and G.Wright,Can. J.Chem.40,2278(1962).

[8]M.Hermann,W.Engel,and N.Eisenreich,Propellants. Explos.Pyrotech.17,190(1992).

[9]M.Hermann,W.Engel,and N.Eisenreich,Zeitschrifffur Kristallographie 204,121(1993).

[10]R.K.Weese and A.K.Burnham,Propellants.Explos. Pyrotech.30,344(2005).

[11]C.K.Saw,Kinetics of HMX and Phase Transitions: effects of Grain Size at Elevated Temperature,UCRLJC-145228,Lawrence Livermore National Laboratory (LLNL),June 13(2002).

[12]C.Xue,J.Sun,B.Kang,Y.Liu,X.F.Liu,G.B.Song, and Q.B.Xue,Propellants.Explos.Pyrotech.35,333 (2010).

[13]J.Deschamps,M.Frish,and D.Parrish,J.Chem.Crystallogr.41,966(2011).

[14]B.L.Weeks,C.M.Ruddle,J.M.Zaug,and D.J.Cook, Ultramicroscopy 93,19(2002).

[15]D.C.Sorescu,B.M.Rice,and D.L.Thompson,J. Phys.Chem.B 102,6692(1998).

[16]D.Bedrov,C.Ayyagari,G.D.Smith,and T.D.Sewell, J.Comput.Aided Mater.Des.8,77(2001).

[17]M.W.Conroy,I.I.Oleynik,S.V.Zybin,and C.T. White,J.Appl.Phys.104,053506(2008).

[18]L.Y.Lu,D.Q.Wei,X.R.Chen,G.F.Ji,X.J.Wang, J.Chang,Q.M.Zhang,and Z.Z.Gong,Mol.Phys. 107,2373(2009).

[19]H.L.Cui,G.F.Ji,X.R.Chen,W.H.Zhu,F.Zhao, Y.Wen,and D.Q.Wei,J.Phys.Chem.A 114,1082 (2010).

[20]N.N.Ge,Y.K.Wei,G.F.Ji,X.R.Chen,F.Zhao, and D.Q.Wei,J.Phys.Chem.B 116,13696(2012).

[21]W.H.Zhu,H.Huang,H.J.Huang,and H.M.Xiao,J. Chem.Phys.136,044516(2012).

[22]T.T.Zhou,S.V.Zybin,Y.Liu,F.L.Huang,and W. A.Goddard,J.Appl.Phys.111,124904(2012).

[23]H.Akkbarzade,G.A.Parsafar,and Y.Bayat,Appl. Surf.Sci.258,2226(2012).

[24]Materials Studio 5.0,,San Diego:Accelys Inc.,CA (2009).

[25]C.S.Choi and H.P.Boutin,Acta Crystallogr.,Sect.B 26,1235(1970).

[26]H.Sun,S.J.Mumby,J.R.Maple,and A.T.Hagler, J.Am.Chem.Soc.116,2978(1994).

[27]H.Sun,Macromolecules 28,701(1995).

[28]H.Sun and D.Rigby,Spectrochim.Acta Part A 53, 1301(1997).

[29]H.Sun,J.Phys.Chem.B 102,7338(1998).

[30]L.Verlet,Phys.Rev.159,98(1967).

[31]P.P.Ewald,Ann.Phys.Leipzig 64,253(1921).

[32]N.Karasawa and W.A.Goddard,J.Phys.Chem.93, 7320(1989).

[33]H.C.Andersen,J.Phys.Chem.72,2384(1980).

[34]M.Parrinello and A.Rahman,J.Appl.Phys.52,7182 (1981).

[35]J.M.Rosen and C.Dickinson,J.Chem.Eng.Data 14, 120(1969).

[36]M.D.Segall,P.J.D.Lindan,M.J.Probert,C.J. Pickard,P.J.Hasnip,S.J.Clark,and M.C.Payne,J. Phys.:Condens.Matter 14,2717(2002).

[37]D.Vanderbilt,Phys.Rev.B 41,7892(1990).

[38]D.M.Ceperley and B.J.Alder,Phys.Rev.Lett.45: 566(1980).

[39]J.P.Perdew and A.Zunger,Phys.Rev.B 23,5048 (1981).

ceived on July 9,2013;Accepted on September 18,2013)

?Author to whom correspondence should be addressed.E-mail：syjfree@sina.com