RDX Compression, α→ γ Phase Transition, and Shock Hugoniot

Aug 26, 2016 - Prediction of the density and lattice compression properties of the α and ... the available experimental Hugoniot shock data for both ...
2 downloads 3 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCC

RDX Compression, α→ γ Phase Transition, and Shock Hugoniot Calculations from Density-Functional-Theory-Based Molecular Dynamics Simulations Dan C. Sorescu*,† and Betsy M. Rice‡ †

National Energy Technology Laboratory, U.S. Department of Energy, Pittsburgh, Pennsylvania 15236, United States U.S. Army Research Laboratory, Aberdeen Proving Ground, Maryland 21005, United States

‡

S Supporting Information *

ABSTRACT: Prediction of the density and lattice compression properties of the α and γ phases of the hexahydro-1,3,5-trinitro1,3,5-s-triazine (RDX) crystal and of the low-pressure α → γ phase transition upon pressure increase are general tests used to assess the accuracy of density-functional-theory- (DFT-) based computational methods and to identify the essential parameters that govern the behavior of this high-energy-density material under extreme conditions. The majority of previous DFT studies have analyzed such issues under static optimization conditions by neglecting the corresponding temperature effects. In this study, we extend previous investigations and analyze the performance of dispersion-corrected density functional theory to predict the compression of RDX in the pressure range of 0−9 GPa and the corresponding α → γ phase transition under realistic temperature and pressure conditions. We demonstrate that, by using static dispersion-corrected density functional theory calculations, direct interconversion between the α and γ phases upon compression is not observed. This limitation can be addressed by using isobaric−isothermal molecular dynamic simulations in conjunction with DFT-D2-calculated potentials, an approach that is shown to provide an accurate description of both the crystallographic RDX lattice parameters and the dynamical effects associated with the α→ γ phase transformation. An even more comprehensive and demanding analysis was done by predicting the corresponding shock Hugoniot curve of RDX in the pressure range of 0−9 GPa. It was found that the theoretical results reproduce reasonably well the available experimental Hugoniot shock data for both the α and γ phases. The results obtained demonstrate that a satisfactory prediction of the shock properties in high-energy-density materials undergoing crystallographic and configurational transformations is possible through the combined use of molecular dynamics simulations in the isobaric− isothermal ensemble with dispersion-corrected density functional theory methods. phases denoted as the α and β phases are known to exist, but of these, only the α phase was found to be stable. A full characterization of the crystallographic parameters of this phase was obtained by single-crystal neutron diffraction.9 The molecular crystal (see Figure 1a) has an orthorhombic unit cell with the Pbca space group and Z = 8 molecules per unit cell (24 C, 48 N, 48 O, and 48 H atoms). The conformation of the RDX molecules in the α phase is that in which two of the nitro groups are in the pseudoaxial (A) position and the third nitro group is in the pseudoequatorial (E) position, denoted as the AAE conformation (see the left panel in Figure 1b). With increasing pressure, a phase transition to the so-called γ phase (see the right panel in Figure 1a) was determined in the early studies of Olinger et al.10 This phase was fully characterized

1. INTRODUCTION Predictions of the crystal density and detonation velocity are essential metrics in computer-assisted methods for the development of new formulations of energetic materials. Such properties are directly dependent on the structure and type of interactions present in energetic materials, which, in turn, can undergo significant modifications under the high-temperature and -pressure conditions characteristic of detonation processes.1,2 When accurate methods for evaluating interatomic interactions are available, computational analysis can be very attractive for obtaining a large number of useful properties including the crystallographic packing,3−5 the lattice energy landscape and stability ranges of different polymorphs,6 and the interplay between the mechanical shock loading and the induced chemistry.7,8 One of the most important nitramine systems widely used for military applications is RDX (hexahydro-1,3,5-trinitro-1,3,5s-triazine).2 Under ambient conditions, two polymorphic © XXXX American Chemical Society

Received: June 24, 2016 Revised: August 12, 2016

A

DOI: 10.1021/acs.jpcc.6b06415 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

molecular potentials for nitramines.21,22 Flexible potentials, even if successful in predicting nonreactive dynamics properties, have been found to have their own specific limitations. For example, using a fully flexible molecular potential, Munday et al.23 showed that the α → γ phase transition was observed only for uniaxial deformation along the crystallographic c axis, whereas no such phase transition could be observed under hydrostatic loading. A significant further step in the development of classical force fields included the possibility of describing chemical reactivity along with the structure of the shock fronts. This was achieved by the development of the reactive force field (ReaxFF),24 a methodology that was successfully applied to a large number of relevant problems including the description of the structure of shock fronts in the RDX crystal prior to detonation,25 thermal decomposition,26 ultrafast chemistry under nonequilibrium conditions,8 and energy-transfer dynamics at a hot-spot interface.27 Improving the accuracy of theoretical predictions for the original ReaxFF24 continues to be an area of interest, achieved by the addition of new London dispersion terms28 (which were missing in the original form of the potential) or as demonstrated more recently by parametrization with additional (SAPT)DFT (where SAPT stands for symmetry-adapted perturbation theory) calculations of RDX and 1,1-diamino-2,2-dinitroethene (FOX-7) dimer interactions.29 Aside from classical force fields, the second major approach used to describe RDX theoretically is based on first-principles calculations using, for computational convenience, density functional theory (DFT) methods. Early RDX studies30 employing standard DFT methods demonstrated the existence of large errors for crystal structure predictions with deviations as high as 9.6% between the calculated lattice parameters and the corresponding experimental values. Such errors were attributed to the lack of long-range dispersion interactions, which were missing in standard DFT methods. It was further demonstrated that the agreement with experimental data improves in the case of simulations performed at pressures in excess of 6 GPa, a fact associated with a diminished contribution of the dispersion interactions relative to the repulsive interactions among the molecules.31 Significant improvement in the accuracy of the DFTpredicted crystallographic parameters under ambient conditions for molecular crystals in general and for the RDX crystal in particular was achieved by considering methods that include long-range dispersion interactions. Shimojo et al.32 demonstrated that good agreement with experimental data can be obtained by using the empirical DFT-D2 method originally introduced by Grimme33 with a minimal additional computational cost. In contradistinction, the nonempirical van der Waals density functional (vdW-DF)34 method was also found to provide an accurate description of crystallographic parameters but at the expense of orders-of-magnitude more computation. Using a similar DFT-D2 method, we demonstrated35 that the crystallographic properties of 10 molecular crystals representative of energetic-materials applications, including the α and γ phases of RDX, can be very well predicted under both ambient and compression conditions. More recently, Hunter et al.36 similarly used the DFT-D2 method to predict a variety of properties of the α, β, γ, and ε forms of RDX, including the lattice energies of the α and β forms of RDX, for which their relative values was found to be in excellent agreement with the corresponding experimental data. Similar good agreement was also obtained for the calculated

