Equations of state for the α and γ polymorphs of the energetic molecular crystal cyclotrimethylene trinitramine (RDX) have been developed from their Helmholtz free energies. The ion motion contribution to the Helmholtz free energy is represented by Debye models with density-dependent Debye temperatures that are parameterized to vibrational densities of states computed from dispersion-corrected density functional theory. By separating the vibrational density of states into low frequency modes of mainly lattice phonon character and high frequency modes of intramolecular character we were able to significantly improve the description of the heat capacity at low temperatures and the thermal contribution to the pressure. The ion motion contribution to the Helmholtz free energy of the high pressure γ polymorph was constructed from that of the α polymorph to reproduce the temperature-independent transformation pressure seen experimentally. The static lattice energies for both polymorphs were constructed to reproduce published isothermal compression data. The equations of state have been applied to the prediction of the path of the principal Hugoniot in the equilibrium phase diagram.

Cyclotrimethylene trinitramine (RDX), with chemical formula C3H6N6O6, is a heavily used and extensively studied energetic molecular crystal. The phase diagram of RDX in temperature and pressure space contains at least four polymorphs.1–5 The α polymorph is the most stable polymorph under ambient conditions and comprises an orthorhombic unit cell in space group Pbca containing eight symmetry-related molecules.6,7 The γ polymorph was discovered in the late 1970s in room temperature Bridgman anvil experiments at hydrostatic pressures in excess of about 4 GPa.8 While the γ polymorph was identified as having orthorhombic symmetry in the original work, its Pca21 space group was not solved until 2008.9 The α and γ polymorphs of RDX are structurally closely related. Both polymorphs contain eight molecules per unit cell, but whereas all molecules in the α polymorph are symmetry related, the γ polymorph contains two independent molecules.

The ϵ polymorph is found at temperatures in excess of about 450 K and pressures in the interval 2 < P < 6 GPa.10 Like the γ polymorph it belongs to space group Pca21 but contains four molecules per unit cell rather than eight. The δ polymorph has been inferred from spectroscopic measurements above about 18 GPa.1 However, structural studies on this polymorph have not been possible since it requires a pressure transmitting medium that remains hydrostatic at those high pressures. The β polymorph, which is highly metastable under ambient conditions, is obtained from crystallization from boiling volatile solvents.6,11 It too adopts space group Pca21 with eight molecules per unit cell. Changes in the conformation of RDX molecules between the polymorphs in addition to changes in space group and molecular packing have been measured. The conformation of the three nitro groups with respect to the triazine ring for molecules in the α polymorph is denoted “AAE” since two groups are oriented axially (A) and one equatorially (E). The conformation of molecules in the β and ϵ polymorphs are denoted AAA since all three N-N bonds are close to perpendicular to the plane of the ring. The conformation of molecules in the γ polymorph is denoted AAI since the orientation of one nitro group is intermediate (I), between the axial and equatorial positions.5 

The dynamic compression of RDX generates thermodynamic states of elevated temperature and pressure. The path taken by the principal (single shock) Hugoniot through the equilibrium pressure/temperature phase diagram is of particular interest for the interpretation of plate impact experiments and for predicting from which crystallographic phase detonation initiates.12 The capability to compute temperature, pressure, and density under mechanical loading via an equation of state is important in the development of thermomechanical models, especially if they are to be applied in the study of high rate phenomena.13–15 For instance, in Ref. 14, an equation of state provides the non-linear elastic part of the stress tensor.

The heat capacity of organic materials such as RDX at room temperature is small compared with the classical limit, and strongly dependent on temperature. Models using a constant heat capacity based on the measured value at room temperature or classical molecular dynamics will overestimate or strongly underestimate shock heating, respectively. Hence, it is desirable that the free energies and equations of state used for energetic materials capture accurately the temperature dependence of the heat capacity.

We have generated equations of state for α and γ RDX based on Debye theory within the quasi-harmonic approximation16,17 using isothermal compression data for both phases, dispersion-corrected density functional theory calculations of the phonon densities of states as a function of specific volume, and the temperature dependence of the α–γ phase transformation pressure. We used two density-dependent Debye temperatures in our models to represent the low frequency lattice phonons and high frequency intramolecular vibrational modes. The resulting free energies give a good representation of the heat capacity and Grüneisen parameter under ambient conditions. The equations of state are computationally simple and straightforward to implement in finite element simulation codes.

The Helmholtz free energy, F(V, T), of a crystal containing N atoms with 3(N − 1) vibrational degrees of freedom is

