Hydrogen Transfer in Energetic Materials from ReaxFF and DFT

Mar 28, 2017 - Energetic materials are characterized by fast and complex chemical reactions. It makes them hardly available for kinetic experiments in...
0 downloads 21 Views 3MB Size
Download PDF
Article pubs.acs.org/JPCA

Hydrogen Transfer in Energetic Materials from ReaxFF and DFT Calculations Oleg V. Sergeev†,‡ and Alexey V. Yanilkin*,†,‡ †

Center for Fundamental and Applied Research, Dukhov Research Institute of Automatics, P.O. Box 918, Moscow 101000, Russia Department of Molecular and Chemical Physics, Moscow Insitute of Physics and Technology, Moscow, Russia

‡

ABSTRACT: Energetic materials are characterized by fast and complex chemical reactions. It makes them hardly available for kinetic experiments in relevant conditions and a good target for reactive molecular dynamics simulations. In this work, unimolecular and condensed-phase thermal decomposition of pentaerythritol tetranitrate (PETN) are investigated by ReaxFF molecular dynamics. It is shown that the decomposition kinetics in condensed phase may be described with the activation barrier lower by a factor of 2 than that for isolated molecules. The effect of the intermolecular hydrogen transfer is revealed in condensed phase. Energetic barriers for hydrogen transfer in two energetic materials (methyl nitrate, which is a nitroester as well as PETN, and o-nitrotoluene) are studied with ReaxFF and DFT using nudged elastic band technique. The results indicate that ReaxFF gives significantly lower activation energy for intermolecular hydrogen transfer in nitroesters than different DFT approximations, which explains the molecular dynamics results for PETN.

■

INTRODUCTION ReaxFF interatomic potential1 is extensively used for molecular dynamics simulations of wide range of materials. In particular, it allows large-scale (up to millions of atoms) simulations of energetic materials with atomistic resolution. In energetic materials research, ReaxFF was successfully applied to such problems as shock-wave propagation,2−4 effects of interfaces,5−7 formation and growth of hot spots,8−12 energy transfer, and burning propagation in crystalline explosives.13−15 The bond order concept lying in the basis of ReaxFF simplifies the analysis of chemical reactions. Chemical information was analyzed in a number of works using ReaxFF2,16−18 to determine, for example, final composition or kinetic curves for product formation. It was shown that ReaxFF was successfully parametrized to describe some known mechanisms of energetic materials decomposition quite well.16 The force field was also used to discover reaction pathways in thermal decomposition process.19−21 In the present work we use molecular dynamics simulations with ReaxFF potential to analyze kinetics of thermal decomposition of pentaerythritol tetranitrate (PETN). PETN is an energetic material used both in pure form and in compounds. It is characterized by high sensitivity compared to other common materials such as RDX and HMX. 22 Experimental studies showed the anisotropic sensitivity of PETN crystals to shock loading.23−26 In principle, this sensitivity depends on both the initial thermal response to shock-wave passage, that is, the initial formation of heated regions (hot spots) where the energy of shock wave dissipates, © 2017 American Chemical Society

and the following interplay of the competing processes of heat production in chemical reactions and heat sink in hydrodynamic and heat conduction processes. Both stages are almost unavailable for experimental investigation in realistic initiation conditions, because they are very fast and occur in condensed phase. Microscopic mechanisms of hot-spot formation were studied with molecular modeling methods for both single-crystal9,25,27 and polycrystalline9 PETN. Modeling of the hot-spot growth was also performed for PETN with atomistic resolution.11,13 However, there is no established model of chemical kinetics derived from reactive molecular dynamics and applicable to mesoscale simulations. Work is underway to construct such models for various energetic materials. The exhaustive study of chemical reactions in TATB was performed recently by Long and Chen,28 where more than 450 species and 5500 elementary reactions were identified, each reaction being described with Arrhenius kinetics. The authors applied the model to simulate fast reaction zone in detonation front. For the mesoscale modeling of the hot-spot growth and the onset of deflagration, however, it is impractical to use many intermediate species. Another approach was utilized in the work of Lee et al.29 concerning RDX. All the species were splitted into four categories based on their molecular weight and the role in decomposition process. Reaction rates for transformations of Received: December 30, 2016 Revised: March 23, 2017 Published: March 28, 2017 3019

DOI: 10.1021/acs.jpca.6b13088 J. Phys. Chem. A 2017, 121, 3019−3027

Article

The Journal of Physical Chemistry A