Figure 1. (a) Pictorial view of the unit cell of RDX crystal in the (left) α and (right) γ phases. (b) Molecular conformations of RDX molecules in the (left) α and (right) γ phases. In the α phase, the crystal has a single molecule in the asymmetric unit with an AAE conformation, whereas in the γ phase, there are two molecules in the asymmetric unit with AAE and AAI conformations. The orientation of the nitro group with respect to the ring plane is given by the angle δ, defined as the angle between the N−N bond and the C−N−C plane11 (as shown in the left panel).

more recently by Pulham and collaborators11 based on neutron powder diffraction and single-crystal X-ray diffraction. Specifically, over the investigated pressure range of 3.90−7.99 GPa, the γ phase was found to remain orthorhombic (Z = 8 molecules per unit cell) with the Pca21 space group and with two independent molecules in the asymmetric unit. The molecular conformations of these molecules were found to be different than that of the α phase and were described as being of the forms AAE and AAI (see the right panel in the Figure 1b), where I represents an intermediate orientation of the NO2 group, midway between the axial and equatorial positions.11 Two other high-pressure phases of RDX (the δ and ε phases) were also observed experimentally,12,13 but we will not discuss them further as they are not the focus of this study. Over the past decade, significant improvements have been achieved in the development of computational methods capable to describe the structural and energetic properties of highenergy-density materials under realistic external conditions. Here, we limit ourselves and discuss briefly such developments only for the case of the RDX crystal considered in this study. Much of the early efforts focused on prediction of the structure and density by ab initio crystal structure predictions using classical force fields.14 In this context, our own investigations15−19 were based on the development of analytic classical force fields fitted to reproduce the structural and energetic properties of molecular crystals under ambient conditions. Such a potential, developed initially for the case of the α-RDX phase,15 was demonstrated to be transferable to a large number of nitramine crystals16−19 and to other non-nitramine compounds from the CHNO class20 common for energeticmaterials applications. A primary limitation of such force fields was the rigid nature of the molecular structure. Such a limitation was addressed subsequently by several groups through the development of nonreactive, fully flexible B

DOI: 10.1021/acs.jpcc.6b06415 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

conserving pseudopotentials of Goedecker, Teter, and Hutter (GTH).46−48 A real-space basis set with a cutoff of 1500 Ry was used for expansion of the electron density. Crystal lattice optimizations were performed by employing a full-dimensional variable-cell optimization at a given imposed external pressure and using an analytic stress tensor. In these cases, the conjugate-gradient (CG) method was used to optimize the cell variables (lattice dimensions and angles), whereas the Broyden−Fletcher−Goldfarb−Shanno (BFGS) optimizer was used for the atomic degrees of freedom. All MD simulations were performed in the constant-temperature and -pressure ensemble using a flexible cell. A Nosé−Hoover thermostat with a time constant of 50 ps was applied to all degrees of freedom, and a barostat time constant of 200 fs was selected. In all simulations, a time step of 0.50 fs was used to integrate the equations of motion Simulations were performed initially for the case of a 1 × 1 × 2 (48 C, 96 N, 96 O, and 96 H atoms) periodic supercell, whereas finite-size effects were determined from simulations employing a 2 × 2 × 2 supercell (192 C, 384 N, 384 O, and 384 H atoms). For pressures larger than 0.5 GPa, all trajectories were integrated for a minimum of 25 ps with averages taken over the last 15 ps. The slowest convergence was observed in the case of trajectories at low pressures (≤0.5 GPa) and ambient temperature. In these cases, trajectories were integrated for at least 40 ps. The lattice parameters were obtained as averages over the last 20 ps of the trajectory. 2.2. Hugoniot Equation of State. The BOMD simulations were also used to determine the Hugoniot equation of state following the procedure originally developed by Erpenbeck.49 In this case, the Hugoniot points satisfy the relation

phonon modes relative to the experimental inelastic neutron scattering spectrum. We note that currently a large number of computational approaches exist that attempt to correct standard DFT methods for dispersion interactions, and a general classification of these methods can be found in the review of Klimeš and Michaelides.37 In our own work,38 we have compared the performances of more than 17 different methods with corrections for long-range dispersion interactions to describe the crystallographic properties of a set of 26 high-nitrogencontent salts relevant for energetic-materials applications. These methods ranged from adding atom−atom dispersion corrections with environment-independent and environmentdependent coefficients to methods that incorporate dispersion effects through dispersion-corrected atom-centered potentials (DCACPs),39 to methods that include nonlocal corrections.34,40,41 For the specific case of the RDX crystal, some of these methods have also been applied.32,42 For example, Balu et al.42 obtained very good agreement with experiment when using DFT methods and dispersion-corrected atom-centered pseudopotentials developed by Rothlisberger and colleagues.39 A common characteristic of previous DFT studies32,35,36,42 related to simulations of RDX is the fact that the predicted lattice parameters were obtained using static optimizations under the assumption that the thermal effects do not make a large contribution. However, the accuracy of such an assumption has not yet been tested, particularly if one recalls that the compression of the RDX crystal leads not only to variations in the lattice dimensions but also to changes in crystal symmetry and molecular conformations. As shown in this work, the lack of thermal effects, as in the case of static calculations, can have a negative impact even on simple prediction of the α → γ phase transition. In this work, we address these limitations by considering computational methods that include both thermal and compression effects. Such an approach allowed us to fully observe the α → γ phase transition and to corroborate the intramolecular conformational changes associated with this solid−solid phase transition. Furthermore, we extend the analysis to include the prediction of shock Hugoniot curves in the pressure range of 0−9 GPa where experimental data for both the α and γ phases are available for comparison.10 The organization of the article is as follows: In section 2, we discuss the computational methods employed in this work. The results of static and dynamic simulations are analyzed in section 3, as are the calculated Hugoniot data. Further comparison of the theoretical predictions with the corresponding experimental data is provided in the same section. The main conclusions of this work are summarized in section 4.

U − U0 =

1 (P + P0)(V − V0) 2

(1)

where U is the internal energy, P is the pressure, V is the volume, and the subscript 0 indicates the unshocked reference state. The points satisfying eq 1 can be obtained following the iterative procedure of Brennan and Rice.50 Briefly, for each pressure, a set of different trajectories were run at different temperatures that produce a set of equation-of-state points. An interpolation fit can then be used to determine the Hugoniot temperature at which eq 1 is satisfied. Essential in the selection of the set of temperatures is that they bracket the corresponding Hugoniot temperature. The calculated Hugoniot data obtained are transformed into a format typically used to report experimental shock Hugoniot information, in terms of the shock velocity (Us)