F(V,T)=Φ0(V)+Fion(V,T),
(1)

where Φ0(V) is the temperature-independent static lattice energy, Fion(V,T) the ion motion free energy, V the specific volume, and T the temperature. The ion motion free energy depends on the vibrational density of states

Fion=12i=13(N1)ωi+kBTi=13(N1)ln(1exp(ωi/kBT)),
(2)

where is the reduced Planck constant, kB Boltzmann's constant, and ωi the angular frequency of the i-th normal mode of the crystal. The first term on the right is the zero point energy of the set of harmonic oscillators. We approximate the ion motion free energy with a Debye model, that is, Fion ≈ FD, where

FD=9N[18kBθD+kBTY(θDT)],
(3)

θD is the Debye temperature that characterizes the set of vibrational modes, and

Y(x)=(1x)30xz2ln(1ez)dz.
(4)

Within the quasi-harmonic approximation the Debye temperature is a function of the specific volume, that is, θD=θD(V).16 We compute the Debye temperature and its dependence on density by numerically fitting the heat capacity derived from the Debye model

CVD(V)=9NkB(TθD(V))30θD(V)/Tz4ez(ez1)2dz,
(5)

to heat capacities computed directly from sets of normal mode frequencies18 that are extracted from density functional theory calculations on crystals at various densities, {ω(V)}, where

CV(V)=kBi=13(N1)(ωi(V)/kBT)exp(ωi(V)/kBT)(1exp(ωi(V)/kBT))2.
(6)

Hence, the Debye temperature at a given specific volume, θD(V), is that which minimizes the objective function

χ2=0(CVD(θD(V))CV(V))2dT.
(7)

The vibrational density of states of organic molecular crystals typically comprise a continuum from zero to several hundred cm−1 that are mainly of lattice phonon character and a set of discrete intramolecular vibrational modes up to about 3200 cm−1. The thermal pressure is dominated by the low frequency lattice phonons whereas the large classical limit and strong temperature dependence of the heat capacity arises from the intramolecular modes. Owing to the different character and dependence on the specific volume of the lattice phonons and intramolecular vibrational modes, we represented the ion motion free energy by two separate Debye models, Fion(V)=FD,1(V)+FD,2(V). We assign the Q normal modes of lowest frequency of mostly lattice phonon character to FD,1 and the remaining modes to FD,2, that is,

FD,1(V)=3Q[18kBθD,1(V)+kBTY(θD,1(V)T)],
(8)
FD,2(V)=3(3NQ)[18kBθD,2(V)+kBTY(θD,2(V)T)],
(9)

where θD,1(V) and θD,2(V) are the density-dependent Debye temperatures for the two populations of normal modes. We evaluate both Debye temperatures and their dependence on density by minimizing the objective function, Eq. (7), for each set of normal modes, 1Q and Q+13(N1), obtained from first principles calculations.

We represent the dependence of the Debye temperatures on specific volume using Greeff's approach17,19,20 where the volume dependence of the Grüneisen parameter, γ=dlnθD/dlnV, is

γ(V)=γ+A(VV0)+B(VV0)2,
(10)

where γ=2/3 is the Grüneisen parameter at infinite compression, V0 is a reference volume, and A and B are adjustable parameters. With this representation of the Grüneisen parameter, the volume dependence of the Debye temperature is

θD(V)=θD(V0)(VV0)γexp[A((VV0)1)B2((VV0)21)].
(11)

Equation (11) is parameterized for both populations of normal modes by a least squares fit to Debye temperatures calculated from the vibrational density of states.

The static lattice energy, Φ0(V), or cold curve, is represented by the Vinet form17,21

Φ0(V)=4V*B*(B1*1)2[1(1+X)exp(X)]+ΔΦ,
(12)

where

X=32(B1*1)[(VV*)1/31],
(13)

and V* is the specific volume where Φ0(V)/V=0,B* the bulk modulus, B1* its pressure derivative, and ΔΦ a constant energy shift. The cold curves for both phases are parameterized to isothermal compression data, that is, density as a function of hydrostatic pressure. For the equations of state to reproduce isothermal compression data from experiment we require

Pref(V,T)=P0(V)+PD(V,T),
(14)

