Molecular dynamics simulations of the decomposition and Us–Up relationship of RDX molecular crystal subjected to high velocity impact

1D shock wave propagation in RDX crystal is simulated using molecular dynamics (MD) with the ReaxFF potential, a forcefield that can simulate chemical reactions in C, H, N, O systems. High velocity impact simulations in RDX molecular crystal along the X , Y and Z axes have been carried out to (i) obtain the relation between the shock velocity ( U s ) versus the particle velocity ( U p ) and (ii) to study the decomposition products of RDX under shock conditions. The U s – U p shows a reasonable match with experiments. It is seen that crystal RDX decomposes only at impact velocities greater than 2.0 km/s and the decomposition products observed favours two of the three different decomposition pathways that have been proposed earlier for RDX.


Introduction
The behaviour of condensed high energy materials under shock loading is determined by the coupling between shock loading, induced mechanical response and chemical reactions. It is desirable to predict coupled mechanical and chemical response of such materials to shock loading. The decomposition of these materials can proceed through a complex set of reaction pathways, yielding a large number of intermediates and final products. The entire process can be studied using multispecies molecular dynamics P. Pahari poonam.pahari@gmail.com; poonamp@barc.gov.in A.D.P. Rao adprao18@gmail.com M. Warrier manoj.warrier@gmail.com; manojwar@barc.gov.in 1 (MD) simulations provided an interatomic potential exists that can permit bond breaking and formation. Ab initio MD codes perform quantum mechanical calculation to determine inter-atomic forces and atomic charge states at each time-step, but at high computational cost [1]. An alternative technique is to use the reactive force field potential such as REBO [2], COMB3 [3] and ReaxFF [4]. In this work, we have used ReaxFF potential. ReaxFF has been parameterised using fits to quantum chemical data obtained from DFT calculations. It is capable of simulating bond breaking and bond formation of C, H, N, O systems. It is known that shock compression and heating of the condensed energetic materials lead to vibrational excitation of energetic molecules, followed by chemical reactions, eventually leading to a final stable set of decomposition products [5].
Experimental efforts to understand the molecular level decomposition mechanisms of RDX have been widely reported. The infrared multiphoton dissociation of RDX molecular beam was studied using time of flight velocity spectra (TOFVS) to find both, the primary and secondary decomposition products of RDX [6]. Zhao et. al have estimated in their experiments that a laser pulse energy deposition of 80 kcal/mol makes the concerted RDX ring breaking energetically accessible. N − NO 2 bond breaking is observed in the ultraviolet photolysis experiments [7][8][9][10]. Thermogravity modulated beam spectroscopy by Behrens and Bulusu states the formation of oxy-striazine by elimination of HONO and HNO [11]. The concerted dissociation of HONO through five member cyclic transition state is observed by Capellos et. al [12] through photo-dissociation of RDX using a 248-nm KrF laser.
The variations of the shock velocity (U s ) as a function of the particle velocity (U p ) are a useful parameter that is measured experimentally [13] and helps in establishing the validity of the equation of state (EOS) of the material [14]. The U s -U p relation depends upon the structural configuration of the atomic system and can be different along different crystallographic directions for molecular crystals like RDX which are not isotropic. Earlier simulations have established the U s -U p for RDX using the ReaxFF potential along X direction of the RDX crystal [15]. The ReaxFF force field was developed for simulation of large-scale chemically reactive systems. It consists of several hundred of parameters derived from quantum chemical calculations. In this paper, we have used the ReaxFF potential, specifically the parameter set in [5], distributed with the LAMMPS code [16] to carry out MD simulations of shocks in α − RDX along the X, Y and Z axes. We obtain the U s -U p relation along the three directions for crystal RDX. We observe shock-induced decomposition for U p > 2 km/s. Chakraborty et al. have reported three possible primary decompostion pathways for RDX [17,18]. Pahari et al [19] have carried out kinetic Monte Carlo (KMC) simulations to study which of these proposed pathways dominate at different temperatures and pressures. The observed dominant pathway of RDX from our MD simulations is presented and compared with the three proposed pathways. In the next section, details of the MD simulation are described. The "Results and discussion" section presents the results and discussions. Finally, the conclusions are presented.