2. COMPUTATIONAL METHODS 2.1. Crystal Lattice Optimizations and Born−Oppenheimer Molecular Dynamics Simulations. The optimizations of crystal configurations of RDX in different phases and the dynamic properties of RDX under external temperature and pressure conditions were investigated using Born−Oppenheimer MD (BOMD) simulations as implemented in the Quickstep module43 of the CP2K program.44 The calculations used the Perdew−Burke−Ernzerhof (PBE) functional45 with the Grimme-D233 correction for the long-range dispersion, hereafter denoted as PBE-D2. The Gaussian triple-ξ TZVP valence basis set augmented with one set of d-type and p-type polarization functions was used in conjunction with the norm-

Us =

ρ (P − P0) ρ0 (ρ − ρ0 )

(2)

and the particle velocity (Up) Up =

(P − P0)

(ρ − ρ0 ) ρρ0

(3)

where ρ and ρ0 are the system densities at finite pressure P and in the reference state at P0, respectively. Rearrangement of eqs 2 and 3 provides the P−V information corresponding to eq 1 C

DOI: 10.1021/acs.jpcc.6b06415 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C V=

V0(Us − Up) Us

(4)

and

P = P0 +

UU s p V0

(5)

3. RESULTS AND DISCUSSIONS 3.1. Static Crystal Optimizations at Ambient Pressure. It was previously documented35 that optimization of the crystal lattice based on the calculated stress tensor requires cutoff energies of the plane-wave basis set larger than those generally used to obtain the convergence of energies. Similarly, the first task in the current study was to determine an optimal cutoff energy for the real-space integration grid that ensured convergence of the stress tensors and implicit of the predicted lattice parameters under different compression conditions. For this purpose, a full-dimensional variable-cell optimization at ambient pressure was performed using an analytical stress tensor in which both the ionic degrees of freedom and the unit cell parameters were optimized. Calculations were performed for the case of the 1 × 1 × 2 and 2 × 2 × 2 supercells. We note that an increase in supercell size was required, as our version of the CP2K program was a gamma-point-only code. The corresponding results illustrating the variations of the predicted lattice dimensions and unit cell volume with cutoff energy are presented only for the case of the 2 × 2 × 2 supercell in the Supporting Information (Figure S1). As shown in this figure, convergence was achieved for cutoff energies of 1500 Ry. For this cutoff energy, the predicted lattice constants of the α-RDX crystals are a = 13.252 Å, b = 11.440 Å, and c = 10.703 Å. These values are in extremely close agreement with the neutron diffraction values obtained by Choi and Prince,9 with a maximum error of −1.16% for the case of the b lattice constant. The predicted unit cell volume was only −0.68% different from the corresponding experimental value. We note that a similar level of agreement was obtained when using the smaller 1 × 1 × 2 supercell, indicating that the lack of k-point sampling in the Brillouin zone did not make a major contribution for the sizes of the simulation boxes used here. 3.2. Static Compression of RDX. Having identified the optimal cutoff energy for the real-space integration grid and demonstrated the accuracy level of the static calculations for the ambient phase of α-RDX, we focused on the prediction of the lattice compression using lattice minimization calculations. These were done for both the α and γ phases of RDX. As demonstrated in previous experimental investigations,11 the α phase is stable in the pressure range of 0−3.8 GPa, whereas the γ phase is stable above 3.9 GPa. Similarly to the calculations performed under ambient conditions, the current compression studies at different pressures were performed by relaxing both the simulation supercell and all of the atoms in this cell by minimizing the forces and the stress tensor subject to the external imposed pressure. The variations in the lattice parameters and unit cell volume with pressure predicted by static optimizations are provided in Figures 2a and 3a, respectively. In the same figures, we include for comparison the experimental values obtained by Olinger et al.10 (open white symbols) and by Hunter et al.36 (solid red symbols) for the α-RDX phase and the data of Davidson et al.11 for the γRDX phase (solid yellow symbols). Also, two sets of calculation

Figure 2. Variations with hydrostatic pressure of the lattice parameters for α-RDX and γ-RDX as predicted by (a) lattice minimization (LM) and (b) BOMD (MD) calculations. For comparison, the experimental data for α-RDX from ref 10 (exp1) and ref 36 (exp2) and for γ-RDX from 11 (exp3) are also included. For the calculated data, the sets denoted with indices I and II correspond to crystal structures whose initial configurations are in the α and γ-RDX phases, respectively.

data are shown in these figures, denoted by I and II. These represent the results for the cases when the initial structure at the start of optimization corresponds to the α phase (set I) and to the γ-RDX phase (set II). For the first set (I), the analysis was done in the pressure range of 0−9 GPa, which obviously extends beyond the stability range observed experimentally for α-RDX. The results from Figures 2a and 3a show two main characteristics. First, within the range of stability observed experimentally for the α (0−3.85 GPa) and γ (3.9−8.0 GPa) phases, the predicted and experimentally measured crystallographic data are in very good agreement. Similar good agreement is observed for the predicted unit cell volumes within the individual ranges of stability of the two phases, such that the experimental and theoretical data are almost superimposed. In the pressure range of 0−3.85 GPa, the maximum percentage deviations of the calculated lattice parameters from experimental data are −1.02% with respect to the Olinger et al.10 data and 1.07% relative to the Hunter et al.36 data set, whereas in the pressure range of 3.90−7.99 GPa for the γ-RDX phase, the maximum error deviation of the predicted lattice parameters is 1.08% relative to the data obtained by Davidson et al.11 The second important characteristic observed is that the calculated lattice parameters and unit cell volumes for set I are in poor agreement with experiments outside the range of stability of phase α. Specifically, from the data in Figures 2a and 3a is it seen that for pressures outside the stability region of the α-RDX phase (i.e., P > 3.9 GPa), simulations started in the α phase evolve on a different curve than those started in the γ phase. In particular, the evolutions of D

DOI: 10.1021/acs.jpcc.6b06415 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 4. Variations in time (represented by gray lines) of the (a) unit cell volume, (b) temperature, (c) pressure, and (d−f) unit cell lattice parameters during a BOMD trajectory in the NPT ensemble for αRDX at ambient pressure and T = 300 K. In each case, the red line represents the running average of the respective quantity.

Figure 3. Variations with hydrostatic pressure of the unit cell volume for α-RDX and γ-RDX as predicted by (a) lattice minimization (LM) and (b) BOMD (MD) calculations. For comparison, the experimental data for α-RDX from ref 10 (exp1) and ref 36 (exp2) and for γ-RDX from 11 (exp3) are also included. For the calculated data, the sets denoted with indices I and II correspond to crystal structures whose initial configurations are in the α and γ-RDX phases, respectively.