where Pref(V,T) are the measured hydrostatic pressure and specific volume at temperature T, PD(V,T) is the pressure derived from the ion motion model, and P0(V) the contribution to the hydrostatic pressure from the cold curve. Hence, by removing the contribution to the pressure along the measured isothermal compression curve from the ion motion free energy, i.e., P0=PrefPD, the cold curve, Eq. (12) is parameterized by a least squares fit of the corresponding pressure as a function of specific volume

P0(V)=4V*B*(B1*1)2XXVexp(X),
(15)

where

XV=B1*12(V2V*)1/3.
(16)

Parameterizing the equations of state to experimental isothermal compression data ensures that the models correctly reproduce the mass density under ambient conditions.

All structural optimizations and calculations of the phonon densities of states were computed using dispersion-corrected density functional theory. Numerous groups have demonstrated the usefulness of empirical dispersion corrections for obtaining accurate specific volumes for systems where weak van der Waals interactions contribute significantly to cohesion.22,23

The normal mode frequencies of α RDX were computed using the D3(BJ) dispersion correction of Grimme et al.24 with the DZVP basis set, Goedecker–Teter–Hutter pseudopotentials,25 and the Perdew–Burke–Ernzerhof exchange-correlational functional26 within the cp2k package of codes.27 A plane wave cut-off of 1200 Ry with a relative cut-off (REL_CUTOFF) of 60 Ry gave both a well-converged total energy and a good distribution of the Gaussian basis functions between the four integration grids. The use of larger basis sets for the evaluation of the normal modes frequencies led to prohibitively slow calculations.

We first optimized the lattice parameters and atomic coordinates of α RDX at zero external pressure while maintaining orthorhombic symmetry using the structure published by Choi and Prince as the starting configuration.7 The relaxation was stopped when the maximum force acting on any atom fell below 5.14 × 10−7 eV/Å. The resulting lattice parameters and specific volume were a = 13.319, b = 11.347, c = 10.703 Å, and V0 = 0.5481834 cm3/g, which are in excellent agreement with the room temperature experimental values of a = 13.182, b = 11.574, and c = 10.709 Å.

The normal mode frequencies were computed starting from the optimized structure by numerically calculating all force constants and solving the dynamical matrix using the native functionality of the cp2k code. The volume dependence of the normal mode frequencies was determined by homogeneously compressing the α RDX lattice in 21 increments between zero and 40% compression, i.e., V/V0=0.6 and re-computing the force constants after the internal coordinates were relaxed until the maximum force acting on any atom was less than 5.14 × 10−7 eV/Å. The homogeneous compression of the lattice leads to non-hydrostatic stress states owing to the anisotropy of the lattice.

In Fig. 1, we plot the smeared, normalized phonon density of states

n(f)=1πj=13(N1)Im(1ffj+iη),
(17)

where f is the frequency, fj the frequency of the j-th normal mode, i=1, and η a positive real number that controls the smearing of the density of states for plotting purposes. The inset of Fig. 1 shows the vibrational density of states in the interval between zero frequency and 400 cm−1. The continuous band of lattice phonons is clearly visible from about 20 to 200 cm−1. At zero applied pressure, the lattice phonon band, f ≤ 200 cm−1, comprises the 93 normal modes of lowest frequency per unit cell. During the subsequent analysis of the dependence of the normal mode frequencies on specific volume, we always assign the first 93 modes to the lattice phonon population and the remaining 408 normal modes to the set of intramolecular vibrations.

FIG. 1.

Vibrational density of states of α RDX at zero external pressure computed from dispersion-corrected density functional theory. The inset highlights the frequency interval from 0 to 400 cm−1. The 200 cm−1 cut-off used to separate the lattice phonons and intramolecular vibrations is highlighted by the broken vertical line.

FIG. 1.

Vibrational density of states of α RDX at zero external pressure computed from dispersion-corrected density functional theory. The inset highlights the frequency interval from 0 to 400 cm−1. The 200 cm−1 cut-off used to separate the lattice phonons and intramolecular vibrations is highlighted by the broken vertical line.

Close modal

The heat capacity, CV, of α RDX at 300 K and zero pressure determined directly from the normal mode frequencies via Eq. (6) is 217.1 J/mol/K. This value compares very well with the calculated value in Ref. 28 of 216.2 J/mol/K. Values for CV at room temperature derived from experimental measurements of the heat capacity at constant pressure, CP, volumetric thermal expansion coefficient, isothermal bulk modulus, and density span 200.1–205.4 J/mol/K.29,30