phase simulation until the first reaction event, and bulk phase thermal decomposition to final products. First, we construct an ensemble of MD cells each of which initially contains one PETN molecule in empty space. The geometry of the molecule is initially optimized with the conjugated gradient technique using ReaxFF potential. Then the velocities are assigned to atoms so that the net linear and angular momenta are zero, and after short equilibration period the velocity distribution in each MD cell corresponds to some particular temperature. Initial velocities distribution is constructed with random number generator, and different seed is used for each MD cell, so that the cells turn out to be independent from each other in phase space. After the velocities are set, the simulation continues in NVE ensemble in each cell independently. At some time moments (random and depending on the temperature), dissociation events occur in some cells. If the temperature is not too high, one particular path of decomposition with the highest reaction rate will dominate this stochastic process. The ratio of the cells with nonreacted molecules to the total initial number of MD cells is described by the kinetics with the same kinetic coefficients as those characterizing unimolecular decomposition of PETN (as there is one molecule per MD cell). Thermal decomposition in condensed phase is considered from two perspectives. The kinetics of the dominating initial reaction is obtained independently by the technique similar to the one used for unimolecular decomposition. For every target temperature in the range from 700 to 1300 K 20 MD cells are created, each containing 128 PETN molecules (4 × 4 × 4 unit cells). All cells are led to zero pressure and the temperature of 300 K in NPT calculation for 20 000 timesteps (1 ps). Then the velocities are ascribed to atoms from scratch such that after the equipartitioning of energy each of 20 MD cells has the temperature in the target range. These initial temperature values in 20 cells are very close to each other, while the velocities of corresponding atoms are different. We again calculate the ratio of the systems in which no chemical reactions occurred to the total number of 20, and that ratio is also described by some kinetics. As there are 128 molecules in each MD cell, the reaction rate for this “system decomposition” reaction is taken to be 128 times higher than the reaction rate for actual reactions between molecules. It may be noted in advance that in both cases described above the kinetics for the decrease of “unreacted cells” (MD cells where no reactions occurred since simulation start) number is first-order. In unimolecular decomposition the reaction events are, again, unimolecular, so the first-order kinetics is expected. However, it will keep first-order for the number of unreacted cells in bulk phase even though the reaction itself may well be bimolecular. The reason is that until any chemical event occurred in a cell, it is in the same configuration (contains the same molecules) as any other unreacted cell. So the first event may occur in any cell at any particular time with equal probability, and this probability does not depend on time. Those are the conditions for the first -order kinetics. Finally, thermal decomposition to final products is simulated. In this part we consider only one larger MD cell (1024 PETN molecules, 8 × 8 × 8 unit cells) for each temperature. The procedure of creating the initial high-temperature state is the same as in the previous part. At each moment during the simulation we identify all molecules in the MD cell, and thus for each chemical species there is a kinetic curve of its

these meta-species were determined by functions constructed so that the continuum modeling reproduced molecular simulation results. The technique implemented in the present work also uses a few-component (here a three-component) kinetic model but employs Arrhenius formulation for rate constants. We pay special attention to the decomposition mechanism in condensed phase, and particularly to the possible intra- or intermolecular hydrogen transfer. The process is not uncommon for other materials: several papers30−36 were devoted to theoretical investigation of intramolecular hydrogen transfer and showed that this mechanism could be the first step in gas-phase decomposition due to its relatively low activation energy. In o-nitrotoluene (o-NT) the activation energy varies from 16933,35,36 to 194 kJ/mol30,31 depending on computational details, in TNT it equals 150 kJ/mol,32 in HMX186 and 182 kJ/mol.34,37 Intramolecular hydrogen transfer in condensed phase and on a surface of β-HMX was considered by Sharia et al.38 The presence of surface drastically affects the activation energy and reduces it to 107 kJ/mol, although in bulk the value is equal to 206 kJ/mol for β-HMX. 34,38 Intermolecular process was investigated by Lewis39 in the tight-binding approximation. The obtained value is much larger than that for the intramolecular transfer in condensed phase ca. 300 kJ/mol for δ-HMX versus 206 kJ/mol in the work of Sharia et al.38 We test the accuracy of ReaxFF description of the discovered pathways against density functional theory (DFT) calculations. To facilitate the DFT calculations, we take two other materials, namely, o-nitrotoluene and methyl nitrate (MN), for this part of work. The former is chosen because there are literature data available for the activation barrier of hydrogen transfer in the molecule.30,31,33,35,36 Methyl nitrate is considered because it is a nitroester and is believed to have the kinetic parameters for the studied reactions similar to those in PETN. The rest of the paper is organized as follows. First, we describe atomistic models used to obtain kinetic parameters for thermal decomposition of PETN and to compare ReaxFF with DFT. Then we present the results of molecular dynamics simulations of thermal decomposition of PETN and show that some of them are related to hydrogen transfer process. Third, energetic profiles of hydrogen transfer reactions in o-NT and MN are calculated, which allows estimation of the accuracy of the description of this process with ReaxFF. Then follow the conclusions.

■

CALCULATION DETAILS Molecular Dynamics Models. Molecular dynamics (MD) simulations for PETN with ReaxFF in this work were mostly performed in NVE ensemble. There are numerous works in literature where chemical reactions are observed in a thermostated system.16,17,19,21,40−43 Such an approach usually makes it more straightforward to extract kinetic coefficients from the simulation results, when the temperature is kept constant. However, energetic materials are characterized by large heat effects of typical reactions. A thermostat capable to keep such systems at constant temperature may significantly influence the microscopic dynamics of atoms. We intended to let the local interactions inside groups of atoms be only determined by the interatomic potential, without any effect of the external heat bath. Kinetics of thermal decomposition of PETN is studied with ReaxFF in three approaches: unimolecular decomposition, bulk 3020

DOI: 10.1021/acs.jpca.6b13088 J. Phys. Chem. A 2017, 121, 3019−3027

Article