imposed external conditions, namely, T = 300 K, P = 1 bar. During this period, there is an important and expected initial expansion of the lattice parameters and unit cell volume, after which there are only small variations in the average crystallographic parameters. The average lattice parameters (a, b, c) extracted from the last 20 ps of this trajectory are 13.398, 11.708, and 10.935 Å, respectively, and the corresponding average unit cell volume is 1715.3 Å. Inspection of these values shows that the average volume is practically identical to the product of the individual average lattice parameters. This result demonstrates that, during this trajectory, the simulation box remains practically orthogonal with very small deviations from 90°. When compared to the experimental neutron diffraction data of Choi and Prince,9 the percentage deviations of the calculated average lattice parameters are found to be relatively small, with values of 1.64%, 1.16%, and 2.12% for the unit cell lengths and 4.98% for the unit cell volume. Overall, inclusion of thermal effects at T = 300 K and ambient pressure leads to a slight expansion of the lattice relative to the configuration obtained in simple lattice minimization calculations. As indicated above, the reliability of the converged averages obtained in this initial trajectory was further tested with the aid of a second trajectory performed at the same temperature and pressure and started from the final configuration of the first trajectory. This new trajectory was run for an additional 35 ps and was fully independent of the first trajectory in terms of the initial set of velocities and calculated averages. As seen in Figure S2, the variation of the calculated lattice averages is very small (except for the initial equilibration period of a few picoseconds) and suggests the presence of a much better equilibrated structure. For this trajectory, the final average lattice dimensions are 13.405, 11.707, and 10.931 Å with a unit cell volume of 1715.4 Å3. These values are very close to those determined from the first trajectory and confirm that the

the lattice dimensions (Figure 2a) and unit cell volumes (Figure 3a) for set I in the pressure range of 3.9−9.0 GPa are clearly distinct from the experimental data set for the γ-RDX phase. These findings point to a major limitation of static optimizations, namely, that the α → γ transition is not obtained upon pressure increase, even for pressures as high as 9 GPa. These results indicate the necessity of using a different approach that includes thermal contributions in addition to compression effects. Such a requirement can be achieved by using isobaric−isothermal molecular dynamic simulations, as described in the next section. 3.3. BOMD Compression Simulations of RDX. In this section, we analyze the compression of the RDX crystal at finite temperatures based on the predictions obtained using molecular dynamics simulations in the isothermal−isobaric ensemble with cells that are flexible in both size and shape. All results presented correspond to the 2 × 2 × 2 supercell unless otherwise indicated. It is instructive to start the theoretical analysis for simulations performed under ambient conditions. The variations in time of the relevant parameters, namely, the unit cell volume, temperature, pressure, and lattice dimensions, are illustrated for two different trajectories in Figure 4 and Figure S2. The first trajectory (shown in Figure 4) was obtained starting from the crystal geometry determined in previous lattice minimizations, whereas the second trajectory (shown in Figure S2) was started from the final configuration in the first trajectory. The purpose of this second trajectory was to further evaluate the convergence of the calculated averages determined in the first trajectory. As seen in Figure 4, both the average temperature and pressure equilibrate in less than 15 ps to the E

DOI: 10.1021/acs.jpcc.6b06415 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C average structures obtained in the two trajectories are similar, independent of the initial starting point of these trajectories. Using a methodology similar to that described above for ambient conditions, we further determined the variations with pressure at ambient temperature of the lattice parameters of RDX for both the α and γ phases. The corresponding results are illustrated in Figures 2b and 3b. In these figures, data sets I and II correspond to structures initially taken to be in the α and γ phases, respectively. From the data in Figure 2b, it is clearly seen that the predicted lattice parameters closely follow the corresponding experimental results. A more precise evaluation of the errors for the predicted lattice parameters and of the pressure dependence of these errors can be obtained from Figure S3. In this case, we used only those pressure results for which experimental data are available. As seen in this figure, over the entire pressure range, the largest error of 2.12% corresponds to predictions under ambient conditions. With increasing pressure, the error level drops below 1.5%. Similarly, the largest error for the unit cell volume of 4.98% was observed for simulations under ambient conditions. With increasing pressure, this error level decreases quickly below 3% for the entire pressure range of 1.75−7.99 GPa. The second major feature of the data given in Figures 2b and 3b is seen in the transition region between the α and γ phases. Specifically, trajectories from set I initiated using crystal configurations in the α phase undergo significant crystallographic changes such that, at pressures above 3.85 GPa, the crystal structure is transformed into the γ phase. More importantly, the equilibrium states reached in these trajectories are very close to the states obtained for set II, which were initiated in the γ phase. These findings are in sharp contrast to those illustrated in panels a of Figures 2 and 3, where the lattice parameters (especially the b and c parameters) for set I do not approach those for set II even for pressures as high as 9 GPa. These results clearly illustrate the important role played by thermal contributions in dynamic simulations in the observation of the pressure-induced phase transition of RDX crystal between the α and γ phases. A more detailed view of the evolution of the lattice parameters during the α → γ phase transition is captured in Figure 5. In this case, the trajectory represents the results obtained at P = 3.9 GPa and T = 300 K, conditions under which the transition from α- to γ-RDX was also observed experimentally.11 In the initial state, the crystal was selected to be in the α phase, having the structure determined from a simple lattice minimization at P = 3.9 GPa. As seen from the time evolution of the lattice parameters, there are important changes in the first 5 ps of the trajectory. The major changes take place for the b and c lattice constants, which exhibit opposite variations, whereas a smaller increase is seen for the a lattice constant. Overall, there is a substantial decrease of the unit cell volume from about 1420 to 1390 Å3, after which the volume stabilizes at an average value of about 1392 Å3. We note that the drop in the unit cell volume at the phase transition was also observed experimentally in the studies by both Olinger et al.10 and Davidson et al.11 One can also raise the question of whether selection of the initial configuration has not already biased the observed phase transition. To address this issue, we considered two other starting configurations: one corresponding to optimized structure at the ambient state of α-RDX and a second corresponding to a compressed state as determined from a MD run at P = 3.36 GPa. Starting from these two different

Figure 5. Variations in time of the (a) unit cell volume; (b) temperature; (c) pressure; and (d−f) lattice dimensions a, b, and c during a BOMD simulation at ambient temperature and P = 3.9 GPa.

