• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Molecular Dynamics Simulations of Crystalline δ-HMX with Void Defect

    2015-05-10 05:43:16LIAONingLIWenpengZHOUXiaoqingDUANXiaohui
    含能材料 2015年12期

    LIAO Ning, LI Wen-peng, ZHOU Xiao-qing, DUAN Xiao-hui

    (1. State Key Laboratory Cultivation Base for Nonmetal Composites and Functional Materials, Southwest University of Science and Technology, Mianyang 621010, China; 2. Institute of Chemical Materials, China Academy of Engineering Physics, Mianyang 621999, China)

    1 Introduction

    As an important energetic material,octahydro-1, 3, 5, 7-tetranitro-1, 3, 5, 7-tetrazocine (C4H8N8O8, HMX) has four known crystalline phases, namely,α,β,δandγ. It has been observed that theδform has 4 coplanar carbon atoms with aC2axis;α-HMX has aC2axis, andβ-HMX has a center of inversion symmetry[1]. The stabilities of the crystalline forms at room temperature areβ>α>γ>δ, the same as the order of the densities[2]. When heated to temperatures above 435 K, the monoclinicβ-HMX converts to the hexagonalδmorphology[3-4], the most reactive phase[5-6]. For the transformation fromβphase toδphase, it involves major changes of both the crystal lattice and the molecular ring conformation[3-6]. Accompanying the transition process is that a significant lattice expansion[7]results in an appearance of cracks and large amount of hot spots in the material and the fast growth of reaction during shock compression[8]. Those can be contributed to sensitivity increase of theδphase[3-4,8]. The decomposition of theδ-HMX practically coincides with theβδphase transition and the possible decomposition of theβ-phase.

    Many experimental reports declare a high sensitivity ofδphase, atomistic theoretical investigations are scarce[9-11]. Many factors, such as density, defects, electronic excitations, and particle sizes, suggest to contribute to the sensitivity of energetic materials. Among them, defects including voids, molecular vacancies, dislocations, surface, and interfaces etc. play a crucial role to initiate explosive decomposition reactions, referred to as sensitivity. At present, the investigations on defect effects mainly focus on theβ-HMX[12-14], and only a literature had studied the surface-accelerated decomposition ofδ-HMX[15]. These investigations have indicated that the surfaces, interfaces, voids, or vacancies can lower the activation barriers and accelerate kinetics of decomposition of HMX. Seeing that the special function ofδphase in the initial reaction of decomposition, dynamics characterization ofδ-HMX containing defects are very important to further understand the reaction mechanism of HMX.

    In our previous work[14], we investigated the void defects contained in crystallineβ-HMX. Void collapse observed with runtime was promoted by the change of the strain field and the increased system pressure due to applying canonical ensemble (NVT). In this paper, aimed at gaining a better understanding of molecular dynamics characterization of void defects, we further researched the void defects in crystallineδ-HMX applying NPT (normal pressure and temperature) simulations, including the activation barriers for an isolated molecule, void evolvement, and molecular conformations change along with simulation time. This investigation may provide some atomistic details for the observed sensitivity increase of theδphase compared to theβphase and the degradation mechanism of HMX.

    2 Methodology

    The structure of crystallineδ-HMX is hexagonal and belongs to the space group P61. The unit cell of the crystal is shown in Fig.1. The initial structure was created using the crystallographic data from Ref. 4. In theδ-phase lattice, each C4H8N8O8molecule has a boat conformation, having all four NO2groups on the same side of the C4N4ring (Fig.1). According to neutron diffraction experiments ofβ-HMX[16], we constructed the unit cell model presented in Fig.1.β-lattice belongs to monoclinic space groupP21/cand contains 2 independent HMX molecules in unit cell. Eachβ-HMX molecule presents chair form and has a symmetry center. Theβ-HMX molecule with central symmetry has lower energy thanδ-HMX molecule withC2symmetry. The super-cells used for simulation of void defects are composed of 5×5×2 unit cells forδcrystal and 5×5×6 forβone (300 molecules, 8400 atoms). After that, cohesive voids were introduced by removing 30 molecules from the 300-molecules super-cell, first one and its neighbors over an increasing distance; this corresponds to a 10% concentration of vacancies in the material.

    Fig.1 Unit cell structures forβandδphases and HMX molecular conformations inδandβphases. Legend: carbon-gray, hydrogen-white, nitrogen-blue, oxygen-red

    The active energy barriers were determined by the QST3 method, which searched for the transition structure using the Synchronous Transit-Guided Quasi-Newton (STQN) method. This method was developed by H. B. Schlegel and coworkers[17-18], and it uses a linear synchronous transit or quadratic synchronous transit approach to get closer to the quadratic region around the transition state and then uses a quasi-Newton or eigenvector-following algorithm to complete the optimization. Density function theory (DFT) method combined with 6-31++G (d, p) basis set and B3LYP functional was applied. Convergence criterion was set to be “tight”. Firstly, the geometry optimization was performed at the B3LYP/6-31++G (d, p) level. Secondly, the search for the transition state was carried out by QST3 method. Finally, the harmonic vibrational frequencies were calculated based on the optimized geometries. All quantum chemical calculations were completed through the Gaussian software[19].

    The dynamics characterization of void defectwas investigated by a MD simulation. A large number of literatures have indicated that the COMPASS (condensed-phase optimized molecular potentials for atomistic simulation studies) force field[20]can be effectively applied to the theoretical simulations for the crystalline HMX[14, 21-24]. Therefore, COMPASS force field was used in our simulations. Experimental,δ-HMX is stable in the temperature range of 433-553 K[25]andβ-HMX under room temperature at an atmospheric pressure[2, 25-26]. Besides, the melting point ofβ-HMX is in the range of 540-550 K due to different determined methods. So, the simulation temperatures were set to be 200, 300 K and 500 K, containing the stable and unstable temperatures for two forms. The Coulombic and van der Waals interactions were calculated by the standard Ewald method with quality “fine”[27]. The charge was calculated using the QEd method[28]. The basis of this method is the equilibration of atomic electrostatic potentials with respect to a local charge distribution. The NPT ensemble was applied with a time step of 1 fs. Temperature was controlled using the Nose thermostat and pressure by Berendsen method. The system was simulated over a total time of 1ns and the trajectory was output every 5000 time steps. The simulation was run with the commercial molecular modeling software package Materials Studio 3.0[29].

    3 Results and Discussion

    3.1 Active Barriers of Conformational Transitions in Vacuum

    Four unique low-energy conformers were found for HMX molecule in vacuum, two whose conformational geometries correspond closely to those found in HMX polymorphs,βandα, and two additional, lower energy conformers that are not seen in the crystalline phases, boat-chair (BC) and boat-boat (BB)[30]. The relative energies of all conformers of an isolated molecule can be found in literature [30], namely,α>β>BB>BC. These conformations are known to exist within a narrow interval of energies.

    In this work, the energies of transition states for the conformational transitions were obtained by the QST3 method at the level of B3LYP/6-31++G(d, p). The calculated active barriers were listed in Table 1 with the correction of quantum zero-point vibrational effects. The frequency calculations indicate that there is one and only one imaginary frequency for all the transition structures, and the four conformers, namely,α,β, BB, and BC, have no any imaginary frequencies. Therefore, we can determine the structures obtained by the QST3 method corresponding to the transition states.

    From the active barriers listed in Table 1, we can find that the transitions from other conformations to BC have lower active barriers. Specially, the barrier fromαto BC is negative, which means this conversion is very favorable in energy. In all possible transitions, the value of the highest energy barrier for theβBB conversion is 33.22 kJ·mol-1. Secondly, the barrier fromβtransits to α is 32.72 kJ·mol-1. The two barriers are approximately 7 times higher thanRTat 500 K. Other barriers are much lower than the two values. Based on the Arrhenius equation, we predict that all conversions may easily take place thermodynamically.

    Table 1 Active barriers of the conformation transitions (in kJ·mol-1)

    transitionsβ?αβ?BBβ?BCα?BBα?BCBB?BCforward32.7233.2210.3321.46-3.644.85backward21.4625.569.5825.565.449.58

    3.2 Dynamics Evolvement of Void

    For the MD simulations at different temperatures,we have observed the different evolvement trends of voids inserted into two different forms,δandβ(see Fig.2). Forδpolymorph, the void has been completely collapsed at three temperatures. At 200 K and 300 K, the molecules collapsed into the void present the characteristics of partial order, but when the temperature increases up to 500 K, not only the molecules near the void but also the whole lattice melt into the liquid state. This can be further verified by the mean square displacement (MSD) analysis later (see Fig.3). Under low temperatures of 200 K and 300 K, except for the molecules located at or near the surface of a void, the whole lattice has also large change, losing strictly periodic structure seen from Fig.2. However, forβphase, at the temperatures of 200 K and 300 K, the void is still hold on, only the shrink of void and some shift of molecules toward the center of void occur due to the change of the strain field. When the temperature is up to 500 K, the case is the same as theδphase, namely, the whole lattice completely collapsed into liquid state.

    Compared to initial model, the density and volume of system have changed after 1ns MD simulations. Under low temperatures of 200 K and 300 K, the lattice shrink ofδ-HMX results in the decrease of volume and the increase of density. For example, the density increases from 1.67 g·cm-3to 1.70 g·cm-3at 300 K and 1.73 g·cm-3at 200 K. The volume change accompanying void collapse is 1.53% and 3.20%, respectively. However, when the simulated temperature goes to 500 K, the lattice expands approximately 10.5%. Accordingly, the density reduces down to 1.51 g·cm-3. Different fromδ-HMX, the volume ofβ-HMX is on the increase with the simulated temperature. When the temperature changes from 200 K to 500 K, the volume grows from 2.1% to 16.4%. At 500 K, despite ofβorδphase, the system completely melts into liquid state and has the same density -1.15 g·cm-3.

    The different evolvement trends result from the different lattice arrangements and intermolecular interactions. The experimental value for heat of sublimationHsubis 184.76 kJ·mol-1and 175.90 kJ·mol-1for theβandδphases, respectively[31-32]. Structural stability and sublimation enthalpy are very useful for understanding the safety of HMX in different crystal phases. More sublimation enthalpy means that the system needs more energy to overcome the lattice bounding and inter- and intra- interactions when the solid energetic materials convert into small gaseous molecules. Besides, the formation energy of a void defect with sizen,Ef(n), is also illustrated the different evolvement trends ofδandβphases.Ef(n) was calculated according to the method in Ref [33]. The value forEf(30) in crystallineδ-HMX is calculated to be 6833.35 kJ·mol-1, and the formation energy per molecule removed,Ef(n)/n, is 227.78 kJ·mol-1. Applying the same method, this value ofEf(30) forβ-HMX is 7009.66 kJ·mol-1,Ef(30)/30 equals to 233.66 kJ·mol-1. These quantitative values indicate the difference in the crystal lattice, the molecular conformation, and the microstructure of two crystalline phases. At the same time, it has also illustrated thatβ-HMX is more stable thanδphase. Our MD simulations show that theβ-HMX containing a void defect is still more stable thanδ-HMX with the same size void. Generally, the so-called “hot spots” in energetic materials are mainly the local regions of a crystal (crystal defects or deformations). Voids, molecular vacancies, dislocations, pore, impurities, and other types of defects play a crucial role in initial reactions. According to hot spots theory, the initiation of combustion in energetic materials is associated with “hot spots”, in which collapse of voids may be one of the most important mechanisms the leads to ignition in pressed explosives. So, the different evolvement trends also give some atomistic details for the observed sensitivity increase of theδphase compared to theβphase and the initial reaction mechanism of HMX.

    a.δ, 200 K b.δ, 300 K c.δ, 500 K

    d.β, 200 K e.β, 300 K f.β, 500 K

    Fig.2 Final geometries forδandβphases containing a void defect after 1 ns NPT simulations at 200, 300 K and 500 K, respectively

    The behavior of MSD as a function of time is presented in Fig.3 and Fig.4. The data in the range of 600-1000 ps were used to plot the MSD curve. The behavior of MSD as a function of time can be used to discriminate between a solid and a liquid phase. The MSD of the solid system oscillates about a mean value. For a fluid, without underlying regular structure, the MSD gradually increases with time. So, from Fig.3, we can see that theδ-HMX crystal with a void defect presents part liquid-state feature at 200 K and 300 K, which is agreement with the geometries shown in Fig.2. When simulated temperature rises to 500 K, the system exhibits normal linear diffusive behavior, and this means that the lattice has completely collapsed into liquid state. Diffusion coefficient (D) is calculated by the well-known Einsteins relation, which obtains diffusion coefficients as a tangent of the linear region of the displacement function. The Einsteins relation is expressed as:

    (1)

    whereri(t) denotes the position vector at timet,ri(t0) the position vector at time origin, the angular bracket denotes averaging over all choices of time origin within a dynamics trajectory, andNis the total number of the diffusion particles. The diffusion coefficient at 500 K was calculated to be 1.2×10-10m2·s-1, significantly above the diffusion coefficient of baseline molecular of the order of 10-12m2·s-1. This value is observed for a molecule in the defect-free lattice and restricted to thermal fluctuations around its equilibrium site.

    Forβphase, the MSD curves are shown in Fig.4. At the temperature of 200 K, MSD oscillates about a mean value, showing an obvious solid-phase characteristic. At 300 K, we can still observe the oscillation phenomenon, but the MSD value increases slowly with runtime, indicating the lattice binding slightly reduces. When temperature rises up to 500 K, the case is very similar with theδ-HMX. Namely, the system finally changes into liquid state after 1 ns MD simulation. The difference is that normal linear region appears only after approximately 800 ps. This illustrates the rate of lattice collapsing is slower than that ofδphase. Besides, the diffusion coefficient is determined using the data of linear region 800-1000 ps, and the value is 1.5×10-10m2·s-1, almost the same as that ofδphase.

    a. 200 K, 300 K b. 500 K

    Fig.3 Calculated MSD forδ-HMX system containing a void defect at different temperatures

    a. 200 K, 300 K b. 500 K

    Fig.4 Calculated MSD forβ-HMX system containing a void defect at different temperatures

    3.3 Conformational Transition of Single Molecule with Time

    For collapsed liquid-state molecules, tracing the single molecule along with runtime, it can be observed the frequent transitions of molecule between different conformations. To illustrate this, we cut out a section of snapshots from MD simulation with a time period of 20 ps (see Fig.5). During 20 ps runtime, we have observed the obvious conformation transitions, namely from BB toα. Besides, we captured the transition state or intermediated structures undergone by this conversion, presenting the reversion of two downward NO2groups and adjustment of ring structure. When the system hasn′t completely melted into liquid state, for example the molecules collapsed into the void inδcrystal lattice under 200 K and 300 K, we can still observe the conformational transition, but the conversion rate is much slower due to the effects of the crystal field.

    The frequent transitions can be explained by the active barriers and the formation energy of the void defect. In section 3.2, we had already calculated the formation energies of voids inserted intoβandδcrystal lattices, and they are 7009.66 and 6833.35 kJ·mol-1, respectively. Combined with the active barriers listed in Table 1, we can predict that the formation energy is enough to induce conformational changes, rotations, and center-of-mass translations of the molecules in system, and even the melting of crystal lattice.

    Fig.5 Change of molecular conformations during 20 ps from the MD simulation for liquid-state system at 500 K

    3.4 Conformational Distribution in Final Simulation System

    At the simulation temperature of 500 K,no matterδorβphase, the system containing a void has collapsed into liquid state going through 1 ns simulation, and has the same density of 1.51 g·cm-3. So, we have thought the conformation distribution of final liquid system to be roughly similar for two crystal phases. And that is true through our conformational analysis. We find that some molecules are in stable conformations(about 40%), but a large number of molecules with transition state or intermediate structures(approximately 60%). In the different conformations, the ratio of BC is the largest(approximately 18%), next BB(approximately 12%), and the proportion ofαandβis relatively low(approximately 4% and 6%, respectively). The reason may be the combination of conversion barrier and the conformation stability. Such as, all the active barriers conversing to BC are relatively lower (see Table 1), and the energy of BC is the global minimum, which result in BC with highest ratio. Besides, a larger number of transition state and intermediate structures also make the collapsed system higher energy and more active. It can provide some evidence for the view that the collapse of voids seems to be one of the most important mechanisms that can lead to ignition in pressed explosives.

    Under lower temperatures of 200 K and 300 K, for the molecules collapsed into the void inδphase, we observe the conformation of most molecules is not change, still keeping theαconformation, owing to the effect of crystal field. Secondly, the transition state and intermediate are dominant in the rest of collapsed molecules. Besides, minority molecules are found with BC or BB structure.

    Forβphase, we have observed that the molecules on the surface of a void shift to the center of void at 200 K and 300 K from the Fig.2. The anisotropic interactions with others in their vicinities lead to variations in numbers and strengths of interactions as well as possible conformational changes. We find the conformational transition is mainly fromβto BC for the shifted molecules. It agrees with the lowest active barrier ofβ→BC among all transitions fromβto other conformations (10.33 kJ·mol-1, Table 1). Of course, the transition state and intermediates structures are also observed on the surface of void. These molecules have higher energy than that ofβ-HMX molecule, and will show more activity to initiate the degradation reaction of HMX. It had been demonstrated by DFT calculations coupled with transition state theory[12].

    4 Conclusions

    MD simulation for crystallineδ-HMX with a void defect was performed to study the dynamics characteristic of void evolvement and molecular conformation change, compared withβphase. The calculations for active energies of conformational transitions for the isolated molecules indicate that the higher barriers are 33.22 and 32.72 kJ·mol-1, and others are much lower than the two values. This means that all conversions can easily take place in energy.

    MD simulation shows that the crystalline system containing 10% vacancy concentration completely collapses into liquid phase after 1 ns NPT simulation at 500 K, no matterδorβphase. Compared to initial models, the lattice expands result in the system density down to 1.51 g·cm-3. For the liquid system, four conformations, namely,α,β, BB, and BC, have been observed. Among them, the proportion of BC is the largest. Besides, a large number of transition state and intermediate molecules are observed. Tracing the single molecule along with the runtime, frequent transitions between conformational states are found due to low energy barrier.

    At 200 K and 300 K, different evolvement trends are observed forδandβphases. Forδphase, the void has completely collapsed, and the whole crystal lattice has also lost the rigid periodicity. However, void and crystal lattice ofβphase have been hold on, and there only occur the shrink of void and some shift of molecules toward the center of void. Contrary to the volume shrink ofδphase, lattice expansion is observed forβphase. Besides, for the molecules collapsed into void or on the surface of void, except for the stable conformations, transition state or intermediate molecules have been found. Different evolvement behaviors of the voids contained inδandβphases may provide some important information to explain the sensitivity difference of two crystalline phases and the initial reaction mechanism of HMX.

    [1] Cobbledick R E, Small R W H. The crystal structure of theδ-form of 1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane (δ-HMX)[J].ActaCrystallographicaSectionB:StructuralCrystallographyandCrystalChemistry, 1974, 30(8): 1918-1922.

    [2] Goetz F, Brill T B, Ferraro J R. Pressure dependence of the raman and infraredspectra ofα-,β-,γ-, andδ- octahydro-1,3,5-7-tetranitro-1,3,5,7-tetrazocine[J].TheJournalofPhysicalChemistry, 1978, 82(17): 1912-1917.

    [3] Cady Howard H, Larson Allen C, Cromer Don T. The crystal structure ofα-HMX and a refinement of the structure ofβ-HMX[J].ActaCrystallographica, 1963,16(7): 617-623.

    [4] Cobbledick R E, Small R W H. The crystal structure of the δ-form of 1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane (δ-HMX)[J].ActaCrystallographicaSectionB:StructuralCrystallographyandCrystalChemistry, 1974,30(8): 1918-1922.

    [5] Brill T B, Reese C O. Analysis of intra- and intermolecular interactions relating to the thermophysical behavior of .alpha.-, .beta.-, and delta, -octahydro-1,3,5,7-tetranitro-1,3,5,7-tetraazocine[J].TheJournalofPhysicalChemistry, 1980, 84(11): 1376-1380.

    [6] Brill T B, Karpowicz R J. Solid phase transition kinetics. The role of intermolecular forces in the condensed-phase decomposition of octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine[J].TheJournalofPhysicalChemistry, 1982,86(21): 4260-4265.

    [7] Henson B F, Asay B W,Sander R K et al. Dynamic measurement of the HMXβ-δPhase transition by second harmonic generation[J].PhysicalReviewLetters, 1999,82(6): 1213-1216.

    [8] Urtiew P A, Forbes J W, Tarver C M, et al. Shock sensitivity of LX-04 containing delta phase HMX at elevated temperatures[C]∥Aip Conference Proceedings, 2004.

    [9] Lewis J P. Energetics of intermolecular HONO formation in condensed-phase octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX)[J].ChemicalPhysicsLetters, 2003,371(5): 588-593.

    [10] Chakraborty D, Muller R P, Dasgupta S, et al. Mechanism for unimolecular decomposition of HMX (1,3,5,7-tetranitro-1,3,5,7-tetrazocine), an ab initio study[J].TheJournalofPhysicalChemistryA, 2001, 105(8): 1302-1314.

    [11] Manaa M R, Fried L E, Melius C F, et al. Decomposition of HMX at extreme conditions: a molecular dynamics simulation[J].TheJournalofPhysicalChemistryA, 2002, 106(39): 9024-9029.

    [12] Sharia O, Kuklja M M. Rapid materials degradation induced by surfaces and voids: ab initio modeling of β-octatetramethylene tetranitramine[J].JournaloftheAmericanChemicalSociety, 2012,134(28): 11815-11820.

    [13] ZHOU Ting-ting, HUANG Feng-lei. Effects of defects on thermal decomposition of HMX via ReaxFF molecular dynamics simulations[J].TheJournalofPhysicalChemistryB, 2010,115(2): 278-287.

    [14] DUAN Xiao-hui, LI Wen-peng, PEI Chong-hua, et al. Molecular dynamics simulations of void defects in the energetic material HMX[J].JournalofMolecularModeling, 2013, 19(9): 3893-3899.

    [15] Sharia O, Tsyshevsky R, Kuklja M M. Surface-accelerated decomposition ofδ-HMX[J].TheJournalofPhysicalChemistryLetters,2013, 4(5): 730-734.

    [16] Choi C S, Boutin H P. A study of the crystal structure ofβ-cyclotetramethene tetranitramine by neutron diffraction[J].ActaCrystallographicaSectionB:StructuralCrystallographyandCrystalChemistry, 1970, 26(9): 1235-1240.

    [17] PENG Chun-yang, Schlegel H B.CombiningsynchronoustransitandQuasi-Newtonmethodstofindtransitionstates[J].IsraelJournalofChemistry, 1993, 33(4): 449-454.

    [18] PENG Chun-yang, Ayala P Y, Schlegel H B, et al. Using redundant internal coordinates to optimize equilibrium geometries and transition states[J].JournalofComputationalChemistry, 1996,(17): 49-56.

    [19] Gaussian 03: Revision B.03. Gaussian, Inc.: Pittsburgh, PA. 2003.

    [20] SUN Hai. An ab initio force-field optimized for condensed-phase applications overview with details on alkane and benzene compounds[J].TheJournalofPhysicalChemistryB, 1998,102(38): 7338-7364.

    [21] DUAN Xiao-hui, WEI Chun-xue, PEI Chong-hua, et al. A molecular dynamics simulation of solvent effects on the crystal morphology of HMX[J].JournalofHazardousMaterials, 2010,174(1): 75-180.

    [22] XIAO Ji-jun, WANG Wen-rui, CHEN Jun, et al. Study on the relations of sensitivity with energy properties for HMX and HMX-based PBXs by molecular dynamics simulation[J].PhysicaB:CondensedMatter, 2012, 407(17): 3504-3509.

    [23] XIAO Ji-jun, LI Song-yuan, CHEN Jun, et al. Molecular dynamics study on the correlation between structure and sensitivity for defective RDX crystals and their PBXs[J].JournalofMolecularModeling, 2013, 19(2): 803-809.

    [24] CUI Hong-ling, JI Guang-fu, CHEN Xiang-rong, et al. Phase transitions and mechanical properties of octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine in different crystal phases by molecular dynamics simulation[J].JournalofChemical&EngineeringData, 2010, 55(9): 3121-3129.

    [25] Cady H H, Smith L C. Studies on the polymorphs of HMX. Los Alamos Scientific Laboratory Report [R], LAMS-2652 TID-4500, Los Alamos National Laboratory, Los Alamos, NM, 1962.

    [26] Choi C S, Boutin H P. A study of the crystal structure ofβ-cyclotetramethylene tetranitramine by neutron diffraction[J].ActaCrystallographicaSectionB:StructuralCrystallographyandCrystalChemistry, 1970,26(9): 1235-1240.

    [27] Ewald P P. Evaluation of optical and electrostatic lattice potentials[J].AnnPhys, 1921,64: 253-287.

    [28] Rappe A K, Goddard W A, Charge equilibration for molecular dynamics simulations[J].TheJournalofPhysicalChemistry, 1991,95(8): 3358-3363.

    [29] Materials Studio 3.0, Accelrys Inc. San Diego, CA, 2004.

    [30] Smith G D, Bharadwaj R K. Quantum chemistry based force field for simulations of HMX[J].TheJournalofPhysicalChemistryB, 1999,103(18): 3570-3575.

    [31] Rosen J M, Dickinson C. Vapor pressures and heats of sublimation of some high-melting organic explosives[J].JournalofChemicalandEngineeringData, 1969,14(1): 120-124.

    [32] Taylor J W, Crookes R J. Vapour pressure and enthalpy of sublimation of 1,3,5,7-tetranitro-1,3,5,7-tetra-azacyclo-octane (HMX)[J].JournalofTheChemicalSociety,FaradayTransactions1:PhysicalChemistryinCondensedPhases, 1976, 72: 723-729.

    [33] Boyd S, Murray J S, Politzer P. Molecular dynamics characterization of void defects in crystalline (1,3,5-trinitro-1,3,5-triazacyclohexane)[J].TheJournalofChemicalPhysics, 2009,131(20): 204903.

    精品午夜福利在线看| 国产成人精品久久久久久| 亚洲国产精品国产精品| 成人漫画全彩无遮挡| 亚洲一区二区三区欧美精品 | 国产在线男女| 青春草国产在线视频| 免费观看的影片在线观看| 色哟哟·www| 日本与韩国留学比较| 白带黄色成豆腐渣| 国产色婷婷99| 欧美3d第一页| 插阴视频在线观看视频| 超碰av人人做人人爽久久| www.色视频.com| 国产女主播在线喷水免费视频网站| 一级毛片黄色毛片免费观看视频| 久久久久久久久久人人人人人人| 少妇熟女欧美另类| 色婷婷久久久亚洲欧美| 国国产精品蜜臀av免费| xxx大片免费视频| 黄片wwwwww| 国产av码专区亚洲av| av天堂中文字幕网| 波多野结衣巨乳人妻| 蜜桃亚洲精品一区二区三区| 能在线免费看毛片的网站| 毛片一级片免费看久久久久| 91aial.com中文字幕在线观看| 大香蕉久久网| 国产淫语在线视频| 夜夜爽夜夜爽视频| 午夜福利视频1000在线观看| 成人国产麻豆网| av网站免费在线观看视频| 男人和女人高潮做爰伦理| 黄色配什么色好看| 精品人妻熟女av久视频| 天美传媒精品一区二区| 亚洲激情五月婷婷啪啪| 少妇人妻精品综合一区二区| 亚洲综合色惰| 人妻一区二区av| 97在线人人人人妻| 亚洲色图综合在线观看| 我的老师免费观看完整版| 成人午夜精彩视频在线观看| 嫩草影院入口| 免费不卡的大黄色大毛片视频在线观看| 新久久久久国产一级毛片| 久久韩国三级中文字幕| 春色校园在线视频观看| 亚洲精品日韩av片在线观看| 日韩一区二区视频免费看| 久久久色成人| 黄色一级大片看看| 久久精品久久久久久久性| 亚洲国产欧美在线一区| 国产精品人妻久久久影院| 亚洲精品亚洲一区二区| 成人午夜精彩视频在线观看| 青春草国产在线视频| 91精品伊人久久大香线蕉| 免费电影在线观看免费观看| 2018国产大陆天天弄谢| 亚洲精品乱码久久久久久按摩| 黄色日韩在线| 国产伦在线观看视频一区| 精品一区二区三卡| 如何舔出高潮| 日日啪夜夜爽| 热re99久久精品国产66热6| 亚洲国产色片| 亚洲综合精品二区| 少妇猛男粗大的猛烈进出视频 | 蜜桃亚洲精品一区二区三区| 日本黄色片子视频| 一区二区三区免费毛片| 在线观看国产h片| 日韩制服骚丝袜av| 久久久久精品久久久久真实原创| 午夜福利高清视频| 在线看a的网站| 在线天堂最新版资源| 波野结衣二区三区在线| 天堂网av新在线| 亚洲人成网站高清观看| 亚洲精品一区蜜桃| 成人美女网站在线观看视频| 日日啪夜夜爽| 极品少妇高潮喷水抽搐| 亚洲人成网站在线播| 欧美丝袜亚洲另类| 偷拍熟女少妇极品色| 99热国产这里只有精品6| 久久99蜜桃精品久久| 久久久久精品久久久久真实原创| 精品久久久久久久久亚洲| 亚洲激情五月婷婷啪啪| 中文字幕av成人在线电影| 国产成人精品久久久久久| 在线a可以看的网站| 午夜福利高清视频| 欧美精品人与动牲交sv欧美| 寂寞人妻少妇视频99o| 国产精品久久久久久久电影| 一级毛片黄色毛片免费观看视频| 色综合色国产| 国产精品99久久久久久久久| 熟女av电影| 国产成人午夜福利电影在线观看| 少妇的逼水好多| 91精品一卡2卡3卡4卡| 老司机影院毛片| 精品视频人人做人人爽| 中文字幕人妻熟人妻熟丝袜美| 欧美3d第一页| 男人狂女人下面高潮的视频| 亚洲精品成人av观看孕妇| 熟女人妻精品中文字幕| 男女边吃奶边做爰视频| 精品午夜福利在线看| 亚洲综合色惰| 中国国产av一级| 麻豆国产97在线/欧美| 免费不卡的大黄色大毛片视频在线观看| 欧美日韩在线观看h| 亚洲真实伦在线观看| 搡女人真爽免费视频火全软件| 一区二区三区精品91| 只有这里有精品99| 国产精品福利在线免费观看| 日本一本二区三区精品| 别揉我奶头 嗯啊视频| 久久精品夜色国产| 精品人妻熟女av久视频| 久久久久精品性色| 如何舔出高潮| 熟妇人妻不卡中文字幕| 国产爽快片一区二区三区| 毛片女人毛片| 国产精品无大码| 欧美亚洲 丝袜 人妻 在线| 亚洲av欧美aⅴ国产| 免费av不卡在线播放| 欧美+日韩+精品| 人妻一区二区av| 人妻夜夜爽99麻豆av| 午夜福利高清视频| 汤姆久久久久久久影院中文字幕| 欧美日韩在线观看h| 久久精品久久久久久噜噜老黄| 777米奇影视久久| 日本猛色少妇xxxxx猛交久久| 国产精品爽爽va在线观看网站| 在线观看国产h片| 最近中文字幕高清免费大全6| av在线天堂中文字幕| 中国国产av一级| 国产老妇伦熟女老妇高清| 丝袜脚勾引网站| 国产中年淑女户外野战色| 三级国产精品片| 免费看av在线观看网站| 极品教师在线视频| 成人亚洲精品av一区二区| 噜噜噜噜噜久久久久久91| 免费大片18禁| 久久国内精品自在自线图片| 热99国产精品久久久久久7| 亚洲最大成人手机在线| 最近2019中文字幕mv第一页| 人妻系列 视频| 热re99久久精品国产66热6| 香蕉精品网在线| 别揉我奶头 嗯啊视频| 熟妇人妻不卡中文字幕| 亚洲精华国产精华液的使用体验| 国产毛片在线视频| 亚洲国产日韩一区二区| 美女脱内裤让男人舔精品视频| 肉色欧美久久久久久久蜜桃 | 大片电影免费在线观看免费| 久久久精品欧美日韩精品| 久久99热这里只频精品6学生| 69av精品久久久久久| 亚洲欧美日韩另类电影网站 | 日日摸夜夜添夜夜爱| 51国产日韩欧美| 18禁动态无遮挡网站| 中文在线观看免费www的网站| 亚洲人成网站高清观看| 一个人看的www免费观看视频| 日本免费在线观看一区| 午夜福利高清视频| 国产精品福利在线免费观看| 日韩精品有码人妻一区| 国产成人免费观看mmmm| 国产成人精品婷婷| 丰满少妇做爰视频| 乱系列少妇在线播放| 国产一级毛片在线| 久热久热在线精品观看| 狠狠精品人妻久久久久久综合| 国产精品人妻久久久久久| 国产精品一区二区三区四区免费观看| 插逼视频在线观看| 日韩人妻高清精品专区| 国产精品无大码| 国产高清国产精品国产三级 | 在线观看三级黄色| 免费大片18禁| 韩国av在线不卡| 久久久久久久国产电影| 国产精品麻豆人妻色哟哟久久| 成人一区二区视频在线观看| 色网站视频免费| 日韩av在线免费看完整版不卡| 日本午夜av视频| 亚洲自偷自拍三级| 一本久久精品| 狂野欧美白嫩少妇大欣赏| 91精品国产九色| 国产一区有黄有色的免费视频| 联通29元200g的流量卡| 成年版毛片免费区| 老司机影院毛片| 成人毛片a级毛片在线播放| 免费高清在线观看视频在线观看| 久久久国产一区二区| 久久精品综合一区二区三区| 一区二区三区四区激情视频| 2021少妇久久久久久久久久久| 韩国av在线不卡| 另类亚洲欧美激情| 91狼人影院| 国产男女超爽视频在线观看| 免费大片黄手机在线观看| 高清午夜精品一区二区三区| 制服丝袜香蕉在线| 精品少妇黑人巨大在线播放| 日韩,欧美,国产一区二区三区| 中文天堂在线官网| 啦啦啦在线观看免费高清www| 亚洲一级一片aⅴ在线观看| 国产精品久久久久久av不卡| 香蕉精品网在线| 久久这里有精品视频免费| 亚洲精品aⅴ在线观看| 国产欧美另类精品又又久久亚洲欧美| 亚洲av一区综合| 成人无遮挡网站| 国产爽快片一区二区三区| 亚洲精品国产av蜜桃| 女的被弄到高潮叫床怎么办| 日韩视频在线欧美| 特级一级黄色大片| 免费观看在线日韩| 在线观看一区二区三区| 菩萨蛮人人尽说江南好唐韦庄| av在线天堂中文字幕| 能在线免费看毛片的网站| 国产亚洲一区二区精品| 亚洲天堂国产精品一区在线| 小蜜桃在线观看免费完整版高清| 日韩欧美 国产精品| 亚洲欧美精品专区久久| 亚洲精品色激情综合| 欧美一区二区亚洲| 在现免费观看毛片| 超碰97精品在线观看| 久久久欧美国产精品| 欧美日韩综合久久久久久| 大陆偷拍与自拍| 国产乱人视频| 久久久久精品久久久久真实原创| 久久精品综合一区二区三区| 好男人视频免费观看在线| 99热全是精品| 中文字幕人妻熟人妻熟丝袜美| 欧美人与善性xxx| 成人毛片60女人毛片免费| 亚洲欧洲日产国产| 国产免费一区二区三区四区乱码| 免费黄色在线免费观看| 真实男女啪啪啪动态图| 亚洲最大成人中文| 男女下面进入的视频免费午夜| 2018国产大陆天天弄谢| 六月丁香七月| 男人狂女人下面高潮的视频| www.av在线官网国产| 久久99精品国语久久久| 大片电影免费在线观看免费| 免费大片18禁| 欧美 日韩 精品 国产| 成人二区视频| 日韩视频在线欧美| 91午夜精品亚洲一区二区三区| av免费在线看不卡| 最近中文字幕高清免费大全6| 国产精品熟女久久久久浪| 日韩一区二区视频免费看| 国产有黄有色有爽视频| 欧美bdsm另类| 噜噜噜噜噜久久久久久91| 成年免费大片在线观看| 赤兔流量卡办理| 国产午夜福利久久久久久| 亚洲欧美成人综合另类久久久| 好男人视频免费观看在线| 国产欧美日韩一区二区三区在线 | 99热全是精品| 免费黄网站久久成人精品| 乱码一卡2卡4卡精品| 尾随美女入室| 最近中文字幕2019免费版| 一级片'在线观看视频| 69人妻影院| 免费在线观看成人毛片| 中文字幕久久专区| 国产精品一区二区在线观看99| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 十八禁网站网址无遮挡 | 亚洲av电影在线观看一区二区三区 | 五月玫瑰六月丁香| 日韩免费高清中文字幕av| 国产色婷婷99| 欧美xxⅹ黑人| 少妇人妻久久综合中文| 国产成人免费无遮挡视频| 久久精品夜色国产| 日韩,欧美,国产一区二区三区| 亚洲欧美精品专区久久| 搡老乐熟女国产| 亚洲成人久久爱视频| 联通29元200g的流量卡| 只有这里有精品99| 99九九线精品视频在线观看视频| av女优亚洲男人天堂| 高清视频免费观看一区二区| 国产爱豆传媒在线观看| 亚洲av福利一区| 18禁在线播放成人免费| 国产精品蜜桃在线观看| 天堂网av新在线| kizo精华| 国产高潮美女av| 国产一区亚洲一区在线观看| 久久国内精品自在自线图片| 水蜜桃什么品种好| 精品人妻一区二区三区麻豆| 成人综合一区亚洲| 国产老妇女一区| 国产精品秋霞免费鲁丝片| 欧美激情在线99| 久热久热在线精品观看| 国产熟女欧美一区二区| 国产成人91sexporn| 亚洲国产精品成人综合色| 又爽又黄无遮挡网站| 免费看不卡的av| 3wmmmm亚洲av在线观看| 91在线精品国自产拍蜜月| 99热这里只有是精品在线观看| 亚洲精品一二三| 亚洲精品自拍成人| 一区二区三区四区激情视频| 免费看日本二区| 人人妻人人爽人人添夜夜欢视频 | 欧美日韩精品成人综合77777| 99久久精品热视频| 一级毛片我不卡| 国产亚洲av片在线观看秒播厂| 插逼视频在线观看| 国产精品久久久久久精品电影| 国产美女午夜福利| 99久久精品一区二区三区| 两个人的视频大全免费| 国产亚洲91精品色在线| 美女cb高潮喷水在线观看| 国产成人精品一,二区| 你懂的网址亚洲精品在线观看| 日韩强制内射视频| 天堂俺去俺来也www色官网| 中国美白少妇内射xxxbb| 街头女战士在线观看网站| 久久久成人免费电影| av在线播放精品| 丝袜美腿在线中文| 精品99又大又爽又粗少妇毛片| 国产白丝娇喘喷水9色精品| 草草在线视频免费看| 亚洲最大成人手机在线| 亚洲av成人精品一区久久| 制服丝袜香蕉在线| 亚洲自偷自拍三级| 亚洲国产色片| 黄色一级大片看看| 精品一区二区免费观看| videossex国产| 丝袜喷水一区| 免费看a级黄色片| 在线观看美女被高潮喷水网站| 国产69精品久久久久777片| 岛国毛片在线播放| 少妇人妻一区二区三区视频| 99久国产av精品国产电影| 在线观看一区二区三区激情| 久久久久久九九精品二区国产| 亚洲欧洲日产国产| 亚洲av在线观看美女高潮| 色网站视频免费| 各种免费的搞黄视频| 国产老妇伦熟女老妇高清| 大香蕉久久网| www.色视频.com| 久久国内精品自在自线图片| 国产淫语在线视频| 在线观看一区二区三区| 卡戴珊不雅视频在线播放| a级毛色黄片| 免费看光身美女| 午夜福利在线在线| 亚洲av电影在线观看一区二区三区 | 久久韩国三级中文字幕| 亚洲图色成人| 久久99精品国语久久久| 日韩成人av中文字幕在线观看| 成人亚洲精品av一区二区| 久久女婷五月综合色啪小说 | 在线观看免费高清a一片| 国产爽快片一区二区三区| 精品国产一区二区三区久久久樱花 | 欧美性感艳星| 成人亚洲精品av一区二区| 亚洲一区二区三区欧美精品 | 国产久久久一区二区三区| 国产毛片a区久久久久| 激情五月婷婷亚洲| 蜜桃亚洲精品一区二区三区| 精品国产露脸久久av麻豆| 舔av片在线| 国产v大片淫在线免费观看| 久久久久久久久久人人人人人人| xxx大片免费视频| 欧美xxxx性猛交bbbb| 亚洲电影在线观看av| 亚洲av男天堂| 精品一区二区三区视频在线| 亚洲性久久影院| 草草在线视频免费看| 男人狂女人下面高潮的视频| 老司机影院毛片| 高清av免费在线| 韩国av在线不卡| 两个人的视频大全免费| 听说在线观看完整版免费高清| 男女无遮挡免费网站观看| 婷婷色麻豆天堂久久| av在线播放精品| 国产精品久久久久久精品古装| 久久精品熟女亚洲av麻豆精品| 日产精品乱码卡一卡2卡三| 亚洲国产色片| 在线精品无人区一区二区三 | 国产极品天堂在线| 亚洲欧美日韩卡通动漫| 最近最新中文字幕大全电影3| 欧美日韩综合久久久久久| 国产成年人精品一区二区| 在线免费观看不下载黄p国产| 九草在线视频观看| 成人漫画全彩无遮挡| 免费黄网站久久成人精品| 两个人的视频大全免费| 青青草视频在线视频观看| 乱码一卡2卡4卡精品| 在线观看一区二区三区| 五月天丁香电影| 久久鲁丝午夜福利片| 亚洲av免费在线观看| 国产中年淑女户外野战色| 男的添女的下面高潮视频| 日本爱情动作片www.在线观看| 高清在线视频一区二区三区| av网站免费在线观看视频| 精品少妇久久久久久888优播| 欧美老熟妇乱子伦牲交| 久久久成人免费电影| 狂野欧美激情性bbbbbb| 精品酒店卫生间| 大码成人一级视频| 天堂俺去俺来也www色官网| 色网站视频免费| 亚洲精品国产色婷婷电影| 成年av动漫网址| 欧美日韩在线观看h| 高清欧美精品videossex| 亚洲无线观看免费| 亚洲av成人精品一二三区| 一本一本综合久久| 国产成人一区二区在线| 成人特级av手机在线观看| 国产久久久一区二区三区| 亚洲精品亚洲一区二区| 日本爱情动作片www.在线观看| 国产成人一区二区在线| 亚洲一级一片aⅴ在线观看| a级毛色黄片| 国产av码专区亚洲av| 美女xxoo啪啪120秒动态图| 久久精品国产a三级三级三级| av在线观看视频网站免费| 熟女电影av网| 欧美日韩视频高清一区二区三区二| 波多野结衣巨乳人妻| 国产国拍精品亚洲av在线观看| 美女视频免费永久观看网站| 亚洲欧美精品自产自拍| 国产黄片视频在线免费观看| 男人添女人高潮全过程视频| 亚洲自偷自拍三级| 国产69精品久久久久777片| 成人毛片60女人毛片免费| 男男h啪啪无遮挡| 一个人看的www免费观看视频| av一本久久久久| 国产久久久一区二区三区| 午夜福利在线在线| 久久人人爽人人片av| 精品国产露脸久久av麻豆| 午夜亚洲福利在线播放| 午夜视频国产福利| 91久久精品国产一区二区三区| 97在线人人人人妻| 一级毛片 在线播放| 免费少妇av软件| 麻豆国产97在线/欧美| 午夜视频国产福利| 久久影院123| 高清视频免费观看一区二区| 国产黄片视频在线免费观看| 伦理电影大哥的女人| 国模一区二区三区四区视频| 成人亚洲欧美一区二区av| 午夜福利视频1000在线观看| 午夜福利在线观看免费完整高清在| 色婷婷久久久亚洲欧美| 中文字幕av成人在线电影| 秋霞在线观看毛片| av播播在线观看一区| 18禁动态无遮挡网站| 日韩制服骚丝袜av| 国产高清三级在线| 久久人人爽人人爽人人片va| 欧美变态另类bdsm刘玥| 国产精品女同一区二区软件| 国产欧美亚洲国产| 国产成人a∨麻豆精品| 熟女电影av网| 国产淫片久久久久久久久| 中文字幕av成人在线电影| a级毛色黄片| 日本黄大片高清| 亚洲成人久久爱视频| av播播在线观看一区| 精品少妇久久久久久888优播| 可以在线观看毛片的网站| 亚洲国产精品专区欧美| 真实男女啪啪啪动态图| 亚洲国产av新网站| 中文字幕人妻熟人妻熟丝袜美| 日本色播在线视频| 好男人在线观看高清免费视频| 国产高潮美女av| 国产在视频线精品| 日本熟妇午夜| 精品久久久久久久末码| 日本熟妇午夜| 一级二级三级毛片免费看| 人妻制服诱惑在线中文字幕| 久久久精品欧美日韩精品| 精品视频人人做人人爽| 狂野欧美激情性xxxx在线观看| 在线观看一区二区三区激情| 菩萨蛮人人尽说江南好唐韦庄| tube8黄色片| 搡女人真爽免费视频火全软件| 国产精品久久久久久久久免| 国产免费一级a男人的天堂| 乱系列少妇在线播放| av.在线天堂| 在线观看免费高清a一片| 高清午夜精品一区二区三区| 午夜福利视频精品| 在线精品无人区一区二区三 | 国产真实伦视频高清在线观看| 男人和女人高潮做爰伦理| 十八禁网站网址无遮挡 | 男男h啪啪无遮挡| 永久网站在线| 最近中文字幕高清免费大全6| 免费黄频网站在线观看国产| 1000部很黄的大片| 青青草视频在线视频观看| 国产 精品1| 亚洲国产高清在线一区二区三| 国产精品久久久久久av不卡| 91狼人影院| 久久久久久久午夜电影| 男人添女人高潮全过程视频| 美女国产视频在线观看|