The Journal of Physical Chemistry A formation and decomposition with time. Some species, when formed, do not take part in further reactions, thus being the final products. We describe the formation of these final products with the chain kinetic model, where the rate of formation of intermediate products is taken from the “first step” modeling described above, while the kinetic parameters for the chain continuation and chain break stages are found independently. In all classical MD simulations in this work we use LAMMPS package44 and the ReaxFF-lg parametrization for energetic materials developed by Liu et al.45 The time step for molecular dynamics is 0.05 fs. In chemical analysis we take instantaneous bond orders into account; two atoms are considered connected if the order of the bond between them is higher than 0.2. Nudged Elastic Band Calculations. The comparison of ReaxFF and DFT descriptions of the hydrogen transfer process is made by energies and barriers of the reactions. The reaction energy ΔrE is calculated as the difference between the energy of the relaxed final and initial states. The activation energy is obtained by the nudged elastic band (NEB) method46,47 with spring constant between parallel replicas of 5 eV/Å2. In NEB calculations we consider two species different than PETN, namely, o-nitrotoluene and methyl nitrate. To model an isolated molecule or two molecules, we constructed a 10 × 10 × 10 Å3 cell containing, correspondingly, a single molecule or two molecules of methyl nitrate. Two molecules of nitrotoluene are considered in 15 × 10 × 15 Å3 cell. Classical NEB calculations are performed by means of LAMMPS44 code. Interatomic potential is ReaxFF-lg.45 Quantum NEB calculations based on DFT are performed by VASP software package.48,49 The projector augmented wave method49,50 is used for electron wave functions. The exchangecorrelation energy is considered in generalized gradient approximation in the PBE form51 and in hybrid functional form of B3LYP.52 Gamma point only is used for k-point mesh. The energy cutoff of plane-wave basis is fixed at 600 eV. The value of 1000 eV provides the same results. For hydrogen dissociation and intermolecular hydrogen transfer spinpolarized calculations are performed.

Figure 1. Part of isolated PETN molecules not dissociated to the moment t. Velocity distribution for each molecule corresponds to 1600 K.

Figure 2. Reaction rate constant dependency on the inverse temperature for unimolecular decomposition.

First-Stage Kinetics in Bulk PETN. In condensed phase, there is no possibility to construct a kinetic model that would include all the possible reactions. A particular mechanism of thermal decomposition is highly variative and gets complicated after only a few initial reactions. However, for the very first stage the kinetic parameters are defined, as the initial geometry and chemical composition are the same for all MD cells. We use the procedure similar to the unimolecular decomposition case to obtain these parameters. Again, the plot of the logarithm of reaction rate constant against inverse temperature gives the activation energy (Figure 3). The slope of the line on the figure corresponds to the activation energy of 82 ± 5 kJ/mol. It is quite different from the unimolecular case and indicates a possibility for a different firststage reaction mechanism. To check for this, we studied the chemical composition after the first reaction in MD cell for each cell independently. In all MD runs, the first reaction appeared to be the intermolecular hydrogen transfer. On Figure 4 there is a visualization of the process. Only molecules participating in the reaction are shown. After the hydrogen transfer, the nitro group of the molecule that lost hydrogen becomes weakly connected to the rest of the molecule, because the oxygen in the same branch is likely to form second bond with the unpaired electron of the carbon. The corresponding decrease of the O−N bond order may be

■

RESULTS Unimolecular Decomposition. We do not count possible isomerisations in the unimolecular case, so the kinetics observed are for the decomposition only. Analysis of the results shows that the most common product for the first step of the decomposition is NO2, with some addition from other pathways. For each particular temperature the decrease in number of nondecomposed molecules with time is exponential. It shows that from the kinetic point of view there is one firstorder process dominating this initial step (Figure 1). The parameter λ of the approximating function on Figure 1 is the reaction rate constant for unimolecular PETN decomposition at the corresponding temperature. The plot of its logarithm against inverse temperature is linear in the considered temperature range of 1210 ÷ 2140 K. It indicates that the decomposition mechanism remains the same over the range (Figure 2). The slope of the line is proportional to the activation energy of the decomposition process. For the obtained data it gives the activation energy of 156 ± 2 kJ/mol, which is in agreement with the known experimental value53 and the results of ab initio calculation of the NO2 detachment barrier.54 3021

DOI: 10.1021/acs.jpca.6b13088 J. Phys. Chem. A 2017, 121, 3019−3027

Article

The Journal of Physical Chemistry A d[E] = −k1[E] − k 2[E][X] dt d[X] = k1[E] + k 2[E][X] − 2k 3[X]2 dt d[P] = 2mk 3[X]2 dt

where [E], [X], [P] are concentrations of the corresponding model species, while k1, k2, and k3 are the reaction rate constants determined by the Arrhenius law ki = Zie−Ea,i /RT

Integration of these equations gives model kinetic curves for E, X, and P. These model curves are determined by the coefficients Zi and Ea,i. For each particular set of the coefficients and the initial temperature T0 we define the approximation error ωset(T0) as

Figure 3. Reaction rate constant dependency on the inverse temperature for the first reaction in bulk phase decomposition. Error bars are within the points.

Nspecies Nsteps

ωset(T0) =

∑ ∑ (Ci ,model(t ) − Ci ,MD(t ))2 i=1

t=0

where Ci,model(t) is the concentration of ith substance at time step t, Ci,MD(t) is the corresponding concentration actually observed in MD simulation. Concentrations are normalized per initial number of PETN molecules. As the simulation is performed in NVE ensemble, the heat produced in the decomposition increases the system temperature. Rate constants depend strongly on temperature. To account for this dependency, we interpolated the temperature calculated in MD with cubic splines and then substituted temperature values from these splines while integrating kinetic equations. The number of knots the splines are interpolated through is much lower than the number of steps where temperature is logged. Knots number is chosen for each MD run independently so that the resulting curve conserves the general shape of the one from molecular dynamics but smoothes the fluctuations. We find the optimal set of kinetic constants by minimizing the total approximation error Ωset = ∑TTmax ωset(T0). Starting 0=Tmin parameters for the first reaction were taken from the results of the first stage simulation (see above). In the process of minimizing Ωset, however, these parameters were allowed to change. The resulting reaction rate constants are

Figure 4. Atomistic picture of hydrogen transfer in bulk PETN. Gray is carbon, bluenitrogen, redoxygen, and whitehydrogen. Bonds have width proportional to the corresponding bond order and are colored accordingly, with pure blue corresponding to the bond order of 1.0 and pure greento the bond order of 2.0. Bond orders are instantaneous and calculated by ReaxFF. Only the molecules participating in the process are shown.