MD simulations of 1D shock
In order to validate the potential, we compute the bulk modulus and shear modulus of RDX crystal using LAMMPS. The bulk and shear modulus computation was carried out using single unit cell of RDX crystal. Table 1 shows the values computed for bulk modulus and shear modulus from the simulations. The average sound speed for the anisotropic crystal can be calculated using the following equation [20].
where V avg is the average sound speed, V l is the longitudinal sound speed, V t is the transverse velocity and B and G are the bulk and shear modulus respectively. Using this formula, we have calculated the average sound speed, we have got the value of 2.51 km/s. The sound speed from published experiments and simulations is shown in Table 2 for comparison. RDX unit cell of size 13.314, 11.86 and10.816Å containing 8 molecules with initial density 1.75 g/cc is considered in the simulations. Figures 1, 2 and 3 show the arrangement of atoms in a unit cell along the X, Y and Z directions using the visual molecular dynamics package (VMD) [29]. In this, C, H, N and O atoms are denoted by cyan, white, blue and red colour balls respectively. Note that the atomic/molecular arrangements are different along the three directions. Hence, it is necessary to separately study impact along the different directions. In the present work, we have chosen to study impact along the X, Y and Z directions. Impact simulations of two similar sized α-RDX single crystals (25 × 1 × 1 and 25 × 2 × 2 unit cells, consisting of 8400 and 33600 atoms respectively) have been carried out. The system was minimised using conjugate gradient method. Hence, the initial temperature is 0 K. Periodic boundary conditions are implemented in the transverse directions to the impact. Shrink-wrapped boundary conditions are specified along the direction of the shock. Shrink-wrapped boundary conditions imply that as atoms cross a boundary, the simulation volume increases to accommodate the atom inside the volume. This allows the shock to reflect from the free surface and a rarefaction wave to proceed back into the crystal. Initially, the two crystals are assigned equal and opposite velocities, along the longitudinal direction, directed towards each other, such that their sum is equal to the desired impact velocity. Impact velocities of 1, 2, 3, 4, 5 and 6 km/s are simulated along the X, Y and Z directions. These impact velocities imply particle  Variation of U s vs U p from our shock simulations along X, Y and Z direction for the 50 × 1 × 1 and its comparison with simulations [5] and experiments [13,14]. The data from 50 × 2 × 2 overlaps with 50 × 1 × 1 and is not plotted To calculate the shock velocity, the system is divided into 50 spatial bins along the impact direction and the density and stress tensor in each bin are computed by sampling over 1 fs and averaging over 10 fs. Due to this spatial and temporal averaging of the pressure, the computed shock velocity will have an error given by δl l + δt t × U s . For 50 unit cells, the bin size is 0.02 for normalised length 1. The uncertainty in time δt = 0.1 fs averaged over t = 10 fs, so the error in our computation of U s is equal to 0.02 This data is post-processed by visualising the shock reaching the free surface ends of the crystal by plotting the pressure in each bin from the stress tensor data. Due to shock compression, the crystal size also shrinks. The shock velocity is calculated by looking at the distance travelled by the shock and the time taken for it to travel at different snap shots. The shock velocity is the slope of the distance travelled by the shock as a function of time.