configurations, new trajectories at P = 3.90 GPa and T = 300 K were run, and the corresponding results are shown in Figures S4 and S5. As can be seen, regardless of the initial starting crystal state, the system undergoes a phase transition with characteristic crystallographic changes similar to those represented in Figure 5, that is, with an opposing expansion and contraction for the b and c lattice dimensions, respectively, and a substantial drop of the unit cell volume. In all cases, the final average volume is close to the 1392 Å3 value also determined from simulations presented in Figure 5. A similar good agreement is seen for the average lattice parameters, with values very close to those obtained in Figure 5 (a = 12.877 Å, b = 11.117 Å, c = 9.724 Å). These results clearly confirm the fact that the observed phase transition is not dependent on the initial starting configuration as long as it is a state belonging to the α phase. 3.4. Intramolecular Structural Analysis. We further tested the accuracy of the current simulations not only to predict the variation of the unit cell with pressure but also to determine the pressure-induced changes of the molecular conformations. For RDX, the most prominent conformational change with increasing pressure corresponds to the specific orientation of the nitro groups relative to the ring of the molecule and can be depicted by the angle δ between the C− N−C ring atoms and the corresponding N−N bond (defined in Figure 1b11). As described in the Introduction, in the α phase, all molecules of the crystal have the same AAE conformation, whereas in the γ phase, there are two groups of molecules with different conformations denoted as AAE and AAI. Consequently, the analysis of the δ angles for each of the three nitro groups of the RDX molecule can be used as a convenient tool for identifying the intramolecular conformation and, hence, the conformational phase of the system. Moreover, as shown in this section, monitoring of the δ angles can also be used to corroborate the transition between RDX crystal phases. F

DOI: 10.1021/acs.jpcc.6b06415 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Time-averaged values for the set of three δ angles for each molecule in the unit cell were also averaged over all unit cells that make up the supercell at different pressures and ambient temperature. The values of the set of three δ angles were then averaged over the two types of molecules within the γ phase, and the results are shown in Figure 6. These molecules were

close to the experimentally observed data. Under ambient conditions, only two of the predicted values are very close to the experimental data, whereas a larger difference of about 7° is seen for the third angle. The source of this difference is not clear at this time. Further confirmation of the good agreement with experiment for the predicted conformations of RDX was obtained from an analysis of the individual atomic coordinates recorded during a trajectory at P = 5.2 GPa and T = 300 K for the case of the γRDX crystal. Specifically, we calculated the temporally and spatially averaged crystal and molecular structures from BOMD simulations as previously described29 and graphically superimposed them on the corresponding experimental structure. The results of this comparison are illustrated in Figure 7,

Figure 6. Variations with pressure of the δ angles describing the orientations of the three nitro groups of the RDX molecule (as defined in Figure 1b). The left panel corresponds to a molecule that maintains the AAE conformation in both α and γ phases, whereas the right panel corresponds to a molecule that has the AAE conformation in the α phase and maintains the AAI conformation in the γ phase.

chosen based on their behavior over the range of pressures studied. The first type is representative of the molecules that maintain the AAE conformation in both the α and γ phases (Figure 6a). The other is representative of the second type of molecule in the γ phase, which transitions from AAE in the α phase to AAI (Figure 6b). In the same figure, we indicate the corresponding angles obtained for experimental structures in the α phase under ambient conditions and for the γ phase at P = 5.2 GPa.9,11 Several general characteristics can be noticed from the data in Figure 6. First, there is a clear distinction between the angles calculated at pressures below 3.9 GPa and above this value. For pressures below 3.9 GPa, the system is characterized by only three δ angles, whereas above this pressure, there are two distinct groups of three δ values, fully consistent with the experimental findings. A second important result is that the variation of these angles shows a very sharp discontinuity at pressures close to 3.9 GPa where the experimental α → γ transition takes place. This confirms our initial hypothesis that monitoring the variation of the δ angles can be used as a tool for corroborating the phase transition. It should also be noted that the pressure range at which the transition from the set of three independent δ angles to the set of six independent angles takes place, namely, around 3.9 GPa, is also consistent with the pressure at which sharp variations in the cell dimensions were observed (see Figures 2b and 3b). It is interesting to observe that, in the γ phase, five of the six δ angles have relatively modest changes in the pressure range of 4−9 GPa. In only one case (represented with solid red symbols in the right panel of Figure 6) is there a pronounced and continuous increase with pressure. In terms of the accuracy of the predicted values, it can be seen that this is indeed the case, particularly for the γ phase, where the calculated values are very

Figure 7. (a) Pictorial view of the temporally and spatially averaged unit cell of the RDX crystal in the γ phase as determined using BOMD simulations at P = 5.2 GPa and T = 300 K superimposed on the corresponding experimental structure from ref 11. (b,c) Pictorial view of the averaged molecular structures of two of the molecules in the unit cell described in panel a. These represent the two independent molecules in the asymmetric unit superimposed on the corresponding experimental structures.

together with individual details for molecules with the AAE and AAI conformations (shown in the bottom two panels of this figure). From both the general crystal view (Figure 7a) and the individual molecular structures (Figure 7b,c), it is clearly seen that the two superimposed structures (i.e., the calculated and experimental ones) are very similar and that there is little difference in the corresponding molecular orientations. These results further confirm the very good agreement of the theoretical predictions with the experimental data. G

DOI: 10.1021/acs.jpcc.6b06415 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C The ensemble of results in this and previous sections demonstrate that BOMD simulations using the PBE-D2 method are able to provide an overall accurate description of the lattice parameters and unit cell volume of the RDX crystal in the 0−9 GPa pressure interval, including the α- → γ-RDX phase transition. The calculated phase transition is seen theoretically to take place at 3.85 GPa, a pressure value very close to the experimentally observed value of 3.9 GPa. Furthermore, the calculated conformations of the individual molecules of the crystal as described by the evolution of the δ angles (see Figure 1) were found to be fully consistent with the corresponding experimental data. Specifically, simulations correctly identified the number of groups of molecules with different conformations, namely, one in the case of the α phase and two in the case of the γ phase. Furthermore, the range of pressures in which the transition from one group of molecules with one conformation to a second group of molecules with different conformations was found to be consistent with the pressure range in which a sharp transition in the unit cell variation was observed and with the experimentally determined transition pressure value. Before we close this section, it is instructive to comment on the system size dependence of the results. As indicated in the Computational Methods section, calculations were also done using a smaller supercell box containing 1 × 1 × 2 cells. This supercell selection did not seem to have a significant impact on the predicted cell values in lattice minimization calculations or in BOMD simulations relative to the values obtained using the larger 2 × 2 × 2 supercell. The major impact, however, is observed in prediction of the α → γ transition pressure, as shown in Figure S6. As seen in this figure, the transition takes place at a lower pressure of 2.75 GPa. The occurrence of this transition at a lower pressure might be due to much larger pressure fluctuations taking place during NPT-MD simulations for the 1 × 1 × 2 supercell than for the 2 × 2 × 2 supercell. These can be compared by inspecting the pressure variations for the two simulation boxes during trajectories under the same conditions, for example, at P = 2.75 GPa and T = 300 K, as shown in Figure S7. As can be seen in this figure, because of the larger pressure fluctuations taking place for the 1 × 1 × 2 supercell (with a maximum deviations of 1.76 GPa from the average value), the instantaneous pressure can reach values significantly in excess of 3.9 GPa, which, in turn, can push the system into an earlier phase transition than the one expected experimentally. This finding clearly indicates that selection of simulation boxes that are too small can have a negative impact in accurate predictions of the solid−solid phase transitions in molecular crystals. This observation also has implications for requirements imposed when accurate prediction of the detonation velocity is of interest, as shown in the next section. 3.5. Hugoniot Equation of State. The performances of BOMD simulations in predicting the properties of RDX under compression under realistic temperature conditions were subjected to an even more demanding test related to the prediction of the shock Hugoniot of the material under compression. As indicated in the Introduction, this involves identification of the Hugoniot states that satisfy eq 1 for each external pressure, from which the corresponding shock (Us) and particle velocity (Up) quantities (see eqs 2 and 3) can be extracted. The calculated Hugoniot results for RDX in the pressure range of 0−9 GPa are provided in Figure 8a, together with the experimental results obtained by Olinger et al.10 and Ilyukhin et al.51 We note that Ilyukhin et al.51 reported shock