The heat capacity derived from a Debye model with a single Debye temperature, θD = 1657 K, fitted to the full set of normal mode frequencies computed at zero pressure is presented in Fig. 2. The heat capacity at 300 K computed from the simple, single θD Debye model is only CVD=161.9J/mol/K. We also present in Fig. 2 the heat capacity obtained from the extended model with separate Debye models for the phonons and intramolecular vibrational modes. It is clear that the extended model gives a better representation of the temperature dependence of the heat capacity, particularly at low temperatures. The description of the heat capacity under ambient conditions is also in better agreement with experiment than that obtained from the simple model with CVD=178.4J/mol/K.

FIG. 2.

Temperature dependence of the heat capacity of RDX derived from the normal mode frequencies computed at zero pressure. The exact result refers to the direct calculation of the heat capacity from the normal mode frequencies, and the simple and extended models refer to the Debye models with one and two Debye temperatures, respectively.

FIG. 2.

Temperature dependence of the heat capacity of RDX derived from the normal mode frequencies computed at zero pressure. The exact result refers to the direct calculation of the heat capacity from the normal mode frequencies, and the simple and extended models refer to the Debye models with one and two Debye temperatures, respectively.

Close modal

The dependences of the Debye temperatures on specific volume for the lattice phonons and intramolecular vibrational modes are presented in Fig. 3. As expected, the Debye temperature derived from the vibrational modes of mainly lattice phonon character depends far more sensitively on the specific volume than the intramolecular modes. The parameters in Eq. (11) as determined by a least squares fit of these data are presented in Table I.

FIG. 3.

Dependence of the Debye temperatures computed from the lattice phonons and intramolecular vibrational modes on specific volume.

FIG. 3.

Dependence of the Debye temperatures computed from the lattice phonons and intramolecular vibrational modes on specific volume.

Close modal
TABLE I.

Parameterization of the dependence of the two Debye temperatures, θD,1 and θD,2, on specific volume via Eq. (11). The reference specific volume, V0, for the α and γ polymorphs are 0.5481834 and 0.5366715 cm3/g, respectively.

θD(V0) (K)AB
θD,1 188.783 2.00878 −0.805669 
θD,2 2073.6 −0.545878 −0.0991197 
θD(V0) (K)AB
θD,1 188.783 2.00878 −0.805669 
θD,2 2073.6 −0.545878 −0.0991197 

The static lattice energy, Φ0(V), for α RDX was fitted to the room temperature isothermal compression data provided by Olinger et al.8 and Davidson et al.9 by subtracting from the measured pressure at each specific volume the contribution to the pressure from the ion motion term, followed by a least squares fit of Eq. (15). The parameterization of the Vinet cold curve for the α polymorph is given in Table II and the 300 K isothermal compression curves from experiment and the equation of state are presented in Fig. 4.

TABLE II.

Parameterization the Vinet cold curves, Φ0(V), for the α and γ polymorphs.

αγ
V* (cm3/g) 0.532138 0.539957 
B1 (GPa) 13.6699 7.73837 
B1* 9.25677 11.3446 
ΔΦ (J/g) 0.0 25.520 
αγ
V* (cm3/g) 0.532138 0.539957 
B1 (GPa) 13.6699 7.73837 
B1* 9.25677 11.3446 
ΔΦ (J/g) 0.0 25.520 
FIG. 4.

Room temperature isothermal compression data for α and γ RDX. The symbols represent experimental data from Refs. 8 and 9, and the solid lines are derived from the equations of state.

FIG. 4.

Room temperature isothermal compression data for α and γ RDX. The symbols represent experimental data from Refs. 8 and 9, and the solid lines are derived from the equations of state.

Close modal

The ion motion contribution to the Helmholtz free energy of the γ polymorph is constrained heavily by the experimental observation3,4 of the independence of the transformation pressure, Pt, between the α and γ polymorphs on temperature, i.e.,

PtT0.
(18)

We found that applying the procedure described earlier to generate the ion motion free energy for the γ polymorph gives a transformation pressure that depends sensitively on temperature. Hence, to ensure that the free energies satisfy Eq. (18), we derived the ionic motion free energy for the γ polymorph from that of the α polymorph.

The phase transformation pressure, Pt, is the pressure at which the Gibbs free energies of the two phases are equal, i.e.

Fα(Vα,t,T)+PtVα,t=Fγ(Vγ,t,T)+PtVγ,t,
(19)

or

Pt=Fγ(Vγ,t,T)Fα(Vα,tT)ΔVt,
(20)

