J. Phys. Chem. B 2008, 112, 11005–11013
11005
Ab Initio Molecular Dynamics Study on the Initial Chemical Events in Nitramines: Thermal Decomposition of CL-20 Olexandr Isayev,† Leonid Gorb,†,‡ Mo Qasim,‡ and Jerzy Leszczynski*,†,‡ Computational Center of Molecular Structure and Interactions, Jackson State UniVersity, Jackson, Mississippi 39217, and U.S. Army Engineer Research and DeVelopment Center, Vicksburg, Mississippi 39180 ReceiVed: May 29, 2008
CL-20 (2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane or HNIW) is a high-energy nitramine explosive. To improve atomistic understanding of the thermal decomposition of CL-20 gas and solid phases, we performed a series of ab initio molecular dynamics simulations. We found that during unimolecular decomposition, unlike other nitramines (e.g., RDX, HMX), CL-20 has only one distinct initial reaction channelshomolysis of the N-NO2 bond. We did not observe any HONO elimination reaction during unimolecular decomposition, whereas the ring-breaking reaction was followed by NO2 fission. Therefore, in spite of limited sampling, that provides a mostly qualitatiVe picture, we proposed here a scheme of unimolecular decomposition of CL-20. The averaged product population over all trajectories was estimated at four HCN, two to four NO2, two to four NO, one CO, and one OH molecule per one CL-20 molecule. Our simulations provide a detailed description of the chemical processes in the initial stages of thermal decomposition of condensed CL-20, allowing elucidation of key features of such processes as composition of primary reaction products, reaction timing, and Arrhenius behavior of the system. The primary reactions leading to NO2, NO, N2O, and N2 occur at very early stages. We also estimated potential activation barriers for the formation of NO2, which essentially determines overall decomposition kinetics and effective rate constants for NO2 and N2. The calculated solid-phase decomposition pathways correlate with available condensed-phase experimental data. 1. Introduction CL-20 (2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane or HNIW; see Figure 1 for structure) is a nitramine explosive belonging to the class of substances possessing R2N-NO2 functional groups. The best known members of this group are RDX (hexogen or cyclotrimethylenetrinitramine), HMX (octogen or cyclotetramethylenetetranitramine), and the caged nitramine, CL-20. Past studies of thermal decomposition mechanisms of explosives and propellants have been designed to determine the detailed chemical pathways by which these large organic molecules are transformed into small reaction products (N2, NO, N2O, HCN, CO, CO2, H2O, etc.). This information is crucial to predicting long-term stability of explosives for purposes of storage, sensitivity to various stimuli such as heat and mechanical impact, as well as risk assessment associated with their usage. Most common explosives and nitro propellants have been identified as toxic to a variety of organisms in sediments and aquatic environments.1-4 Furthermore, some degradation products may be more toxic than the parent compounds.5-7 Unfortunately, without detailed mechanistic information, the behavior of explosives remains disconcertingly elusive. In addition, it is hoped that the development of better chemical models and detailed kinetics schemes will aid in the design of new energetic materials that are safer to manufacture, store, and handle. Depending on conditions, two general pathways of thermal decomposition of nitramine energetic compounds, such as CL* Corresponding author. E-mail:
[email protected]. † Jackson State University. ‡ U.S. Army Engineer Research and Development Center.
Figure 1. Structure of the CL-20 molecule.
20, can be distinguished. The unimolecular pathway could involve an initiated gas-phase thermolysis, photolysis, or ionization in molecular beams,8-11 whereas thermal decomposition in the condensed phase is usually associated with complex multimolecular processes.12 The identification of both pathways represents a challenging experimental problem, because the rapid release of heat during decomposition and the fast reaction rate make it difficult to isolate obtained species and identify individual steps in the mechanism, especially during the early stages of the reaction. Therefore, current experimental studies are mainly limited to analysis of postreaction residues and thermochemical studies.2,13-17 Nevertheless, they provide valuable but indirect insight on global decomposition kinetics, activation energy barriers, and final gaseous decomposition products. For example, it was established that RDX primarily decomposes in four different pathways: concerted breaking of ring C-N bonds,8 N-NO2 homolysis,13 successive HONO elimination,18 and oxygen migration from NO2 to the neighboring C atom.13 In contrast, no oxygen migration from NO2 has been described for HMX.14 Experimental data on the initial pathways
10.1021/jp804765m CCC: $40.75 2008 American Chemical Society Published on Web 08/08/2008
11006 J. Phys. Chem. B, Vol. 112, No. 35, 2008 of CL-20 thermal decomposition is less than for RDX and HMX in that only one channel of decomposition was identified (N-NO2 homolysis).13 Unlike studies of RDX and most other nitramines, those devoted to CL-20 are limited to condensedphase decomposition.13-17 Moreover, detailed kinetics for most primary and secondary reactions as well as formation time of key intermediates are not known. There are reported activation barriers (Ea) that vary from 120 to 200 kJ/mol, depending on such conditions as heating rates, stage of decomposition, and sample size. Also, the experimental pre-exponential factor (A) deviates as much as 5 orders of magnitude. Theoretical (both quantum mechanical [QM] and molecular mechanics [MM] based) atomic-scale simulations, which are both capable of providing the subpicosecond time resolution needed to study such fast events, provide alternatives to experimental methods. In comparison to experimental studies, other advantages of computational methods include safety and, also, significant savings in time and cost. However, current theoretical methods are also known to have limitations. Specifically, the MM-based simulations are not able to properly describe the bond-breaking processes (despite recent attempts to overcome this problem via of ReaxFF force field development19-24). The main challenges of canonical ab initio quantum chemistry arise from the rapid branching of reaction channels. This number grows factorially with system size, effectively preventing individual calculations. Nevertheless, most progress in theoretical nitro explosives analysis has been achieved using these two approximations. For example, comprehensive studies devoted to theoretical analysis of different aspects of RDX and HMX thermolysis conducted by the Goddard group25,26 allowed construction of a unified scheme for unimolecular decomposition of RDX and HMX, including various secondary channels calculated at the density functional theory (DFT) level. Rate constant values and temperature dependence of RDX and HMX decomposition have been also evaluated.27 Another important study contributed by Truong and co-workers,28,29 includes temperature and pressure dependence of rate constants for HMX N-NO2 homolysis and HONO elimination (tunneling effect was also included in considered HONO elimination). Just recently Naoto et al. studied decomposition and chemisorption of RDX on Al(111) surface.30 State-of-the-art work in computational modeling of thermolysis of explosives was exemplified in the development of the ReaxFF force field by Goddard and co-workers,23,27 making the study of thermal decomposition and shock-induced chemistry of crystalline RDX possible. Further extension of this reactive force field has also allowed analysis of thermal initiation of crystalline triacetoneperoxide.24 Although first-principle simulations of materials under extreme conditions have been conducted for simple systems, such as nitromethane31 and NTO32 (5-nitro-2,4-dihydro-3H-1,2,4triazol-3-one), most large-scale simulations of species of practical interest have to rely on classical interatomic potentials33 or tight-binding approximations, particularly SCC-DFTB.34 The growing demand for ab initio molecular dynamics (AIMD) simulations has driven the development of massively parallel high-performance computer systems.35,36 Common availability of the teraflop-level computer environment and recent petaflop computational platform initiative from the National Science Foundation and other governmental agencies provides the ultimate framework for further progress in simulation models.37 Altogether, this, coupled to the advances in AIMD methodology,38,39 and, particularly, the Car-Parrinello frame-
Isayev et al. work,40 makes it possible to study such complex chemical phenomena fully ab initio. The purposes of this study are to investigate the initial chain of chemical events during CL-20 decomposition of isolated and condensed phases and, additionally, to analyze the nature, energetic, and kinetic patterns of the primary and secondary products. We utilized AIMD to address these issues. The elementary reactions involved in the complicated decomposition process were obtained directly from the MD trajectories. First, a set of trajectories was collected by simulating an isolated CL20 molecule at various temperatures. In addition to unimolecular decomposition, significant effort was also spent on analysis of condensed-phase reactions. Consequently, the unimolecular decomposition process was analyzed only qualitatiVely. In the second part of the research, we produced an overview of atomistic mechanisms of decomposition of crystalline CL-20 and estimated the rate constants of the primary and secondary reactions under constant high-temperature conditions. 2. Methods and Models 2.1. Molecular Dynamics. All calculations were performed within the Kohn-Sham DFT method in its plane wave/ pseudopotential formulation, using the CPMD code.41 The BLYP exchange-correlation functional was chosen, and core electrons were taken into account using the Troullier-Martins norm-conserving pseudopotentials.42 For the set of applied pseudopotentials, we utilized a plane wave cutoff of 85 Ry for molecular dynamics and 120 Ry for geometry and cell optimizations. These values represent an optimum tradeoff between accuracy and computational feasibility of simulations and were sufficient to obtain convergence in cell volume, MD forces, bond length values, and other molecular properties. Constant temperature AIMD simulations were performed using the Car-Parrinello scheme40 along with Nose´-Hoover chain thermostats43 for nuclei. Each system was carefully thermostatted to 300 K for at least 5 ps and heated to target temperature for an additional 0.5 ps, keeping the overall structure rigid. Time steps of 4 and 5 au (0.097 and 0.121 fs) were used for the integration of the equations of motion during production and equilibration runs, respectively. The assumed mass for the orbitals was 700 au, and the hydrogen mass was substituted by the deuterium mass. All initial structures were taken from the Cambridge Structural Database (CSD). 2.2. Isolated System. For calculation of the gas-phase molecule, an 18 Å cubic supercell was used with a cell length chosen to be large enough to prevent interactions between periodic images. The effect of the unit cell length was checked by expanding the box size to approximately 20 Å. Only a minor variation in the total energy was found. The geometry of the system was minimized and subsequently equilibrated at the NVT ensemble. After equilibration, five independent production runs with target temperatures of 1000, 1500, 2000, 2500, and 3000 K were initiated. To improve observed statistics, each simulation was repeated at least one more time, starting from a different configuration. The length of simulations varied from 5 to 10 ps, and actual temperature fluctuations were within 50-150 K from those desired. 2.3. Crystal Phase. ε-CL-20 belongs to the monoclinic space group P21/n with experimental lattice vectors a ) 13.696, b ) 12.254, c ) 8.833 Å and angle β ) 111.18°.44 One unit cell consists of four molecules (144 atoms). Both geometry and cell volume were optimized using standard approaches within the CPMD package. The space group of the cell was restricted during optimization, allowing only isotropic changes. Optimized
Thermal Decomposition of CL-20
J. Phys. Chem. B, Vol. 112, No. 35, 2008 11007
Figure 2. Proposed scheme of unimolecular decomposition of CL-20, based on series AIMD simulations at 1500-3000 K. Initial intermediates are summarized by the shadowed rectangle. Please see Figure 3 for structures and more explanations in text.
lattice parameters were a ) 13.741, b ) 12.607, c ) 9.095 Å and angle β ) 111.36°. Calculated values averaged 1% larger than the experimental structural parameters. The largest deviation was found for lattice vector c at 0.26 Å (∼3%). The obtained lattice parameters resulted in a density of ε-CL-20 at 1.88 g/cm3, 6% smaller than the experimental density. These numbers are in perfect agreement with earlier works employing PBE, rPBE, or PW91 functionals, utilizing ultrasoft pseudopotentials.45,46 After 5 ps of equilibration, four independent production runs with target temperatures of 1500, 2000, 2500, and 3000 K were initiated. The actual temperatures were 1488, 1969, 2547, and 3124, respectively. The length of simulations varied from 8 to 10 ps. In order to obtain better statistical sampling and test system size dependence, we performed simulations with larger supercells (1 × 1 × 2 unit cells, eight CL-20 molecules, 288 atoms). The system was constructed by merging two equilibrated unit cells along the c-axis (smallest axis). After a 1 ps quick reequilibration, production runs were initiated at 2000 and 3000 K. The lengths of simulations were 5.1 (aborted due to hardware failure) and 15.4 ps, respectively; the actual temperatures were 1947 and 2882 K, respectively. 2.4. Analysis and Postprocessing. In order to follow the chemical processes throughout the MD trajectories, we had to define target molecules. This is commonly determined via configuration space: when two atoms are closer than a given cutoff distance, they belong to the same molecule. This cutoff distance is driven by the difference covalent (∼1.5 Å) and van der Waals bond distances (∼3 Å). However, under such extreme conditions, this approach has two drawbacks: (a) two atoms may be close in configurational space for times shorter than a vibrational period23 and (b) some transition states can be incorrectly attributed to the molecular species. Therefore, these problems affected time averaged concentrations only insignificantly. These postprocessing procedures, accounting for periodic boundary conditions, were programmed using the MATLAB suite.
3. Results and Discussion 3.1. Mechanism of Unimolecular Decomposition of CL20. On the basis of a series of CPMD simulations of isolated CL-20 in the 1500-3000 K temperature range, we proposed pathways of its unimolecular decomposition, depicted in Figure 2. Unlike other nitramines (e.g., RDX/HMX25,26,47), we found only one distinct initial decomposition channelshomolysis of the N-NO2 bond (NO2 fission). In particular, we did not observe any HONO elimination reaction under these circumstances, whereas NO2 fission was followed by ring-breaking reaction. Our findings completely agree with experimental studies showing dominance of NO2 fission in the early stages of thermal decomposition. This reaction also determines the overall activation barrier.43 Homolysis of the N-NO2 bond occurs on the 0.1-0.5 ps scale, depending on temperature. As expected, the higher the temperature, the faster the reaction proceeds. We found no clear temperature dependence of the number of broken N-NO2 bonds. Seven out of 10 trajectories (including the most and least extreme) resulted in fission of all six nitro groups. This is in perfect agreement with the mass spectral pattern, where peaks corresponding to elimination of three to six NO2 groups were observed.16,48 We did not observe any distinct differences between axial and equatorial nitro groups, which is probably due to the fact that, under such extreme conditions, the activation barrier difference of 10-20 kJ/mol49 becomes negligible. One of the few possible intermediates is produced by the departure of the NO2 radicals from CL-20. Structures of intermediates are depicted in Figure 3. It is quite probable that, due to the limited statistical sampling, this is not an exhaustive list of all possible geometrical variations. However, 7 out of 10 trajectories resulted in 12-atomic intermediates (Figure 3a). This species oscillates around planar geometry and is stabilized by adoption of an aromatic structure with typical C-N bond lengths, ranging from 1.2 to 1.45 Å, and having a lifetime of approximately 0.5 ps at 3000 K and quickly increasing to 9 ps at 1500 K.
11008 J. Phys. Chem. B, Vol. 112, No. 35, 2008
Isayev et al.
Figure 3. Typical snapshots of CL-20 intermediate after homolysis of N-NO2 bonds. All neighbor molecules are omitted for clarity. (a) Most abundant 12-atomic cyclic intermediate, observed in the 1000-3000 K temperature range; (b) tailed intermediate, observed in two simulations; (c) linear intermediate, found in one simulation, T ) 1500 K.
Figure 4. Schematic representation of several observed secondary reactions during unimolecular decomposition of CL-20.
The 3b intermediate has been found in two simulations. It exhibits behavior similar to the one presented in Figure 3a but has a slightly smaller lifetime at the given temperature. The intermediate, 3c, was found in one trajectory at T ) 1500 K. This is the only example of a species obtained by departure of four nitro groups. Its lifetime was found to be 7 ps. We found three major distinctive channels for subsequent decomposition of CL-20 intermediates (Figure 2). These secondary reactions include (a) fragmentation, (b) oxidation by a NO2 radical, and (c) hydrogen abstraction by a NO2 radical. Fragmentation of the intermediate results in stepwise breaking into HCN molecules. Oxidation by NO2 yields N-methylenformamide derivatives and NO (Figure 2), leading to further decomposition. Methylenformamide derivatives undergo fast decomposition, eliminating one to two HCN molecules (depending on molecular length) and yielding the simplest methylenformamide carbenes. Finally, the CN-COH molecule degrades by the scheme displayed in Figure 4a. The third pathway represents hydrogen abstraction by NO2, resulting in HONO molecules. This is an interesting variation of HONO molecular formation, formed during the later stages of the reaction (τ > 2 ps) by interaction of the NO2 free radical and the intermediate. On the basis of our knowledge, such variation of hydrogen abstraction has not before been described in nitramines, neither experimentally nor theoretically. This finding is in contrast to the suggested concerted intramolecular mechanism of HONO formation in RDX/HMX,25,26,50 where departing NO2 groups capture hydrogen radical. Our results strongly favor the stepwise redox intermolecular process. Our observations are also consistent with the recent experimental study on the interaction of NO2 with hydrocarbon soot51 and the high-level CCSD(T) study on hydrogen abstraction from monomethylhydrazine by NO2.52 Under such extreme conditions, the HONO molecule is not stable and dissociates into NO and OH radicals. We also found evidence related to the possibility of formation of an NO3 (Figure 4b). Similarly, its formation has been also
predicted for RDX by Irikura and Johnson.53 The authors pointed out that, due to several factors, this species has never been detected experimentally. In our case, from 10 simulations, NO3 occurred only once, at T ) 1500 K. Two NO2 radicals combined form NO and NO3. For a short time (