Figure 8. Variation of the shock velocity (Us) as a function of the particle velocity (Up) for RDX. (a) Results obtained from simulations using the 2 × 2 × 2 supercells for both the α and γ phases compared to the experimental isothermal static compression data (variables Ust and Upt) and shock Hugoniot data (variables Us and Up) of Olinger et al.10 (exp1), as well as to the shock Hugoniot data of Ilyukhin et al.51 (exp2). The indices I−IV in the legend refer to the following equations: (I) Ust = 2.68 + 1.9Upt, (II) Ust = 2.49 + 1.8Upt, (III) Us = 2.78 + 1.9Up, and (IV) Us = 2.87 + 1.61Up. (b) Results from simulations obtained using the 1 × 1 × 2 and 2 × 2 × 2 supercells showing an early pressure transition from the α to the γ phase in the case of the smaller supercell.

Hugoniot data for RDX in the pressure range of 6.7−15.5 GPa, whereas Olinger et al.10 reported both isothermal volume compression data in a form similar to the shock Hugoniot equation (denoted using variables Ust and Upt), as well as a predicted shock Hugoniot for the α-RDX phase (indicated using variables Us and Up) obtained using thermal expansion and heat capacity values. Simulations in this case were done for the 2 × 2 × 2 supercell and contain two sets of data. These are indicated with round yellow symbols and solid black triangles and represent the results in which the system was taken to be in the α- and γ-RDX phases, respectively, in the initial state. These sets of results are denoted by ph(I,calc.) and ph(II,calc.), where indices I and II refer to the initial phase from which the simulation was started. During the trajectory, the system might remain in the same phase, as is the case at lower pressures, or the system might transition from the α phase to the γ phase. The results in Figure 8a indicate that the observed variation is practically linear with a corresponding inflection region at Up = 0.6 km/s when the transition from one phase to another is observed. Furthermore, it can be seen that the last data point from the ph(I,calc.) set (at Up = 0.64 km/s and corresponding to the pressure P = 3.90 GPa) practically coincides with the first data point from the ph(II,calc.) set. This demonstrates the independence of the predicted shock velocity data from the initial state and phase used in these calculations. The second H

DOI: 10.1021/acs.jpcc.6b06415 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

This is a major limitation of simple optimization methods and points to the need to incorporate thermal effects into simulations. (b) In contradistinction to simple lattice optimizations, it was demonstrated that, by using BOMD simulations in the isothermal−isobaric ensemble, it is possible to describe accurately the variations with pressure of the lattice parameters under ambient conditions and to predict accurately the corresponding α → γ phase transition pressure at about 3.85 GPa, a value that agrees well with the experimentally determined value (3.9 GPa).11 (c) Modification of the intramolecular conformation of RDX in the α and γ phases during the trajectory can be observed by monitoring the δ angles, which describe the orientation of the nitro groups of the molecule relative to its ring. Dynamic evolution of these angles was found to be a quite sensitive indicator describing the type of crystal phase and the transition pressure between phases, in good agreement with available experimental results.9,11 (d) A similar very good agreement with experimental data11 was found for the calculation of the Hugoniot curve and the analysis of the relationship between the shock velocity and particle velocity. The use of larger 2 × 2 × 2 supercells was demonstrated to be necessary for an accurate prediction of the transition pressure, whereas the results obtained using a smaller 1 × 1 × 2 supercell were significantly lower than the measured result. Taken in total, the ensemble of results obtained in the current study based on BOMD simulations with PBE-D2 potentials points to the availability of an accurate computational capability for describing not only the long-range dispersion interactions existing between the molecules in molecular crystals but also the dynamical crystallographic changes related to the solid−solid α- → γ-RDX phase transition, as well as the ample conformational changes taking place in these phases.

main observation is that the calculated results are in reasonably close agreement with the experimental findings of Olinger et al.,10 whereas somewhat larger differences are observed for the γ-RDX phase relative to the older data set obtained by Ilyukhin et al.51 A linear fit of the shock velocities and particle velocities for α-RDX gives a bulk sound velocity of 2.43 km/s and a slope of 2.0. These data compare very well with the corresponding values of 2.68 km/s and 1.9 obtained in ref 10. Similarly, for the γ-RDX phase, a linear fit of the form Us = u0 + sUp leads to an intercept of u0 = 2.28 km/s and a slope of s = 1.90, which also compare very well to the values of 2.49 km/s and 1.8 determined by Olinger et al.10 Both of these sets of results reinforce the accuracy of the computational method used. The variation of the calculated P−T Hugoniot data determined in this study was also compared to the corresponding experimental phase diagram of RDX in Figure S8. As seen in this figure, for the range of pressures and temperatures analyzed, the Hugoniot path crosses from the α phase to the γ phase but not to the ε phase at high pressures. As in the prediction of the lattice parameters, we observed a notable dependency of the Hugoniot data on the system size. To illustrate this dependency, in Figure 8b, we provide a comparison of the Hugoniot data obtained for the 1 × 1 × 2 and 2 × 2 × 2 supercells. As in the previous calculations, two sets of points, I and II, were considered having the same connotation as described above. The data for the 1 × 1 × 2 and 2 × 2 × 2 supercells were found to have the same linear dependence and practically overlapping values both at low pressures (P < 2.75 GPa, Up < 0.33 km/s) for the α-RDX phase and at higher pressures (P > 3.9 GPa, Up > 0.61 GPa) for the γRDX phase. The main differences were observed in the intermediate regime of pressures, 2.75 < P < 3.9 GPa, or alternatively for particle velocities in the range 0.33 < Up < 0.61 km/s. As mentioned in the previous section, this is due to a deficiency in the simulations performed using the smaller 1 × 1 × 2 simulation box, in that it predicts a phase transition at a lower pressure value (P = 2.75 GPa) than observed experimentally (P = 3.9 GPa). In turn, simulations for the 2 × 2 × 2 supercell identify the phase transition at a pressure much closer to the experimental value of 3.85 GPa. These findings suggest that selection of supercell sizes that are too small, preferred in DFT simulations because of their significantly lower computational costs, should be carefully tested to ensure that no detrimental effect on the accuracy of the corresponding predictions will be incurred, particularly when related to the dynamic effects associated with structural phase transitions and molecular reorganization, as is the case in the present study.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b06415. Additional results related to the dependence of the lattice parameters on the cutoff energy, lattice parameters predicted from BOMD simulations under ambient conditions and at specific pressures, percentage errors of the predicted lattice parameters and unit-cell volumes at different pressures and ambient temperature, and the P−T Hugoniot data overimposed on the phase diagram of RDX (PDF)