where Fα and Fγ are the Helmholtz free energies of the α and γ phases, respectively. Vα,t and Vγ,t are the specific volumes of the α and γ phases at the transition, respectively, and ΔVt=Vα,tVγ,t. Since we require the transformation pressure to satisfy Eq. (18), then

1ΔVt(Fγ(Vγ,t,T)TFα(Vα,t,T)T)+(Fγ(Vγ,t,T)Fα(Vα,t,T))TΔVt1=0.
(21)

Assuming that the change in specific volume during the transformation does not depend on the temperature, i.e.,

TΔVt1=0,
(22)

we require

Fγ(Vγ,t,T)T=Fα(Vα,t,T)T.
(23)

Since the cold curves, Φ0(V), do not depend on the temperature, we require that the temperature derivatives of the ion motion terms for the two polymorphs are equal at Pt.

We satisfied Eq. (23) to achieve a temperature-independent transformation pressure by adapting the Debye models developed for α RDX to the γ polymorph. Inspecting the room temperature isothermal compression data of Olinger et al.8 and Davidson et al.,9 we find that at the transformation pressure, Pt = 3.8 GPa,

Vγ,tVα,t0.979.
(24)

Hence, the two-temperature Debye model [Eqs. (8), (9), and (11)] for the γ polymorph is identical to that of the α polymorph (same θD, A, B, and Q) except for a 2.1% reduction of the value of V0 in Eq. (11). Therefore, the parameters given in Table I for α RDX, and the number of normal modes assigned to lattice phonons and intramolecular vibrations, Q = 93, are unchanged for γ RDX. The Debye model for the γ phase instead uses V0 = 0.5366715 cm3/g. Since the α and γ polymorphs are structurally and chemically closely related, we anticipate that the thermal model for the γ phase provides a good description of the temperature dependence of its physical properties, despite not being explicitly parameterized for that phase.

There is some ambiguity in the literature regarding the value of the transformation pressure. For instance, Dreger and Gupta give Pt = 3.7 GPa while Davidson et al. give 3.9 GPa.4,9 A value of 3.8 GPa is consistent with Ref. 4 since a Raman spectrum that unambiguously corresponds to the γ phase is obtained at 3.8 GPa. The neutron powder diffraction patterns presented by Davidson et al. show α RDX at 3.62 GPa and γ RDX at 3.9 GPa. Hence, we have opted to construct the equation of state of the γ polymorph based on the transformation pressure Pt = 3.8 GPa.

The static lattice energy, Φ0(V), for the γ polymorph was parameterized to the isothermal compression data of Olinger et al.8 and Davidson et al.9 following the removal of the contribution from the ion motion free energy. Owing to the absence of experimental data on the γ polymorph for P < Pt, we used dispersion-corrected density functional theory calculations to constrain the parameterization of the Vinet form. Zero temperature hydrostatic compression curves were computed using the same dispersion correction and exchange-correlation functional used in our calculations of normal mode frequencies. However, since static compression curves are much less computationally expensive than normal mode calculations, we instead used the larger TZV2P basis set with a converged plane wave cut-off of 1500 Ry. The structure determined by Davidson et al. at 5.2 GPa was used as the starting point for the structural optimization of the γ polymorph.9 The zero temperature compression curves presented in Fig. 5 were determined by optimizing the unit cell parameters and internal coordinates under external hydrostatic pressures in the interval 0 ≤ P ≤ 10 GPa. The potential energies, E, presented in Fig. 5 are relative to that of the α polymorph at zero pressure, E0. The crystal structure of the γ polymorph was found to be at least metastable until a hydrostatic pressure of 2 GPa below which it transformed to the α phase during the relaxation. The geometry optimizations were terminated when the maximum force acting on any atom fell below 5.14 × 10−7 eV/Å. Least squares fits of the Vinet equation of state, Eq. (12), to the computed zero temperature compression curves infer that the equilibrium volume of the γ polymorph is about 1.47% larger than that of the α polymorph. Hence, we constrained V* = 0.539957 cm3/g in the static lattice energy term, Φ0(V), for the γ polymorph when determining by a least squares fit of the B* and B1* parameters. The energy shift, ΔΦ, in Eq. (12) was determined by requiring Gα=Gγ at Pt = 3.8 GPa. The resulting parameterization is presented in Table II and the computed 300 K isotherm in Fig. 4.

FIG. 5.

Energy as a function of specific volume at zero temperature for α and γ RDX computed using dispersion-corrected density functional theory with TZV2P basis sets. The compression curves are used to estimate the specific volume of the γ polymorph at zero pressure.