seen in the figure, as the bond becomes “thinner” and eventually breaks. Free nitro group, formed with activation barrier different from the unimolecular case, may then participate in further reactions. The next step is to check whether the observed activation energy for the hydrogen transfer influences the overall kinetics of PETN decomposition. Kinetic Model for Products Formation in Bulk PETN. As was said before, the precise scheme of chemical reactions in bulk energetic material is very complex. Instead, we consider a simplified three-component model with the explosive E, intermediate products X, and final products P. Thermal decomposition of the explosive goes by chain mechanism

4

k1 = 1.5 × 1013e−7.8·10

s −1

4

J/mol/RT

molecules·s−1

5

J/mol/RT

molecules·s−1

k 2 = 4.5 × 1024 × 1013e−9.0·10 k 3 = 3.5 × 1024 × 1013e−1.7·10

J/mol/RT

m = 1.47. The number 1024 is the initial number of PETN molecules in the MD cell. It appears in the reaction rate constants for the second-order reactions because concentrations are written in units of [Nmolecules/1024]. Figure 5 shows the comparison of the model and MD concentrations of PETN, intermediate species, and a final product (N2). The results of the approximation show that the proposed chainlike model can successfully describe the decomposition of PETN and formation of final products. The fact that the kinetic curve for X concentration from the model goes close to the total concentration of NO2 and HONO is not provided by the minimization procedure, as X concentration is not included in

E→X E + X → 2X X + X → 2mP

where m is the number of product molecules per one initial explosive molecule. The corresponding kinetic equations for model components are 3022

DOI: 10.1021/acs.jpca.6b13088 J. Phys. Chem. A 2017, 121, 3019−3027

Article

The Journal of Physical Chemistry A

Figure 5. Bulk-phase PETN decomposition described with two-stage model. All x axes are time in picoseconds; all y axes are the number of molecules divided by the initial number of PETN molecules in the cell. Model concentrations are shown with thick lines, while MD concentrations with thinner ones. Final product is N2 (blue lines). For intermediate products in MD the sum concentration of NO2 + HONO is shown (dotted purple lines). The X amount (thick purple lines) was not fitted to these MD dependencies. Fit is made against PETN (red lines) and molecular nitrogen concentrations. Green lines show the contributions to ωset (upscaled for visibility).

the calculation of Ωset. This coincidence may indicate that the two species indeed serve as reaction propagators in N2 formation. Another result is that in the overall process of thermal decomposition of PETN simulated with ReaxFF, the first step shows very similar kinetics to the hydrogen transfer reaction. We did not follow in detail each particular N2 formation event in condensed phase, but the kinetics suggests that intermolecular hydrogen transfer is the common initiation reaction in thermal decomposition of PETN. Thus, accurate description of this reaction is necessary for the reactive MD simulations of PETN with the ReaxFF-lg parametrization to correctly describe the whole thermal decomposition process. To check it we compare energetic barriers for several hydrogen transfer reactions described with ReaxFF and DFT using NEB technique. Intramolecular Hydrogen Transfer. Intramolecular hydrogen transfer is a straightforward way to compare ReaxFF with DFT, as it does not involve degrees of freedom related to the mutual configuration of several molecules. Recent DFT results show that this process can be the first step in gas-phase thermal decomposition of o-NT.31,36 So, first, we compare our results with published data. Then, methyl nitrate (MN) is considered as a simpler example of a nitroester energetic material, containing R-CH2−O−NO2 part and similar in that to PETN. For methyl nitrate “R” in the formula is a hydrogen atom. o-Nitrotoluene. The atomic configurations of initial (IS), transition (TS), and final (FS) states, obtained with DFT, are shown in Figure 6. Our PBE and B3LYP results are close to each other and to the literature data.30,31,36 The difference

Figure 6. Energy surface profiles along reaction pathway for the isomerization (intramolecular hydrogen transfer) of o-NT. Atomic structures of IS, TS, and FS are shown together with isosurfaces of electron density at the level of 0.09 el/Å3. The stars are the results from Khrapkovskii et al.,36 while the cross is from Chen et al.31

between the results is ∼0.005 Å for bond lengths in IS, TS, and FS. The configurations obtained with ReaxFF are quite different, especially for transition and final states. The FS has the dihedral angle ∠(CNOH) ≈ 90°, while in our DFT calculations it is equal to 13°. 3023

DOI: 10.1021/acs.jpca.6b13088 J. Phys. Chem. A 2017, 121, 3019−3027

Article