4. CONCLUSIONS In this work, we investigated the compression, the low-pressure α → γ phase transition, and the shock Hugoniot properties of RDX in the pressure range of 0−9 GPa using BOMD simulations in the isotheral-isobaric ensemble and complementary static optimizations using the dispersion-corrected PBE-D2 method. The main findings of this study can be summarized as follows: (a) Simple static optimizations provide an accurate description of the dependence on pressure of the lattice parameters and unit cell volume for both the α- and γ-RDX phases. However, the transition between these two phases observed experimentally at about 3.9 GPa is not predicted in simple lattice optimizations even for pressures as high as 9 GPa.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: 412-386-4827. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS All calculations were performed at various DOD High Performance Computing Program Major Shared Resource Center Sites. I

DOI: 10.1021/acs.jpcc.6b06415 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

■

Nitramine Intermolecular Potential to Non-Nitramine Crystals. J. Phys. Chem. A 1999, 103, 989−998. (21) Smith, G. D.; Bharadwaj, R. K. Quantum Chemistry Based Force Field for Simulations of HMX. J. Phys. Chem. B 1999, 103, 3570−3575. (22) Cawkwell, M. J.; Sewell, T. D.; Zheng, L. Q.; Thompson, D. L. Shock-Induced Shear Bands in an Energetic Molecular Crystal: Application of Shock-Front Absorbing Boundary Conditions to Molecular Dynamics Simulations. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 014107. (23) Munday, L. B.; Chung, P. W.; Rice, B. M.; Solares, S. D. Simulations of High-Pressure Phases in RDX. J. Phys. Chem. B 2011, 115, 4378−4386. (24) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (25) Nomura, K.-i.; Kalia, R. K.; Nakano, A.; Vashishta, P.; van Duin, A. C. T.; Goddard, W. A. Dynamic Transition in the Structure of an Energetic Crystal during Chemical Reactions at Shock Front Prior to Detonation. Phys. Rev. Lett. 2007, 99, 148303. (26) Strachan, A.; Kober, E. M.; van Duin, A. C. T.; Oxgaard, J.; Goddard, W. A. Thermal Decomposition of RDX from Reactive Molecular Dynamics. J. Chem. Phys. 2005, 122, 054502. (27) 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. (28) Liu, L. C.; Liu, Y.; Zybin, S. V.; Sun, H.; Goddard, W. A. 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. (29) Larentzos, J. P.; Rice, B. M.; Byrd, E. F. C.; Weingarten, N. S.; Lill, J. V. Parameterizing Complex Reactive Force Fields Using Multiple Objective Evolutionary Strategies (MOES). Part 1: ReaxFF Models for Cyclotrimethylene Trinitramine (RDX) and 1,1-Diamino2,2-Dinitroethene (FOX-7). J. Chem. Theory Comput. 2015, 11, 381− 391. (30) Byrd, E. F. C.; Scuseria, G. E.; Chabalowski, C. F. An Ab Initio Study of Solid Nitromethane, HMX, RDX, and CL20: Successes and Failures of DFT. J. Phys. Chem. B 2004, 108, 13100−13106. (31) Byrd, E. F. C.; Rice, B. M. Ab Initio Study of Compressed 1,3,5,7-Tetranitro-1,3,5,7-Tetraazacyclooctane (HMX), Cyclotrimethylenetrinitramine (RDX), 2,4,6,8,10,12-Hexanitrohexaaza-isowurzitane (CL-20), 2,4,6-Trinitro-1,3,5-Benzenetriamine (TATB), and Pentaerythritol Tetranitrate (PETN). J. Phys. Chem. C 2007, 111, 2787−2796. (32) Shimojo, F.; Wu, Z. Q.; Nakano, A.; Kalia, R. K.; Vashishta, P. Density Functional Study of 1,3,5-Trinitro-1,3,5-triazine Molecular Crystal with van der Waals Interactions. J. Chem. Phys. 2010, 132, 094106. (33) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (34) Dion, M.; Rydberg, H.; Schroder, E.; Langreth, D. C.; Lundqvist, B. I. Van der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92, 246401. (35) Sorescu, D. C.; Rice, B. M. Theoretical Predictions of Energetic Molecular Crystals at Ambient and Hydrostatic Compression Conditions Using Dispersion Corrections to Conventional Density Functionals (DFT-D). J. Phys. Chem. C 2010, 114, 6734−6748. (36) Hunter, S.; Sutinen, T.; Parker, S. F.; Morrison, C. A.; Williamson, D. M.; Thompson, S.; Gould, P. J.; Pulham, C. R. Experimental and DFT-D Studies of the Molecular Organic Energetic Material RDX. J. Phys. Chem. C 2013, 117, 8062−8071. (37) Klimeš, J.; Michaelides, A. Perspective: Advances and Challenges in Treating van Der Waals Dispersion Forces in Density Functional Theory. J. Chem. Phys. 2012, 137, 120901. (38) Sorescu, D. C.; Byrd, E. F. C.; Rice, B. M.; Jordan, K. D. Assessing the Performances of Dispersion-Corrected Density Functional Methods for Predicting the Crystallographic Properties of High

REFERENCES