FIG. 5.

Energy as a function of specific volume at zero temperature for α and γ RDX computed using dispersion-corrected density functional theory with TZV2P basis sets. The compression curves are used to estimate the specific volume of the γ polymorph at zero pressure.

Close modal

The dependence of the transformation pressure, Pt, on temperature is presented in Fig. 6 where we show the difference in Gibbs free energy between the two polymorphs, ΔG=GαGγ, as a function of temperature and pressure. Here, the difference in the Gibbs free energy is zero at P = 3.8 GPa for all curves in the interval 50 ≤ T ≤ 450 K, indicating that Eq. (18) is fully satisfied in this temperature interval.

FIG. 6.

Difference between the Gibbs free energies, ΔG, of the α and γ polymorphs as a function of pressure for temperatures in the interval 50 ≤ T ≤ 450 K.

FIG. 6.

Difference between the Gibbs free energies, ΔG, of the α and γ polymorphs as a function of pressure for temperatures in the interval 50 ≤ T ≤ 450 K.

Close modal

The Grüneisen parameter of α RDX under ambient conditions derived from experimental data, γ=αVBV/CV, where αV is the volumetric thermal expansion coefficient, and B the isothermal bulk modulus, is γ = 1.34.31–33 The two Debye models describing the ion motion free energy give very different contributions to the total Grüneisen parameter. Using Eq. (10) and Table I, the contributions to the Grüneisen parameter, γ=V(P/F)V, from the lattice phonons and intramolecular modes at 300 K and zero pressure are γ = 1.875 and 0.01122, respectively. Under these conditions, the total free energy (Eqs. (8) and (9)) gives a total Grüneisen parameter of 1.0, i.e., a value roughly midway between the contributions from FD,1 and FD,2. The dependence of the Grüneisen parameter on temperature and specific volume is presented in Fig. 7. We have plotted the Grüneisen parameter from the RDX polymorph with the lowest Gibbs free energy so that the values for specific volumes less than about 0.465 cm3/g correspond to γ RDX.

FIG. 7.

Dependence of the Grüneisen parameter on specific volume and temperature for the RDX polymorph with the lowest Gibbs free energy under those conditions.

FIG. 7.

Dependence of the Grüneisen parameter on specific volume and temperature for the RDX polymorph with the lowest Gibbs free energy under those conditions.

Close modal

The origin of the low value of the Grüneisen parameter obtained from the free energy is its underestimation of the volumetric thermal expansion coefficient, αV. The experimental value for αV under ambient conditions is 1.935 × 10−4 K−1 whereas the free energy gives 1.271 × 10−4 K−1.33 Nevertheless, the prediction from the free energy is in good agreement with the coefficient of thermal expansion derived from Monte Carlo simulations in the isothermal-isobaric ensemble by Sewell and Bennett, αV = 1.25 × 10−4 K−1.34 Figure 8 depicts the dependence of the coefficient of volumetric thermal expansion on specific volume and temperature derived from the free energy model.

FIG. 8.

Dependence of the volumetric thermal expansion coefficient, αV, of α and γ RDX on specific volume and temperature.

FIG. 8.

Dependence of the volumetric thermal expansion coefficient, αV, of α and γ RDX on specific volume and temperature.

Close modal

The free energy provides accurate descriptions of the bulk modulus, B=V(2F/V2)T, and the pressure dependence of the bulk modulus at room temperature since it was parameterized to high quality isothermal compression data. The temperature dependence of the bulk modulus, depicted in Fig. 9, at zero pressure from the free energy is ∂B/∂T = −0.0106 GPa/K. Experimental data provided by Haussühl enables the thermal softening predicted by our free energy model to be compared with measurements.31 The Voigt average bulk modulus of an orthorhombic crystal is

BV=19(C11+C22+C33)+29(C12+C13+C23),
(25)

where Cij are single crystal elastic constants in Voigt notation. The temperature derivative of BV following Haussühl's notation where Tij=lnCij/T is

BVT=19(C11T11+C22T22+C33T33)+29(C12T12+C13T13+C23T23).
(26)

Using the values tabulated by Haussühl for Cij and Tij, we find BV = 11.36 GPa and ∂BV/∂T = −0.0148 GPa/K.31 Hence, our model for the Helmholtz free energy provides a reasonable description of the temperature dependence of the bulk modulus.

FIG. 9.

The bulk modulus of α RDX as a function of temperature at zero pressure computed from the free energy.