The Journal of Physical Chemistry A Figure 6 shows the potential energy as a function of the C− H bond length along the reaction pathways, obtained by NEB. There are a lot of points with similar energy and bond length at the beginning of the curve. These points correspond to rotation of the methyl group around C−C bond. Then DFT points drastically go up with almost fixed value of C−H bond length. The increase of energy is due to an oxygen getting closer to hydrogen atom. In TS the hydrogen is weakly bonded with carbon, as it is shown by distribution of the electron density on Figure 6. ReaxFF points are only close to those obtained with DFT at the beginning. Then the behavior is quite different as well as the mechanism. The hydrogen atom moves to the oxygen just by stretching of C−H bond without displacement of the oxygen. It can be seen from bond lengths in TS configurations: in DFT d(O−H) equals 1.13 Å when d(C−H) ≈ 1.5 Å, while in ReaxFF d(O−H) ≈ 1.6 Å. One possible reason for the behavior is the long-range interaction between hydrogen and oxygen present in ReaxFF. The reaction energies, calculated at different levels of theory, are the following: DFT gives 139 kJ/mol (PBE) and 152 kJ/ mol (B3LYP), while ReaxFF72 kJ/mol. The ReaxFF value is very different due to another configuration in final state. Activation energies are closer to each other in spite of different mechanisms given by ReaxFF and DFT: 155 kJ/mol (PBE) and 182 kJ/mol (B3LYP), 141 kJ/mol (ReaxFF). Our B3LYP values are between different numbers present in literature.30,31,33,35,36 Methyl Nitrate. The main mechanism of gas-phase thermal decomposition of methyl nitate is O−NO2 bond breaking with activation energy of ∼150−170 kJ/mol.55 However, hydrogen transfer from methyl to nitro group is likely to have a barrier similar to that in o-nitrotoluene. Initial, transition, and final states for the mechanism are shown in Figure 7. It is worth noting that the final state in DFT (both PBE and B3LYP) is not stable and dissociates to give formaldehyde and trans form of nitrous acid. The dissociation already takes place in TS, which can be seen by the absence of electron density between oxygen

and nitrogen atoms (Figure 7). In ReaxFF the aci form of methyl nitrate is stable. Energy profiles along the reaction pathway are shown in Figure 7. As in the case of o-NT, the first steps of the reaction are rotation of methyl group and oxygen displacement. In TS the hydrogen atom is quite close to the oxygen: d(O−H) = 1.25 Å. TS from ReaxFF calculations is characterized by strongly stretched C−H bond. Such behavior is close to the ReaxFF results for NT. In DFT the final state dissociates into CH2O and HNO2. If we calculate the energy of the same products with ReaxFF, it is reduced by 30 kJ/mol compared to ReaxFF final state, but is still larger than PBE value by 47 kJ/ mol. As for NT, in spite of different mechanisms and reaction energies, the activation energies are close: from ReaxFF we have 146 kJ/mol, from DFT 155 kJ/mol (PBE) and 186 kJ/ mol (B3LYP). The reaction energy can be calculated based on the 0 experimental values of heats of formation: Δ rH 298 = 0 0 0 ΔfH298 (CH2O) + ΔfH298(HNO2) − ΔfH298(CH3ONO2). Heats of formation of gas phases are the following: ΔfH0298(CH2O) = −116 kJ/mol, ΔfH0298(HNO2) = −79 kJ/ mol, ΔfH0298(CH3ONO2) = −122 kJ/mol.56 The resulting value ΔrH0298 = −73 kJ/mol is close to our B3LYP result ΔrE = −50 kJ/mol. Different DFT descriptions (PBE and B3LYP) give both the same mechanism and negative energy difference for intramolecular hydrogen transfer. As opposed to this the use of ReaxFF results in a different mechanism and different final state. The value of activation energy, however, is similar to that of PBE and differs from B3LYP result by less than 40 kJ/mol. Intermolecular Hydrogen Transfer. Thermal decomposition of nitro compounds in condensed phase has much more complex mechanism due to the interaction of a large number of molecules. It is in general more complicated to calculate activation energies in such systems. We consider a simple example of condensed-phase reactions, namely, intermolecular hydrogen transfer between two molecules in vacuum. The positions of molecules are chosen to provide close location of methyl and nitro groups. The initial state is obtained by relaxation of such configuration. Potential energy of such system is sligthly lower than doubled value of that of a single molecule due to weak binding. The final state is produced from the initial one by moving a hydrogen atom to the oxygen of nitro group and then relaxing the system. The choice of initial configuration of two molecules can affect the activation energy. It is possible that there is another reaction pathway with smaller barrier. However, in our test calculations other configurations provided similar activation energies. In our opinion the possible variation should not be large and does not influence our conclusions. o-Nitrotoluene. Figure 8 shows energy profiles and consequent snapshots of intermolecular hydrogen transfer in NT. As opposed to the intramolecular process in NT, the mechanism in DFT and ReaxFF calculations is the same. Hydrogen comes close to oxygen by C−H bond stretching. The energy profiles are, however, very different. Transition state as well as final state in ReaxFF have energy lower than that in DFT by 100−120 kJ/mol. It is expected from chemical intuition that in comparison to the intramolecular transfer the activation and reaction energies should be higher due to radicals formation in products. It is so for DFT calculation, but in ReaxFF the final state with two radicals has less energy than two isomers of NT. The difference is ∼86 kJ/mol and is too

Figure 7. Energy surface profiles along reaction pathway for the intramolecular hydrogen transfer in methyl nitrate. Atomic structures of IS, TS, and FS are shown together with isosurfaces of electron density at the level of 0.1 el/Å3. 3024

DOI: 10.1021/acs.jpca.6b13088 J. Phys. Chem. A 2017, 121, 3019−3027

Article

The Journal of Physical Chemistry A

Figure 9. Energy surface profiles along reaction pathway for the intermolecular hydrogen transfer in methyl nitrate. Atomic structures of IS, TS, and FS are shown together with isosurfaces of electron density at the level of 0.08 el/Å3.

Figure 8. Energy surface profiles along reaction pathway for the intermolecular hydrogen transfer in NT. Atomic structures of IS, TS, and FS are shown together with isosurfaces of electron density at the level of 0.1 el/Å3.

The reason for the difference between ReaxFF and DFT is in the description of the C−H and O−H bond energies, like in NT. Calculated DFT values of the C−H bond energy are 320 kJ/mol (PBE) and 300 kJ/mol (B3LYP), from ReaxFF264 kJ/mol; the O−H bond energy in DFT is 165 kJ/mol (PBE) and 187 kJ/mol (B3LYP), in ReaxFF the O−H bond energy is equal to 273 kJ/mol, which is ∼100 kJ/mol larger than DFT results.