Results for the U s -U p relation
Typical propagation of a shock in RDX crystal for 6 km/s impact velocity is shown in Fig. 4a. The six frames from top to bottom show the system configuration at 0.5, 1, 1.5, 2, 2.5 and 3 ps respectively. In Fig. 4b, corresponding shock pressure as a function of normalised system length is plotted. In Fig. 5, the U s -U p results from the 50 × 1 × 1 shock simulations along the X, Y and Z directions are compared with published simulation results [5] and from experiments (computed from isothermal volume compression data) [13,14]. We have fitted the shock velocity data from our impact simulations and the results are presented in Table 3. The fit for U s from our X direction impact simulations matches reasonably well with [13] except it predicts higher value of intercept of about 12% and lower value of slope about 3.1%. The results of the 50 × 1 × 1 and 50 × 2 × 2 shock simulations are similar implying that finite size effects along the transverse direction do not affect the simulation results.
It is seen from Fig. 5 that the U s -U p relation depends on the direction of impact. This can be expected for anisotropic molecular crystals due to the different structural arrangements of the atoms along the different directions. Note that in Fig. 5, for impacts along the Y and Z directions respectively, shows two different slopes, for U p > 2 km/s and U p < 2 km/s. The slope of the U s − U p relation depends on the compression in the direction of shock propagation and the arrangement of atoms along that direction. Structural phase transitions change the slope of the U s − U p relation. This is because the slope depends on the compression in the direction of shock propagation and thereby the arrangement  of atoms along that direction. In case of molecular crystal like RDX, above a certain pressure, the bonds at the periphery of the CN ring, for example the N − NO 2 bond can orient itself in different directions [30] allowing more compression, thereby affecting the shock velocity. Dynamic structural changes are seen during shock propagation along the three directions in our simulations as shown in Fig. 6a, b and c. The figures show the structural changes due to the shock compression at particle velocities 3 km/s and 1 km/s in the X, Y and Z directions. The left half of the figures are the shocked region. Carbon atoms are shown in red, hydrogen in blue, nitrogen in magenta and oxygen in yellow. Note from the figures that we have dynamic conditions wherein the N − NO 2 bonds of the RDX molecule changes orientation dynamically thereby showing α, γ and various transient structural phases. Therefore, the slope of the U s − U p relation will be different in X, Y and Z direction depending upon the orientation of the bonds.
Moreover, due to an extra compression being made possible by re-orientation of molecular bonds, it may be possible that we do not see the change in slope in the X direction whereas in the Y and Z directions, we see a change in the slope of U s − U p . It is also seen that RDX partially decomposes, within the time scale of the MD simulations, at U p greater than 2 km/s for impacts in all the three directions. We describe the decomposition in more detail in the next sub-section.