(1) Wescott, B. L.; Stewart, D. S.; Davis, W. C. Equation of State and Reaction Rate for Condensed-Phase Explosives. J. Appl. Phys. 2005, 98, 053514. (2) Akhavan, J. The Chemistry of Explosives, 3rd ed.; Royal Society of Chemistry: Cambridge, U.K., 2011. (3) Pertsin, A. J.; Kitaigorodsky, A. I. The Atom−Atom Potential Method: Applications to Organic Molecular Solids; Springer Series in Chemical Physics; Springer-Verlag: Berlin, 1987. (4) Mendoza-Cortes, J. L.; An, Q.; Goddard, W. A.; Ye, C. C.; Zybin, S. Prediction of the Crystal Packing of Di-Tetrazine-Tetroxide (DTTO) Energetic Material. J. Comput. Chem. 2016, 37, 163−167. (5) Politzer, P.; Murray, J. S. Perspectives on the Crystal Densities and Packing Coefficients of Explosive Compounds. Struct. Chem. 2016, 27, 401−408. (6) Price, S. L. Predicting Crystal Structures of Organic Compounds. Chem. Soc. Rev. 2014, 43, 2098−2111. (7) Zhou, T. T.; Zybin, S. V.; Liu, Y.; Huang, F. L.; Goddard, W. A. Anisotropic Shock Sensitivity for β-Octahydro-1,3,5,7-Tetranitro1,3,5,7-Tetrazocine Energetic Material under Compressive-Shear Loading from ReaxFF-lg Reactive Dynamics Simulations. J. Appl. Phys. 2012, 111, 124904. (8) Wood, M. A.; Cherukara, M. J.; Kober, E. M.; Strachan, A. Ultrafast Chemistry under Nonequilibrium Conditions and the Shock to Deflagration Transition at the Nanoscale. J. Phys. Chem. C 2015, 119, 22008−22015. (9) Choi, C. S.; Prince, E. The Crystal Structure of Cyclotrimethylenetrinitramine. Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1972, 28, 2857−2862. (10) Olinger, B.; Roof, B.; Cady, H. H. The Linear and Volume Compresssion of β-HMX and RDX to 9 GPa (90 kilobar). In Actes du Symposium International sur le Comportement des Milieux Denses sous Hautes Pressions Dynamiques; Commmissariat à l’Energie Atomique: Paris, 1978; pp 3−8. (11) Davidson, A. J.; Oswald, I. D. H.; Francis, D. J.; Lennie, A. R.; Marshall, W. G.; Millar, D. I. A.; Pulham, C. R.; Warren, J. E.; Cumming, A. S. Explosives under PressureThe Crystal Structure of γ-RDX as Determined by High-Pressure X-ray and Neutron Diffraction. CrystEngComm 2008, 10, 162−165. (12) Ciezak, J. A.; Jenkins, T. A. The Low-Temperature HighPressure Phase Diagram of Energetic Materials: I. Hexahydro-1,3,5Trinitro-S-Triazine. Propellants, Explos., Pyrotech. 2008, 33, 390−395. (13) Millar, D. I. A.; Oswald, I. D. H.; Barry, C.; Francis, D. J.; Marshall, W. G.; Pulham, C. R.; Cumming, A. S. Pressure-Cooking of ExplosivesThe Crystal Structure of ε-RDX as Determined by X-ray and Neutron Diffraction. Chem. Commun. 2010, 46, 5662−5664. (14) Dzyabchenko, A. V.; Pivina, T. S.; Arnautova, E. A. Prediction of Structure and Density for Organic Nitramines. J. Mol. Struct.: THEOCHEM 1996, 378, 67−82. (15) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. Intermolecular Potential for the Hexahydro-1,3,5-Trinitro-1,3,5-S-Triazine Crystal (RDX): A Crystal Packing, Monte Carlo, and Molecular Dynamics Study. J. Phys. Chem. B 1997, 101, 798−808. (16) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. A Transferable Intermolecular Potential for Nitramine Crystals. J. Phys. Chem. A 1998, 102, 8386−8392. (17) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. IsothermalIsobaric Molecular Dynamics Simulations of 1,3,5,7-Tetranitro-1,3,5,7Tetraazacyclooctane (HMX) Crystals. J. Phys. Chem. B 1998, 102, 6692−6695. (18) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. Molecular Packing and NPT Molecular Dynamics Investigation of the Transferability of the RDX Intermolecular Potential to 2,3,6,8,10,12-Hexanitrohexaazaisowurtzitane. J. Phys. Chem. B 1998, 102, 948−952. (19) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. Theoretical Studies of the Hydrostatic Compression of RDX, HMX, HNIW, and PETN Crystals. J. Phys. Chem. B 1999, 103, 6783−6790. (20) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. Molecular Packing and Molecular Dynamics Study of the Transferability of a Generalized J

DOI: 10.1021/acs.jpcc.6b06415 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Nitrogen Energetic Salts. J. Chem. Theory Comput. 2014, 10, 4982− 4994. (39) von Lilienfeld, O. A.; Tavernelli, I.; Rothlisberger, U.; Sebastiani, D. Optimization of Effective Atom Centered Potentials for London Dispersion Forces in Density Functional Theory. Phys. Rev. Lett. 2004, 93, 153004. (40) Klimeš, J.; Bowler, D. R.; Michaelides, A. Van der Waals Density Functionals Applied to Solids. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 195131. (41) Vydrov, O. A.; Van Voorhis, T. Nonlocal van Der Waals Density Functional: The Simpler the Better. J. Chem. Phys. 2010, 133, 244103. (42) Balu, R.; Byrd, E. F. C.; Rice, B. M. Assessment of Dispersion Corrected Atom Centered Pseudopotentials: Application to Energetic Molecular Crystals. J. Phys. Chem. B 2011, 115, 803−810. (43) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. Comput. Phys. Commun. 2005, 167, 103−128. (44) Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. CP2K: Atomistic Simulations of Condensed Matter Systems. WIREs Comput. Mol. Sci. 2014, 4, 15−25. (45) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (46) Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703−1710. (47) Krack, M. Pseudopotentials for H to Kr Optimized for GradientCorrected Exchange-Correlation Functionals. Theor. Chem. Acc. 2005, 114, 145−152. (48) Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic Separable Dual-Space Gaussian Pseudopotentials from H to Rn. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 3641−3662. (49) Erpenbeck, J. J. Molecular-Dynamics of Detonation. I. Equation of State and Hugoniot Curve for a Simple Reactive Fluid. Phys. Rev. A: At., Mol., Opt. Phys. 1992, 46, 6406−6416. (50) Brennan, J. K.; Rice, B. M. Efficient Determination of Hugoniot States Using Classical Molecular Simulation Techniques. Mol. Phys. 2003, 101, 3309−3322. (51) Ilyukhin, V. S.; Pokhil, P. F.; Rozanov, O. K.; Shvedova, N. S. Measurement of the Shock Adiabatics of Cast Trinitrotoluene, Crystalline Hexogene, and Nitromethane. Doklady Akad. Nauk SSSR 1960, 131, 793−796.

K

DOI: 10.1021/acs.jpcc.6b06415 J. Phys. Chem. C XXXX, XXX, XXX−XXX