high to be due to the attractive interaction between two radicals. To understand the reason for such a difference between ReaxFF and DFT, let us consider the intermolecular hydrogen transfer in details. The whole process can be splitted in two steps. The first one is the dissociation of a hydrogen from methyl group, and the second one is its association with nitro group. The reaction energy of the whole process is the sum of reaction energies of the two steps. The reaction energy of the first step is the C−H bond energy in the methyl group of NT. Calculated DFT values are 385 kJ/mol (PBE) and 392 kJ/mol (B3LYP), while ReaxFF gives 335 kJ/mol. The second step is the association of the hydrogen to the nitro group of NT, which results in the additional O−H bond formation. The O−H bond energy calculated with DFT is 200 kJ/mol (PBE) and 218 kJ/ mol (B3LYP), in ReaxFF calculation280 kJ/mol, which is ∼60−80 kJ/mol larger than DFT results. Summing up, in ReaxFF the total reaction energy is lower by 120 kJ/mol, and the activation energy is lower by 100 kJ/mol (186 kJ/mol (PBE), 205 kJ/mol (B3LYP), 105 kJ/mol (ReaxFF)) compared to DFT. Methyl Nitrate. Intermolecular hydrogen transfer in MN has similar mechanism to the intramolecular one (Figure 9). The final state of MN without hydrogen dissociates into formaldehyde and nitrogen dioxide. In summary in DFT the intermolecular hydrogen transfer is endothermic with the reaction energy of 113 kJ/mol (B3LYP) and 145 kJ/mol (PBE), while in ReaxFF the reaction is exothermic with the reaction energy −3 kJ/mol. Such a difference in reaction energies is also reflected in activation barriers. In DFT the values are quite high: 240 kJ/mol (PBE) and 270 kJ/mol (B3LYP), while in ReaxFF the activation energy is only 78 kJ/ mol. The last value is close to the activation energy obtained in molecular dynamics simulations for the first reaction of PETN thermal decomposition in condensed phase.

■

CONCLUSIONS Kinetics of thermal decomposition of PETN is determined from ReaxFF molecular dynamics simulations of single molecules and bulk material at elevated temperatures. The activation energies of the first step reactions for monomolecular and condensed-phase processes differ by the factor of 2: they equal, correspondingly, 156 and 82 kJ/mol. To explain the discrepancy we determine what is the initial thermal decomposition reaction in condensed-phase PETN, and it proves to be the intermolecular hydrogen transfer. The hydrogen transfer process is investigated in detail with ReaxFF and DFT in o-nitrotoluene and methyl nitrate to test ReaxFF accuracy. For the intramolecular hydrogen transfer, different DFT approaches (PBE and B3LYP) provide the same mechanism and close energies. Activation barriers calculated with DFT for o-nitrotoluene are in agreement with published data. The use of ReaxFF results in a different mechanism and different final state. The values of activation energy calculated with ReaxFF for these reactions are similar to the ones given by PBE (the difference is less than 15 kJ/mol) and lower than the ones in B3LYP approximation (by ∼40 kJ/mol), still being close to 150 kJ/mol for methyl nitrate. Very different results are obtained for the intermolecular hydrogen transfer. Activation and reaction energies in ReaxFF are lower than those in DFT by ∼100 kJ/mol. The reason is the incorrect description of bond energies: the C−H bond in methyl group and the O−H bond in radicals. Such low 3025

DOI: 10.1021/acs.jpca.6b13088 J. Phys. Chem. A 2017, 121, 3019−3027

Article

The Journal of Physical Chemistry A