Results for the RDX decomposition pathways
It is seen from Fig. 7 that decomposition of RDX occurs at particle velocity greater than 2 km/s for impact along the X direction. Similar behaviour of RDX decomposition is observed for impacts along the Y and Z directions. For example, see Fig. 8 for decomposition of RDX for impact along the Y direction. As expected, the rate of RDX decomposition is slower for U p = 2.5km/s than for U p = 3km/s. Earlier ab initio studies by Chakraborty et al. [17,18] have shown that the RDX decomposes along three primary pathways as illustrated in Figs. 9, 10 and 11. They have carried out detailed mechanistic study of thermal decomposition of RDX using ab initio DFT methodology. The optimisation of reactants, products, intermediates and transition states was done using B3LYP flavour of the density functional theory and transition energies were obtained for the various decomposition products. Figure 9 shows the decomposition pathway 1 (DP1) in which there is formation of MN molecules (CH 2 NNO 2 ) due to the concerted ring fission, wherein the C-N bond in the ring structure of RDX breaks.
The second decomposition pathway (DP2) is the homolytic fission of the N-N bond leading to detachment of the NO 2 molecule from an N atom in the ring of RDX, forming a RDR molecule and NO 2 as the first step (Fig. 10). The RDR molecule can subsequently decompose along three different sub-pathways (DP2-SP1, DP2-SP2 and DP2-SP3) as shown in Fig. 10. The third decomposition pathway (DP3) is when HONO gets detached from RDX molecule as shown in Fig. 11. The shock in RDX leads to compression of the molecule, which may bring the H atom and NO 2 molecule close to each other leading to HONO formation. The remaining molecule after HONO detachment is called INT175 molecule on which further HONO detachment occurs forming INT128 molecule. INT128 can then decompose along two subpathways, DP3-SP1 and DP3-SP2, as shown in the figure.
Pahari et al. [19] have carried out kinetic Monte Carlo (KMC) simulations of the decomposition of RDX at different pressures and temperatures and studied the decomposition pathways proposed by Chakraborty [17,18]. They incorporate the pressure and temperature depend based on the Chemkin-based chemical kinetics [31] in the KMC code. The KMC simulations have shown that at low temperatures DP2 mechanism with subsequent decomposition along the sub-pathway SP3 is the dominant decomposition pathway and at higher temperatures (above 4500 K), the concerted ring fission of RDX into three MN molecules (DP1) dominates. From Figs. 12 and 13, for impact along X and Y directions at particle velocity 2.5 km/s, it is seen that RDR and INT175 molecules are formed in the MD simulations. MN molecules do not form implying that decomposition pathways DP2 and DP3 are preferred. This is at odds with the KMC simulations of Pahari et al., where DP2 is the preferred pathway. The reason for this could be that in KMC simulations only the rates of interactions are taken into account, whereas in a MD shock simulation, the explicit coming together of the H and NO 2 atoms occurs leading to formation of the INT175 and HONO molecules as in decomposition pathway DP3.
From Figs. 14 and 15, formation of a few MN molecules is observed in impacts with particle velocity 3 km/s in the X and Y direction. Therefore, at higher impact velocities, the possibility of decomposition pathway DP1 is possible. For 3.0 km/s particle velocity, the temperature in the shocked region is 1238 K, as computed by fitting a Maxwellian distribution to the atomic velocities. Our earlier work [19] shows that mostly DP2 and DP3 dominate for the temperature below 2000 K. Production of few MN molecules by breaking of CN ring is observed in shock simulation at particle velocity of 3 km/s, which means that the contribution from DP1, which dominates at 4500 K, is minor at 1238 K. Figures 16 and 17 show formation of NO 2 molecule during X and Y direction impact at particle velocities 2.5 and 3.0 km/s. Formation of O 2 , INT128 and NO for both X and Y direction impact have also been observed. Table 4 summarises these results.

Conclusions
Molecular dynamics simulations of 1D shock propagation have been carried out in three different crystallographic directions in single crystal RDX using ReaxFF, a reactive force field that allows bond breaking and formation. The U s U p relation for these three directions is obtained and compared with the earlier simulations [5] and experimental results [14]. These MD simulations overestimate the intercept values (by 12%) for X-direction impact as compared with the experimental results. The slope values for X-direction impact shows reasonable match within 3.1% of the reported data. It is seen that the slope and intercept values change for different impact directions. Moreover, for impact along the Y and Z directions, two slopes are   observed above and below 2 km/s particle velocities. It was shown that dynamic structural changes occurred during the shock propagation along the three directions. These dynamic changes were due to rearrangement of orientation of molecular bonds during the impact. We have shown that these rearrangements are different for high and low impact velocities (pressures). These rearrangements will be different for impacts along X, Y and Z direction depending upon the orientation of the bonds which is the reason for different U s − U p behaviour.
The decomposition of RDX during the impact simulations has been observed for U p > 2 km/s. The decomposition products from the MD simulations are compared with three proposed possible decomposition pathways of RDX obtained from ab initio simulations [17,18], the rates of which were used in KMC simulations to obtain the preferred decomposition pathway [19]. It is seen that the MD simulations show similar decomposition trends as predicted by the KMC simulations implying that the ReaxFF potential results from MD simulations matches the ab initio predictions. The formation of RDR − NO 2 and INT175 signifies that RDX decomposition is mainly through decomposition pathways DP2 and DP3. The formation of small amount of MN molecules at particle velocity 3 km/s and for the X and Y direction impact indicates that at higher impact velocities, the possibility of DP1, which is the concerted breaking of