FIG. 9.

The bulk modulus of α RDX as a function of temperature at zero pressure computed from the free energy.

Close modal

The principal Hugoniot calculated from the equations of state and the Rankine–Hugoniot jump conditions is overlaid on phase diagrams measured by Dreger and Gupta4 and Baer et al.3 in Fig. 10. Unfortunately, the close to 100 K discrepancy between the triple point among the α, γ, and ϵ phases in the published phase diagrams prevents us from discerning, based on volumetric strain, which of the γ or ϵ phases play important roles in initiation in RDX. Recent work by Dreger and Gupta12 clearly identified the decomposition of the γ polymorph under static high temperature, high pressure conditions and also the role of kinetics and heating rate on the decomposition onset temperature. Since Baer et al. do not report the heating rate used in their experiments, we cannot be sure whether or not the discrepancies between the two phase diagrams arise from kinetic effects. Plastic work, i.e., shear strain at constant volume, also contributes to heating during dynamic compression. Given the potential location of the triple point between the α, γ, and ϵ phases in P-T space, especially in the phase diagram of Baer et al., the contribution of plasticity to the temperature during shock experiments should also be assessed in detail. The relatively small contribution to the heating arising from collapsing porosity in real materials could also significantly affect the phase response. Hence, in situ X-ray diffraction experiments during shock compression are likely to be invaluable for identifying transient high pressure phases in complex materials.35,36

FIG. 10.

Pressure–temperature phase diagram of RDX from Dreger and Gupta4 and Baer et al.3 Note that we have reassigned the β phase in the phase diagram from Baer et al. to the ϵ phase. The path followed by the principle Hugoniot determined from the equation of state is overlaid.

FIG. 10.

Pressure–temperature phase diagram of RDX from Dreger and Gupta4 and Baer et al.3 Note that we have reassigned the β phase in the phase diagram from Baer et al. to the ϵ phase. The path followed by the principle Hugoniot determined from the equation of state is overlaid.

Close modal

We have developed computationally simple expressions for the equations of state for the α and γ polymorphs of the energetic molecular crystal RDX based on parameterized Helmholtz free energies. The free energies are parameterized to normal mode frequencies calculated from dispersion-corrected density functional theory, experimental isotherms, and the RDX phase diagram. They provide good descriptions of the thermal properties of RDX, in particular, the temperature dependence of the heat capacity. The equations of state were applied in predicting the path of the principal Hugoniot through the equilibrium phase diagram. Unfortunately, the large discrepancies between the published phase diagrams prevent us from assessing, based solely on volumetric strain, to what extent the γ and/or ϵ phases influence initiation under impact.

This work was supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory. M.J.C. thanks Carl Greeff for helpful discussions.