(14) Zhou, T.; Song, H.; Liu, Y.; Huang, F. Shock initiated thermal and chemical responses of HMX crystal from ReaxFF molecular dynamics simulation. Phys. Chem. Chem. Phys. 2014, 16, 13914− 13931. (15) Joshi, K.; Losada, M.; Chaudhuri, S. Intermolecular energy transfer dynamics at a hot-spot interface in RDX crystals. J. Phys. Chem. A 2016, 120, 477−489. (16) Strachan, A.; Kober, E. M.; van Duin, A. C.; Oxgaard, J.; Goddard, W. A. Thermal decomposition of RDX from reactive molecular dynamics. J. Chem. Phys. 2005, 122, 4502. (17) Zhou, T.-T.; Huang, F.-L. Effects of defects on thermal decomposition of HMX via ReaxFF molecular dynamics simulations. J. Phys. Chem. B 2011, 115, 278−287. (18) Furman, D.; Kosloff, R.; Dubnikova, F.; Zybin, S. V.; Goddard, W. A., III; Rom, N.; Hirshberg, B.; Zeiri, Y. Decomposition of condensed phase energetic materials: Interplay between uni-and bimolecular mechanisms. J. Am. Chem. Soc. 2014, 136, 4192−4200. (19) van Duin, A. C.; Zeiri, Y.; Dubnikova, F.; Kosloff, R.; Goddard, W. A. Atomistic-scale simulations of the initial chemical events in the thermal initiation of triacetonetriperoxide. J. Am. Chem. Soc. 2005, 127, 11053−11062. (20) Li, Y.; Kalia, R. K.; Nakano, A.; Nomura, K.-I.; Vashishta, P. Multistage reaction pathways in detonating high explosives. Appl. Phys. Lett. 2014, 105, 204103. (21) Yan, Q.-L.; Zeman, S.; Sánchez Jiménez, P.; Zhang, T.-L.; PérezMaqueda, L.; Elbeih, A. The mitigation effect of synthetic polymers on initiation reactivity of CL-20: Physical models and chemical pathways of thermolysis. J. Phys. Chem. C 2014, 118, 22881−22895. (22) Licht, H.-H. Performance and sensitivity of explosives. Propellants, Explos., Pyrotech. 2000, 25, 126−132. (23) Dick, J. Effect of crystal orientation on shock initiation sensitivity of pentaerythritol tetranitrate explosive. Appl. Phys. Lett. 1984, 44, 859−861. (24) Dick, J.; Mulford, R.; Spencer, W.; Pettit, D.; Garcia, E.; Shaw, D. Shock response of pentaerythritol tetranitrate single crystals. J. Appl. Phys. 1991, 70, 3572−3587. (25) Dick, J.; Ritchie, J. Molecular mechanics modeling of shear and the crystal orientation dependence of the elastic precursor shock strength in pentaerythritol tetranitrate. J. Appl. Phys. 1994, 76, 2726− 2737. (26) Yoo, C.; Holmes, N.; Souers, P.; Wu, C.; Ree, F.; Dick, J. Anisotropic shock sensitivity and detonation temperature of pentaerythritol tetranitrate single crystal. J. Appl. Phys. 2000, 88, 70−75. (27) Zybin, S. V.; Goddard, W. A., III; Xu, P.; van Duin, A. C.; Thompson, A. P. Physical mechanism of anisotropic sensitivity in pentaerythritol tetranitrate from compressive-shear reaction dynamics simulations. Appl. Phys. Lett. 2010, 96, 081918. (28) Long, Y.; Chen, J. Theoretical study of the reaction kinetics and the detonation wave profile for 1, 3, 5-triamino-2, 4, 6-trinitrobenzene. J. Appl. Phys. 2016, 120, 185902. (29) Lee, K.; Joshi, K.; Chaudhuri, S.; Stewart, D. S. Mirrored continuum and molecular scale simulations of the ignition of highpressure phases of RDX. J. Chem. Phys. 2016, 144, 184111. (30) Il’ichev, Y. V.; Wirz, J. Rearrangements of 2-nitrobenzyl compounds. 1. Potential energy surface of 2-nitrotoluene and its isomers explored with ab initio and density functional theory methods. J. Phys. Chem. A 2000, 104, 7856−7870. (31) Chen, S. C.; Xu, S. C.; Diau, E.; Lin, M. C. A computational study on the kinetics and mechanism for the unimolecular decomposition of o-nitrotoluene. J. Phys. Chem. A 2006, 110, 10130−10134. (32) Cohen, R.; Zeiri, Y.; Wurzberg, E.; Kosloff, R. Mechanism of thermal unimolecular decomposition of TNT (2,4,6-trinitrotoluene): A DFT study. J. Phys. Chem. A 2007, 111, 11074−11083. (33) Fayet, G.; Joubert, L.; Rotureau, P.; Adamo, C. A theoretical study of the decomposition mechanisms in substituted o-nitrotoluenes. J. Phys. Chem. A 2009, 113, 13621−13627.

activation energies strongly influence the kinetics of thermal decomposition of the materials in condensed phase. As an example, in PETN the first reaction is the intermolecular hydrogen transfer with the activation energy of only 82 kJ/mol. The same effect can be expected in other energetic materials containing nitro groups and C−H bonds.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Oleg V. Sergeev: 0000-0002-1416-6357 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank A. G. Shamov, Prof. G. M. Khrapkovskii (KNRTU, Kazan, Russia), and Dr. R. V. Tsyshevsky (University of Maryland) for helpful discussions.

■

REFERENCES

(1) Van Duin, A. C.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A reactive force field for hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (2) Strachan, A.; van Duin, A. C.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A., III Shock waves in high-energy materials: The initial chemical events in nitramine RDX. Phys. Rev. Lett. 2003, 91, 098301. (3) Nomura, K.-I.; Kalia, R. K.; Nakano, A.; Vashishta, P.; van Duin, A. C.; Goddard, W. A., III Dynamic transition in the structure of an energetic crystal during chemical reactions at shock front prior to detonation. Phys. Rev. Lett. 2007, 99, 148303. (4) Budzien, J.; Thompson, A. P.; Zybin, S. V. Reactive molecular dynamics simulations of shock through a single crystal of pentaerythritol tetranitrate. J. Phys. Chem. B 2009, 113, 13142−13151. (5) An, Q.; Zybin, S. V.; Goddard, W. A., III; Jaramillo-Botero, A.; Blanco, M.; Luo, S.-N. Elucidation of the dynamics for hot-spot initiation at nonuniform interfaces of highly shocked materials. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 220101. (6) An, Q.; Goddard, W. A., III; Zybin, S. V.; Jaramillo-Botero, A.; Zhou, T. Highly shocked polymer bonded explosives at a nonplanar interface: Hot-spot formation leading to detonation. J. Phys. Chem. C 2013, 117, 26551−26561. (7) An, Q.; Goddard, W. A., III; Zybin, S. V.; Luo, S.-N. Inhibition of hotspot formation in polymer bonded explosives using an interface matching low density polymer coating at the polymer-explosive interface. J. Phys. Chem. C 2014, 118, 19918−19928. (8) Nomura, K.-I.; Kalia, R. K.; Nakano, A.; Vashishta, P. Reactive nanojets: Nanostructure-enhanced chemical reactions in a defected energetic crystal. Appl. Phys. Lett. 2007, 91, 183109. (9) Cai, Y.; Zhao, F.; An, Q.; Wu, H.; Goddard, W., III; Luo, S. Shock response of single crystal and nanocrystalline pentaerythritol tetranitrate: Implications to hotspot formation in energetic materials. J. Chem. Phys. 2013, 139, 164704. (10) Joshi, K. L.; Chaudhuri, S. Reactive simulation of the chemistry behind the condensed-phase ignition of RDX from hot spots. Phys. Chem. Chem. Phys. 2015, 17, 18790−18801. (11) Shan, T.-R.; Wixom, R. R.; Thompson, A. P. Extended asymmetric hot region formation due to shockwave interactions following void collapse in shocked high explosive. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 94, 054308. (12) Furman, D.; Kosloff, R.; Zeiri, Y. Effects of nanoscale heterogeneities on the reactivity of shocked erythritol tetranitrate. J. Phys. Chem. C 2016, 120, 11306. (13) Sergeev, O. V.; Yanilkin, A. V. Molecular dynamics simulation of combustion front propagation in a PETN single crystal. Combust., Explos. Shock Waves 2014, 50, 323−332. 3026

