• <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.

    深爱激情五月婷婷| 深爱激情五月婷婷| 国产高清视频在线观看网站| 国内精品一区二区在线观看| 精品日产1卡2卡| 国产精品嫩草影院av在线观看| 国产91av在线免费观看| 午夜老司机福利剧场| 热99在线观看视频| 最近手机中文字幕大全| 国产伦在线观看视频一区| 女人被狂操c到高潮| 免费观看人在逋| 婷婷六月久久综合丁香| 免费观看在线日韩| 99热这里只有是精品50| 婷婷六月久久综合丁香| 三级国产精品欧美在线观看| 欧美日韩国产亚洲二区| 91久久精品电影网| 天堂网av新在线| 人人妻人人澡欧美一区二区| 国产av在哪里看| 日韩欧美精品v在线| 色av中文字幕| 欧美一级a爱片免费观看看| 亚洲精品国产av成人精品 | 久久精品久久久久久噜噜老黄 | 麻豆乱淫一区二区| 22中文网久久字幕| 欧美又色又爽又黄视频| 午夜激情欧美在线| 成年女人永久免费观看视频| 亚洲av免费在线观看| 丰满人妻一区二区三区视频av| 国产亚洲精品久久久久久毛片| 亚洲av中文字字幕乱码综合| 亚洲人成网站在线播| 最近的中文字幕免费完整| 插阴视频在线观看视频| 两个人视频免费观看高清| 精品乱码久久久久久99久播| 变态另类丝袜制服| 国产不卡一卡二| 国产成人精品久久久久久| 尾随美女入室| 淫妇啪啪啪对白视频| 国产精品嫩草影院av在线观看| 天美传媒精品一区二区| avwww免费| 99久久久亚洲精品蜜臀av| 国产黄a三级三级三级人| 一级黄片播放器| 欧美一区二区亚洲| 精品一区二区免费观看| 麻豆久久精品国产亚洲av| 男女那种视频在线观看| 免费不卡的大黄色大毛片视频在线观看 | 午夜视频国产福利| 我要看日韩黄色一级片| 啦啦啦观看免费观看视频高清| 成年女人永久免费观看视频| 亚洲丝袜综合中文字幕| 噜噜噜噜噜久久久久久91| 91午夜精品亚洲一区二区三区| 欧美潮喷喷水| 欧美另类亚洲清纯唯美| 1024手机看黄色片| 国产黄色小视频在线观看| 午夜福利在线观看吧| 国产真实伦视频高清在线观看| 老女人水多毛片| 最近中文字幕高清免费大全6| 两性午夜刺激爽爽歪歪视频在线观看| 国产一区二区亚洲精品在线观看| 日本 av在线| 久久热精品热| 变态另类成人亚洲欧美熟女| 伦精品一区二区三区| 日本熟妇午夜| 亚洲精品一卡2卡三卡4卡5卡| 91狼人影院| 午夜精品国产一区二区电影 | 欧美激情国产日韩精品一区| 久久精品国产亚洲av香蕉五月| 日韩欧美精品v在线| 麻豆乱淫一区二区| 1000部很黄的大片| 欧美成人免费av一区二区三区| 久久亚洲精品不卡| 亚洲精品久久国产高清桃花| 搡老妇女老女人老熟妇| 欧美bdsm另类| 成人三级黄色视频| 男人舔奶头视频| 欧美+亚洲+日韩+国产| 免费观看在线日韩| 在线a可以看的网站| 成人高潮视频无遮挡免费网站| 日韩精品中文字幕看吧| 亚洲在线自拍视频| 欧美中文日本在线观看视频| 久久久国产成人免费| 亚洲四区av| 极品教师在线视频| 一本精品99久久精品77| 日韩成人伦理影院| 老女人水多毛片| 国产在视频线在精品| 国产在视频线在精品| 国产精品野战在线观看| 午夜精品一区二区三区免费看| 天堂影院成人在线观看| 精品人妻一区二区三区麻豆 | 国产欧美日韩一区二区精品| 嫩草影院入口| 免费观看人在逋| 美女xxoo啪啪120秒动态图| 国产精品女同一区二区软件| 淫秽高清视频在线观看| 日本精品一区二区三区蜜桃| 久久久久久久久大av| 午夜影院日韩av| 亚洲av中文av极速乱| av中文乱码字幕在线| 久久精品人妻少妇| 国产三级在线视频| 亚洲精品在线观看二区| 国产 一区精品| 国内揄拍国产精品人妻在线| 国内少妇人妻偷人精品xxx网站| 日本精品一区二区三区蜜桃| 亚洲人与动物交配视频| 天堂av国产一区二区熟女人妻| 成熟少妇高潮喷水视频| 不卡一级毛片| 亚洲第一区二区三区不卡| 亚洲中文字幕日韩| 搡老岳熟女国产| 国产精品人妻久久久影院| av视频在线观看入口| 春色校园在线视频观看| 99riav亚洲国产免费| 人人妻人人看人人澡| 欧美最黄视频在线播放免费| aaaaa片日本免费| 亚洲人与动物交配视频| 国内精品宾馆在线| 日本一二三区视频观看| 亚洲aⅴ乱码一区二区在线播放| 久久久久久久久大av| 热99在线观看视频| 2021天堂中文幕一二区在线观| 国产精品一区二区性色av| 身体一侧抽搐| 一级黄色大片毛片| 欧美不卡视频在线免费观看| 国产精品99久久久久久久久| 国产亚洲精品综合一区在线观看| 美女被艹到高潮喷水动态| 亚洲中文日韩欧美视频| 亚洲性久久影院| 亚洲久久久久久中文字幕| 久久精品国产亚洲av天美| 中文字幕免费在线视频6| 国产av一区在线观看免费| 中文字幕人妻熟人妻熟丝袜美| 别揉我奶头 嗯啊视频| 国产精品永久免费网站| 国产 一区 欧美 日韩| 国产高清三级在线| 亚洲av中文av极速乱| 亚洲精品成人久久久久久| 国产一区二区亚洲精品在线观看| 99热精品在线国产| 欧美高清性xxxxhd video| 直男gayav资源| 午夜久久久久精精品| 国产伦一二天堂av在线观看| 国产成人精品久久久久久| 精品国产三级普通话版| 国产男人的电影天堂91| 色播亚洲综合网| 国产精品国产高清国产av| 欧美中文日本在线观看视频| 亚洲精品色激情综合| 一级a爱片免费观看的视频| 日产精品乱码卡一卡2卡三| 日韩精品有码人妻一区| 亚洲国产欧洲综合997久久,| 噜噜噜噜噜久久久久久91| 国产精品久久久久久av不卡| 两个人视频免费观看高清| 亚洲精品久久国产高清桃花| 禁无遮挡网站| 禁无遮挡网站| 精品一区二区三区视频在线| 欧美成人精品欧美一级黄| 亚洲最大成人av| 亚洲五月天丁香| 在线看三级毛片| 日本一二三区视频观看| 精品人妻一区二区三区麻豆 | 久久久久精品国产欧美久久久| 亚洲天堂国产精品一区在线| 成人美女网站在线观看视频| 亚洲av.av天堂| 国产成年人精品一区二区| 有码 亚洲区| 我要搜黄色片| 1024手机看黄色片| 少妇熟女欧美另类| 美女高潮的动态| 91精品国产九色| 欧美高清成人免费视频www| 精品人妻一区二区三区麻豆 | 成人亚洲精品av一区二区| 午夜爱爱视频在线播放| 露出奶头的视频| 国产精品一区www在线观看| 亚洲精品影视一区二区三区av| 国产高清有码在线观看视频| 深爱激情五月婷婷| 亚洲精品影视一区二区三区av| 日韩av在线大香蕉| 麻豆成人午夜福利视频| 男插女下体视频免费在线播放| 99久久精品国产国产毛片| 久久久久国内视频| 欧美不卡视频在线免费观看| 九色成人免费人妻av| 久久久久精品国产欧美久久久| 女的被弄到高潮叫床怎么办| 亚洲成人精品中文字幕电影| 亚洲真实伦在线观看| 国产精品久久久久久亚洲av鲁大| 在线观看av片永久免费下载| 精品人妻偷拍中文字幕| 男女啪啪激烈高潮av片| 一进一出抽搐gif免费好疼| 12—13女人毛片做爰片一| 久久久久精品国产欧美久久久| 久久久久久久久中文| 国产一区二区在线观看日韩| 午夜福利18| 亚洲成人中文字幕在线播放| 综合色丁香网| 国产美女午夜福利| 国产在视频线在精品| 精品久久久久久久久久免费视频| 男人和女人高潮做爰伦理| 99在线人妻在线中文字幕| 亚洲成人久久爱视频| 国产精品一二三区在线看| 又爽又黄a免费视频| 免费一级毛片在线播放高清视频| 熟妇人妻久久中文字幕3abv| 久久久久久国产a免费观看| 国产欧美日韩一区二区精品| 亚洲一区高清亚洲精品| 最近在线观看免费完整版| 97人妻精品一区二区三区麻豆| 免费av毛片视频| 亚洲内射少妇av| 超碰av人人做人人爽久久| 秋霞在线观看毛片| 赤兔流量卡办理| 久久热精品热| 欧美性猛交黑人性爽| 欧美丝袜亚洲另类| 一卡2卡三卡四卡精品乱码亚洲| 别揉我奶头 嗯啊视频| 99精品在免费线老司机午夜| 日本熟妇午夜| 麻豆久久精品国产亚洲av| 午夜福利18| 丰满的人妻完整版| 国产精品三级大全| 卡戴珊不雅视频在线播放| 老司机影院成人| 男人狂女人下面高潮的视频| 国产探花在线观看一区二区| 少妇高潮的动态图| 少妇高潮的动态图| 婷婷精品国产亚洲av| 麻豆久久精品国产亚洲av| 午夜福利在线在线| 天堂动漫精品| 免费大片18禁| 中国美白少妇内射xxxbb| 伦精品一区二区三区| 日本免费一区二区三区高清不卡| 亚洲自拍偷在线| 精品一区二区免费观看| 美女cb高潮喷水在线观看| 嫩草影院新地址| or卡值多少钱| 亚洲精品久久国产高清桃花| 秋霞在线观看毛片| 搡老熟女国产l中国老女人| 日本免费一区二区三区高清不卡| 黄片wwwwww| 亚洲av成人av| 日韩成人av中文字幕在线观看 | 非洲黑人性xxxx精品又粗又长| 国产aⅴ精品一区二区三区波| 国产黄片美女视频| 91av网一区二区| 一个人免费在线观看电影| 亚洲成a人片在线一区二区| 一级a爱片免费观看的视频| 国产av不卡久久| 男人的好看免费观看在线视频| 国产又黄又爽又无遮挡在线| 成年版毛片免费区| 日韩精品有码人妻一区| 国产成人a∨麻豆精品| 亚洲国产精品久久男人天堂| 人妻制服诱惑在线中文字幕| 国产精品无大码| 免费观看人在逋| 人妻久久中文字幕网| 看免费成人av毛片| 亚洲精品日韩av片在线观看| 一区二区三区四区激情视频 | 最近2019中文字幕mv第一页| 美女被艹到高潮喷水动态| 日本熟妇午夜| 亚洲成人中文字幕在线播放| 免费不卡的大黄色大毛片视频在线观看 | 国产真实乱freesex| 国产aⅴ精品一区二区三区波| 在线观看免费视频日本深夜| 菩萨蛮人人尽说江南好唐韦庄 | 如何舔出高潮| 日韩欧美在线乱码| 国产亚洲欧美98| 在线播放国产精品三级| 特级一级黄色大片| 亚洲精品一区av在线观看| 美女xxoo啪啪120秒动态图| 免费av毛片视频| 成人午夜高清在线视频| av在线老鸭窝| 看十八女毛片水多多多| 神马国产精品三级电影在线观看| 日韩欧美免费精品| 女的被弄到高潮叫床怎么办| 国产黄片美女视频| 日韩三级伦理在线观看| 欧美xxxx性猛交bbbb| 欧美人与善性xxx| 啦啦啦啦在线视频资源| h日本视频在线播放| 国产精品,欧美在线| videossex国产| 亚洲国产精品国产精品| 亚洲成人久久性| 亚洲电影在线观看av| 欧美bdsm另类| 精品久久久久久久末码| 国产精品电影一区二区三区| 久久亚洲精品不卡| videossex国产| 精品人妻一区二区三区麻豆 | 欧美不卡视频在线免费观看| 精品一区二区免费观看| 欧美在线一区亚洲| 国产精品久久电影中文字幕| 亚洲图色成人| 亚洲精品国产av成人精品 | 99热网站在线观看| 禁无遮挡网站| 在线观看一区二区三区| 欧美成人a在线观看| 欧美国产日韩亚洲一区| 小说图片视频综合网站| 日韩欧美三级三区| 亚洲欧美日韩高清专用| 精品久久久噜噜| 久久欧美精品欧美久久欧美| 国产日本99.免费观看| 又黄又爽又免费观看的视频| 三级毛片av免费| 国产又黄又爽又无遮挡在线| 国产精品无大码| av在线天堂中文字幕| 美女 人体艺术 gogo| 特级一级黄色大片| 国产一区二区在线观看日韩| 秋霞在线观看毛片| 久久久久久久久久成人| av在线播放精品| 午夜福利成人在线免费观看| 亚洲国产精品合色在线| 在线观看一区二区三区| 男人的好看免费观看在线视频| 国产欧美日韩精品一区二区| 白带黄色成豆腐渣| 亚洲真实伦在线观看| 成人国产麻豆网| 嫩草影院新地址| 精品人妻熟女av久视频| 精品日产1卡2卡| 午夜视频国产福利| 一区二区三区四区激情视频 | 中文在线观看免费www的网站| 身体一侧抽搐| 亚洲内射少妇av| 亚洲精品亚洲一区二区| 精华霜和精华液先用哪个| 性色avwww在线观看| 全区人妻精品视频| 一级毛片电影观看 | 97超碰精品成人国产| 波多野结衣巨乳人妻| 国产av在哪里看| 69人妻影院| 欧洲精品卡2卡3卡4卡5卡区| 欧美不卡视频在线免费观看| 五月玫瑰六月丁香| 一级黄片播放器| 老熟妇乱子伦视频在线观看| 五月玫瑰六月丁香| 一级av片app| 熟女电影av网| 久久精品综合一区二区三区| 色av中文字幕| 中文字幕av成人在线电影| 久久国内精品自在自线图片| 国产伦一二天堂av在线观看| 亚洲av成人av| 人妻少妇偷人精品九色| 桃色一区二区三区在线观看| 国产在线男女| 亚洲精品日韩av片在线观看| 99riav亚洲国产免费| 真人做人爱边吃奶动态| 国产男靠女视频免费网站| 午夜福利视频1000在线观看| 成人欧美大片| 亚洲成人精品中文字幕电影| 老司机影院成人| 欧美激情久久久久久爽电影| 啦啦啦观看免费观看视频高清| 91午夜精品亚洲一区二区三区| av中文乱码字幕在线| 一夜夜www| 国产精品一区二区免费欧美| 日韩精品青青久久久久久| 久久6这里有精品| 啦啦啦韩国在线观看视频| 亚洲精品在线观看二区| 男人的好看免费观看在线视频| 午夜精品一区二区三区免费看| 日本色播在线视频| 亚洲四区av| 国产精品国产高清国产av| 亚洲国产欧洲综合997久久,| 久久精品夜色国产| 亚洲专区国产一区二区| 成人特级av手机在线观看| 亚洲丝袜综合中文字幕| 国产欧美日韩精品一区二区| 永久网站在线| 美女xxoo啪啪120秒动态图| 亚洲在线自拍视频| 俄罗斯特黄特色一大片| 国产精品av视频在线免费观看| 在现免费观看毛片| 啦啦啦观看免费观看视频高清| 久久久久国产精品人妻aⅴ院| 人人妻人人澡欧美一区二区| 亚洲av一区综合| 久久人妻av系列| www日本黄色视频网| 免费观看的影片在线观看| 久久久久久久久久成人| 国产男人的电影天堂91| 最近最新中文字幕大全电影3| 此物有八面人人有两片| 日韩高清综合在线| 少妇的逼水好多| 我的女老师完整版在线观看| 精品久久久噜噜| av国产免费在线观看| 久久精品夜色国产| 欧美日韩综合久久久久久| 最近视频中文字幕2019在线8| 国产大屁股一区二区在线视频| 欧美日韩在线观看h| 美女cb高潮喷水在线观看| 一个人看视频在线观看www免费| 国产成人精品久久久久久| 男女之事视频高清在线观看| 青春草视频在线免费观看| 欧美bdsm另类| 国产精品一区www在线观看| 精品少妇黑人巨大在线播放 | 国产精品女同一区二区软件| 美女被艹到高潮喷水动态| 91麻豆精品激情在线观看国产| 欧美另类亚洲清纯唯美| 美女免费视频网站| 99久国产av精品| av天堂在线播放| 午夜精品在线福利| 一卡2卡三卡四卡精品乱码亚洲| 寂寞人妻少妇视频99o| 久久久久九九精品影院| 免费av不卡在线播放| 性插视频无遮挡在线免费观看| 特级一级黄色大片| 18禁在线播放成人免费| 国产色爽女视频免费观看| 热99在线观看视频| 色视频www国产| 国产免费一级a男人的天堂| 中国美女看黄片| 久久久久久久久久成人| 国产国拍精品亚洲av在线观看| 久久欧美精品欧美久久欧美| videossex国产| 精品少妇黑人巨大在线播放 | 国产一区二区三区在线臀色熟女| 欧美中文日本在线观看视频| 我的老师免费观看完整版| 久久天躁狠狠躁夜夜2o2o| 搞女人的毛片| 舔av片在线| 国产激情偷乱视频一区二区| 少妇的逼好多水| 午夜精品一区二区三区免费看| 麻豆国产97在线/欧美| 乱码一卡2卡4卡精品| 亚洲精品一区av在线观看| 日韩 亚洲 欧美在线| 精品福利观看| 18禁在线无遮挡免费观看视频 | 婷婷精品国产亚洲av在线| 久久精品国产亚洲网站| 亚洲精品在线观看二区| 日本撒尿小便嘘嘘汇集6| 日本五十路高清| 蜜桃亚洲精品一区二区三区| 一个人免费在线观看电影| 国产精品不卡视频一区二区| 超碰av人人做人人爽久久| 久久久久国产网址| 少妇裸体淫交视频免费看高清| 免费av观看视频| 欧美性感艳星| 日本熟妇午夜| 亚洲精品影视一区二区三区av| 国产亚洲精品久久久久久毛片| 午夜福利在线观看吧| 国产精品久久久久久久电影| 成人特级av手机在线观看| 51国产日韩欧美| 在线观看66精品国产| 日本爱情动作片www.在线观看 | 热99在线观看视频| 色尼玛亚洲综合影院| 日韩精品有码人妻一区| 亚洲成a人片在线一区二区| 久久欧美精品欧美久久欧美| 亚洲性久久影院| 成人三级黄色视频| 99热这里只有是精品在线观看| 日本黄大片高清| 国产亚洲av嫩草精品影院| 小说图片视频综合网站| 99久久九九国产精品国产免费| 免费大片18禁| 国产成人一区二区在线| or卡值多少钱| 久久久久久久久中文| 久久久久久久午夜电影| 成人一区二区视频在线观看| 国产一区亚洲一区在线观看| 最近最新中文字幕大全电影3| 色综合亚洲欧美另类图片| 日本三级黄在线观看| 男女下面进入的视频免费午夜| 国产精品久久久久久久久免| 免费看a级黄色片| 国产日本99.免费观看| 国产欧美日韩精品亚洲av| 观看免费一级毛片| 99热网站在线观看| 别揉我奶头 嗯啊视频| 综合色丁香网| 婷婷亚洲欧美| videossex国产| 国产精品一区二区性色av| 一本一本综合久久| 又爽又黄无遮挡网站| 此物有八面人人有两片| 美女内射精品一级片tv| 变态另类成人亚洲欧美熟女| 亚洲精品国产av成人精品 | 校园人妻丝袜中文字幕| 搡女人真爽免费视频火全软件 | 高清午夜精品一区二区三区 | 1024手机看黄色片| 大型黄色视频在线免费观看| 日韩 亚洲 欧美在线| 神马国产精品三级电影在线观看| 97超级碰碰碰精品色视频在线观看| 亚洲自拍偷在线| 晚上一个人看的免费电影| 婷婷色综合大香蕉| 一本精品99久久精品77| 国国产精品蜜臀av免费| 特级一级黄色大片|