1.
J. A.
Ciezak
,
T. A.
Jenkins
,
Z.
Liu
, and
R. J.
Hemley
,
J. Phys. Chem. A
111
,
59
(
2007
).
2.
J. A.
Ciezak
and
T. A.
Jenkins
,
Propellants, Explos., Pyrotech.
33
,
390
(
2008
).
3.
B. J.
Baer
,
J.
Oxley
, and
M.
Nicol
,
High Pressure Res.
2
,
99
(
1990
).
4.
Z. A.
Dreger
and
Y. M.
Gupta
,
J. Phys. Chem. A
114
,
8099
(
2010
).
5.
D. I. A.
Millar
, “
Energetic materials at extreme conditions
,” Ph.D. thesis (
University of Edinburgh
,
2010
).
6.
W. C.
McCrone
,
Anal. Chem.
22
,
954
(
1950
).
7.
C. S.
Choi
and
E.
Prince
,
Acta Cryst. B
28
,
2857
(
1972
).
8.
B.
Olinger
,
B.
Roof
, and
H.
Cady
,
Behavior of Dense Media under High Dynamic Pressures
(
Commissariat á l'Energie Atomique
,
1978
) pp.
3
8
.
9.
A. J.
Davidson
,
I. D. H.
Oswald
,
D. J.
Francis
,
A. R.
Lennie
,
W. G.
Marshall
,
D. I. A.
Millar
,
C. R.
Pulham
,
J. E.
Warren
, and
A. S.
Cumming
,
CrystEngComm
10
,
162
(
2008
).
10.
D. I. A.
Millar
,
I. D. H.
Oswald
,
C.
Barry
,
D. J.
Francis
,
W. G.
Marshall
,
C. R.
Pulham
, and
A. S.
Cumming
,
Chem. Commun.
46
,
5662
(
2010
).
11.
D. I. A.
Millar
,
I. D. H.
Oswald
,
D. J.
Francis
,
W. G.
Marshall
,
C. R.
Pulham
, and
A. S.
Cumming
,
Chem. Commun.
562
(
2009
).
12.
Z. A.
Dreger
and
Y. M.
Gupta
,
J. Phys. Chem. A
116
,
8713
(
2012
).
13.
J. M.
Winey
and
Y. M.
Gupta
,
J. Appl. Phys.
107
,
103505
(
2010
).
14.
S.
De
,
A. R.
Zamiri
, and
Rahul
,
J. Mech. Phys. Solids
64
,
287
(
2014
).
15.
J. D.
Clayton
and
R.
Becker
,
J. Appl. Phys.
111
,
063512
(
2012
).
16.
D. C.
Wallace
,
Thermodynamics of Crystals
(
Dover
,
1998
).
17.
C. W.
Greeff
and
M. J.
Graf
,
Phys. Rev. B: Condens. Matter Mater. Phys.
69
,
054107
(
2004
).
18.
D. A.
McQuarrie
,
Statistical Mechanics
(
University Science Books
,
2000
).
19.
C. W.
Greeff
,
Modell. Simul. Mater. Sci. Eng.
13
,
1015
(
2005
).
20.
C. W.
Greeff
,
J. C.
Boettger
,
M. J.
Graf
, and
J. D.
Johnson
,
J. Phys. Chem. Solids
67
,
2033
(
2006
).
21.
P.
Vinet
,
J.
Ferrante
,
J. R.
Smith
, and
J. H.
Rose
,
J. Phys. C
19
,
L467
(
1986
).
22.
S.
Grimme
,
J.
Antony
,
T.
Schwabe
, and
C.
Muck-Lichtenfeld
,
Org. Biomol. Chem.
5
,
741
(
2007
).
23.
D. C.
Sorescu
and
B. M.
Rice
,
J. Phys. Chem. C
114
,
6734
(
2010
).
24.
S.
Grimme
,
S.
Ehrlich
, and
L.
Goerigk
,
J. Comput. Chem.
32
,
1456
(
2011
).
25.
S.
Goedecker
,
M.
Teter
, and
J.
Hutter
,
Phys. Rev. B
54
,
1703
(
1996
).
26.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
27.
J.
VandeVondele
,
M.
Krack
,
F.
Mohamed
,
M.
Parrinello
,
T.
Chassaing
, and
J.
Hutter
,
Comput. Phys. Commun.
167
,
103
(
2005
).
28.
S.
Hunter
,
T.
Sutinen
,
S. F.
Parker
,
C. A.
Morrison
,
D. M.
Williamson
,
S.
Thompson
,
P. J.
Gould
, and
C. R.
Pulham
,
J. Phys. Chem. C
117
,
8062
(
2013
).
29.
LASL Explosives Property Data
, edited by
T. R.
Gibbs
and
A.
Popolato
(
University of California Press
,
1980
).
30.
M. S.
Miller
,
J. Thermophys. Heat Transfer
8
,
803
(
1994
).
31.
S.
Haussühl
,
Z. Kristallogr.
216
,
339
(
2001
).
32.
C. A.
Bolme
and
K. J.
Ramos
,
J. Appl. Phys.
116
,
183503
(
2014
).
33.
H. H.
Cady
,
J. Chem. Eng. Data
17
,
369
(
1972
).
34.
T. D.
Sewell
and
C. M.
Bennett
,
J. Appl. Phys.
88
,
88
(
2000
).
35.
K. J.
Ramos
,
B. J.
Jensen
,
A. J.
Iverson
,
J. D.
Yeager
,
C. A.
Carlson
,
D. S.
Montgomery
,
D. G.
Thompson
,
K.
Fezzaa
, and
D. E.
Hooks
,
J. Phys.: Conf. Ser.
500
,
142028
(
2014
).
36.
A. E.
Gleason
,
C. A.
Bolme
,
H. J.
Lee
,
B.
Nagler
,
E.
Galtier
,
D.
Milathianaki
,
J.
Hawreliak
,
R. G.
Kraus
,
J. H.
Eggert
,
D. E.
Fratanduono
,
G. W.
Collins
,
R.
Sandberg
,
W.
Yang
, and
W. L.
Mao
,
Nat. Commun.
6
,
8191
(
2015
).