DOI: 10.1021/acs.jpca.6b13088 J. Phys. Chem. A 2017, 121, 3019−3027

Article

The Journal of Physical Chemistry A

(56) Ray, J.; Ogg, R. The heat of formation of methyl nitrate. J. Phys. Chem. 1959, 63, 1522−1523.

(34) Sharia, O.; Kuklja, M. M. Modeling thermal decomposition mechanisms in gaseous and crystalline molecular materials: Application to β-HMX. J. Phys. Chem. B 2011, 115, 12677−12686. (35) Nikolaeva, E. V.; Shamov, A. G.; Khrapkovskii, G. M. Secondary processes in the mechanism of gas-phase monomolecular destruction of o-nitrotoluene. Russ. J. Gen. Chem. 2014, 84, 2076−2078. (36) Khrapkovskii, G.; Nikolaeva, E.; Egorov, D.; Chachkov, D.; Shamov, A. Energy barriers to gas-phase unimolecular decomposition of mono- and dinitrotoluenes. Russ. J. Org. Chem. 2016, 52, 791−805. (37) Chakraborty, D.; Muller, R. P.; Dasgupta, S.; Goddard, W. A., III Mechanism for unimolecular decomposition of HMX (1,3,5,7tetranitro-1,3,5,7-tetrazocine), an ab initio study. J. Phys. Chem. A 2001, 105, 1302−1314. (38) Sharia, O.; Tsyshevsky, R.; Kuklja, M. M. Surface-accelerated decomposition of δ-HMX. J. Phys. Chem. Lett. 2013, 4, 730−734. (39) Lewis, J. P. Energetics of intermolecular HONO formation in condensed-phase octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX). Chem. Phys. Lett. 2003, 371, 588−593. (40) Rom, N.; Zybin, S. V.; van Duin, A. C.; Goddard, W. A., III; Zeiri, Y.; Katz, G.; Kosloff, R. Density-dependent liquid nitromethane decomposition: Molecular dynamics simulations based on ReaxFF. J. Phys. Chem. A 2011, 115, 10181−10202. (41) Rom, N.; Hirshberg, B.; Zeiri, Y.; Furman, D.; Zybin, S. V.; Goddard, W. A., III; Kosloff, R. First-principles-based reaction kinetics for decomposition of hot, dense liquid TNT from ReaxFF multiscale reactive dynamics simulations. J. Phys. Chem. C 2013, 117, 21043− 21054. (42) Han, S.-P.; van Duin, A. C.; Goddard, W. A., III; Strachan, A. Thermal decomposition of condensed-phase nitromethane from molecular dynamics from ReaxFF reactive dynamics. J. Phys. Chem. B 2011, 115, 6534−6540. (43) Zhou, T.; Liu, L.; Goddard, W. A., III; Zybin, S. V.; Huang, F. ReaxFF reactive molecular dynamics on silicon pentaerythritol tetranitrate crystal validates the mechanism for the colossal sensitivity. Phys. Chem. Chem. Phys. 2014, 16, 23779−23791. (44) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 1995, 117, 1−19. (45) Liu, L.; Liu, Y.; Zybin, S. V.; Sun, H.; Goddard, W. A., III ReaxFF-lg: Correction of the ReaxFF reactive force field for London dispersion, with applications to the equations of state for energetic materials. J. Phys. Chem. A 2011, 115, 11016−11022. (46) Mills, G.; Jonsson, H.; Schenter, G. K. Reversible work transition state theory: Application to dissociative of hydrogen. Surf. Sci. 1995, 324, 305. (47) Henkelman, G.; Jónsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J. Chem. Phys. 2000, 113, 9978−9985. (48) Kresse, G.; Furthmuller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169. (49) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 1758. (50) Blochl, P. E. Projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953. (51) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (52) Becke, A. Density functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (53) Lee, J.-S.; Hsu, C.-K.; Chang, C.-L. A study on the thermal decomposition behaviors of PETN, RDX, HNS and HMX. Thermochim. Acta 2002, 392, 173−176. (54) Tsyshevsky, R. V.; Sharia, O.; Kuklja, M. M. Molecular theory of detonation initiation: Insight from first principles modeling of the decomposition mechanisms of organic nitro energetic materials. Molecules 2016, 21, 236. (55) Manelis, G.; Nazin, G.; Rubtsov, Y. I.; Strunin, V. Thermal Decomposition and Combustion of Explosives and Powders; Nauka: Moscow, Russia, 1996. 3027

DOI: 10.1021/acs.jpca.6b13088 J. Phys. Chem. A 2017, 121, 3019−3027