# Pressure, temperature, and orientation dependent thermal conductivity
of $\alpha$-1,3,5-trinitro-1,3,5-triazinane ($\alpha$-RDX)

Romain Perriot<EMAIL_ADDRESS>Theoretical Division, Los Alamos National
Laboratory, Los Alamos, NM 87545, USA  Michael S. Powell  High Explosives
Science and Technology, Los Alamos National Laboratory, Los Alamos, NM, 87545,
USA  John D. Lazarz  Shock and Detonation Physics, Los Alamos National
Laboratory, Los Alamos, NM, 87545, USA  C. A. Bolme  Shock and Detonation
Physics, Los Alamos National Laboratory, Los Alamos, NM, 87545, USA  Shawn D.
McGrane  Shock and Detonation Physics, Los Alamos National Laboratory, Los
Alamos, NM, 87545, USA  David S. Moore  Shock and Detonation Physics, Los
Alamos National Laboratory, Los Alamos, NM, 87545, USA  M. J. Cawkwell
Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545,
USA  Kyle J. Ramos  High Explosives Science and Technology, Los Alamos
National Laboratory, Los Alamos, NM, 87545, USA

###### Abstract

We use reverse non-equilibrium molecular dynamics ($R$NEMD) simulations to
determine the thermal conductivity in $\alpha$-RDX in the $<$100$>$,
$<$010$>$, and $<$001$>$ crystallographic directions. Simulations are carried
out with the Smith-Bharadwaj non-reactive empirical interatomic potential
[Smith & Bharadwaj, J. Phys. Chem. B 103, 3570 (1999)], which represents the
thermo-elastic properties of RDX with good accuracy. As an illustration, we
report the temperature and pressure dependence of lattice constants of
$\alpha$-RDX, which compare well with experimental and $ab\ initio$ results,
as do linear and volume thermal expansion coefficients, which we also
calculate. We find that the thermal conductivity depends linearly on the
inverse temperature in the 200-400 K regime due to the decrease in the phonon
mean free path. The thermal conductivity also exhibits anisotropy, with a
maximum difference at 300 K of 24% between the $<$001$>$ and $<$010$>$
directions, an effect that remains when temperature increases. Thermal
conductivity in the $<$100$>$ direction is mostly between the two other
directions, although crossovers are predicted with $<$001$>$ at high
temperature, and $<$010$>$ at low temperature under pressure. We observe that
the thermal conductivity varies linearly with pressure up to 4 GPa. The data
are fitted to analytical functions for interpolation/extrapolation and use in
continuum simulations. MD results are validated against experiments using
impulsive stimulated thermal scattering (ISTS) on RDX single crystals at 293 K
and ambient pressure, showing good qualitative and quantitative agreement:
same ordering between the three principal orientations, and an average error
of 10% between the experiments and the model. These results provide confidence
that the extracted analytical functions using the $R$NEMD methodology and the
Smith-Bharadwaj potential can be applied to model the thermal conductivity of
$\alpha$-RDX.

††preprint: AIP/123-QED

##  I Introduction

The accidental initiation of an explosive material subject to an insult is a
major concern to both manufacturers and users. Bowden and Yoffe (2009) In
order to mitigate this type of risk, substantial effort is devoted to
predicting the response of high explosives (HE) to weak shock stimuli, efforts
which require accurate models that notably account for the role of the
material’s microstructure. Indeed, most explosives are characterized by a non-
uniform microstructure, where internal interfaces play critical roles in the
reaction process. Storm, Stine, and Kramer (1990) While the length scale
associated with these HE characteristics (ten nm to a few $\mu$m) invites the
use of mesoscale models, those require accurate parametrizations that are
usually obtained from smaller scale simulations (i.e. atomistic).

One of the key inputs to a mesoscale model is the thermal conductivity of the
material, because heat localization is a necessary precursor to chemical
reactions. Bowden and Yoffe (2009) Anisotropy, pressure, and temperature, are
effects that all must be introduced in order to yield meaningful predictions
of the material’s behavior. Experimentally, McGuire and Tarver (1981);
Fedoroff and Clift (1960); Rogers (1975); Loftus and Gross (1959); Dobratz and
Crawford (1985); Zinn and Mader (1960); Miller (1997); Faubion (1976);
Shoemaker, Stark, and Taylor (1985); Hanson-Parr and Parr (1999); Lawless  _et
al._ (2020) most studies of the thermal conductivity of explosives used
polycrystalline and sometimes porous samples, such that detailed information
on orientation dependences are not readily available. However, since bulk
thermal transport is an atomistic characteristic, Kittel (1966) molecular
dynamics (MD) simulations are well suited to provide insights into the thermal
properties that inform mesoscale models. Classical MD with empirical force
fields was for instance used on two HE crystals, to highlight the anisotropy
and the impact of defects on the thermal conductivity of
2,4,6-triamino-1,3,5-trinitrobenzene (TATB), Kroonblawd and Sewell (2016a, b)
and the thermal conductivity in liquid and crystalline
1,3,5,7-tetranitro-1,3,5,7-tetrazoctane (HMX). Bedrov, Smith, and Sewell
(2000); Chitsazi  _et al._ (2020)

In this work, we present results from MD calculations that predict the thermal
conductivity, $\kappa$, of the heavily used secondary explosive
1,3,5-trinitro-1,3,5-triazinane (RDX). The effects of the crystal orientation,
temperature, and pressure were determined by using the $reverse$ non-
equilibrium molecular dynamics ($R$NEMD) technique, Müller-Plathe (1997);
Zhang  _et al._ (2005) with an accurate non-reactive empirical interatomic
potential that has been shown to reproduce well the thermomechanical response
of RDX. Munday  _et al._ (2011); Josyula, Rahul, and De (2014); Hooks  _et
al._ (2015) We focused on the $\alpha$ polymorph of RDX, which is stable under
ambient conditions and up to $\sim$4 GPa, above which the $\gamma$ phase is
favored.Olinger, Roof, and Cady (1978); Davidson  _et al._ (2008); Dreger
(2012); Cawkwell  _et al._ (2016) We find mostly linear responses to pressure
and the inverse temperature in the regime considered here, corresponding to
moderate insults and operating conditions (200–400 K, 0–4 GPa). We also find
that $\alpha$–RDX displays a notable anisotropy between the $<$100$>$,
$<$001$>$, and $<$001$>$ directions (24% between $<$001$>$ and $<$010$>$ at
300 K and 0 GPa). Analytical expressions for $\kappa$(T,P) are extracted for
each orientation, allowing for the full determination of the crystal thermal
conductivity. Additionally, we show results from impulsive stimulated thermal
scattering (ISTS) to validate the $R$NEMD results. This is a critical
contribution, since, as mentioned above, most experimental results use
polycrystalline samples with varying porosity, which prevents the direct
comparison of MD simulations and experiments. The results, which were obtained
using oriented single crystals, agree quantitatively with the $R$NEMD.

The paper is organized as follows: computational details about the empirical
potential used in the study, the $R$NEMD technique, and simulation geometry
preparation, and the ISTS methodology and setup are given in section II; the
lattice constants and linear and volume thermal expansion coefficients
obtained from the MD simulations are discussed in section III.1; the effect of
the temperature on the thermal conductivity is discussed in section III.2,
while the effect of pressure on the thermal conductivity is presented in the
next section (III.3). In section III.4, we provide a model and coefficients
for $\kappa(T,P)$ for the three principal directions of the crystal. ISTS
results are presented in section III.5. Discussion (sec. IV) and conclusions
(sec. V) follow.

##  II Methods

###  II.1 Smith-Bharadwaj empirical potential

The MD simulations are performed with the Smith-Bharadwaj (SB) empirical
potential (SB-FF) for nitramine explosives that was parameterized to high
level quantum chemical calculations. Smith and Bharadwaj (1999) It is a non-
reactive force field, where the total energy is the sum of five terms:

  * •

bond stretches: $U_{1}=\frac{1}{2}k^{s}_{ab}(r_{ij}-r_{ij}^{0})^{2}$;

  * •

valence bends:
$U_{2}=\frac{1}{2}k^{b}_{abc}(\theta_{ijk}-\theta_{ij}^{0})^{2}$;

  * •

torsions:
$U_{3}=\frac{1}{2}k^{t}_{abcd}\left[1-\textrm{cos}(n\phi_{ijk})\right]$;

  * •

out-of-plane bends: $U_{4}=\frac{1}{2}k^{d}_{abcd}\delta_{ijkl}^{2}$;

  * •

nonbonded:
$U_{5}=A_{ab}\times\textrm{exp}(-B_{ab}r_{ij})-\frac{C_{ab}}{r_{ij}^{6}}+\frac{q_{i}q_{j}}{4\pi\epsilon_{0}r_{rij}}$;

where $k^{s}$, $k^{b}$, $k^{t}$, and $k^{d}$ are force constants, $r$,
interatomic distances, $\theta$, $\phi$, and $\delta$ are angles (dihedral and
out-of-plane for the latter two), and $A$, $B$ and $C$ are constants for the
Buckingham-type non-bonded interactions; indices $a,b,c,d$ refer to the atom
type of atoms $i,j,k,l$, respectively. For the value of each parameter, the
reader is referred to the original paper by Smith and Bharadwaj. Smith and
Bharadwaj (1999) Electrostatic interactions are calculated with the modified
set of partial charges $q_{i}$ proposed by Bedrov et al.. Bedrov  _et al._
(2001)

The force field was successfully applied to reproduce the unit cell
parameters, coefficients of thermal expansion, and heat of sublimation of
several HMX polymorphs, Bedrov  _et al._ (2001) and was also used to predict
the thermal conductivity of HMX. Bedrov, Smith, and Sewell (2000); Chitsazi
_et al._ (2020) The SB-FF was also applied to RDX, and shown to reproduce
lattice and elastic constants, and the equation of state of the $\alpha$ and
$\gamma$ polymorphs. Munday  _et al._ (2011); Josyula, Rahul, and De (2014);
Hooks  _et al._ (2015) Beyond this, studies focused on the behavior of
$\alpha$-RDX under indentation, Weingarten and Sausa (2015) the mechanisms of
the $\alpha\rightarrow\gamma$ phase transformation, Josyula, Rahul, and De
(2019) and the response of RDX to shock compression, Cawkwell  _et al._
(2010); Ramos  _et al._ (2010); Bidault and Pineau (2018) were also performed
with the SB-FF. The good performance of the SB-FF for RDX is not too
surprising, considering its accuracy for HMX, for which it was designed, and
the similarities between the two nitramine molecules.

###  II.2 $R$NEMD simulations

The thermal conductivity $\kappa$ of a material is defined by Fourier’s law
as:

$\kappa=\frac{\bm{J}}{\nabla T},$ (1)

where $\bm{J}$ is the heat flux response to a temperature gradient $\nabla T$,
and vice versa. Following, there are two main paths to calculate $\kappa$ from
MD simulations:

1) $Impose$ $\nabla T$, $measure$ $\bm{J}$
This is similar to an experimental protocol, and is commonly referred to as
the “direct” method; the temperature gradient is imposed by thermostating
distant regions of the sample at different temperatures (namely, $T_{hot}$ and
$T_{cold}$), and measuring the response heat flux $\bm{J}$; the latter can be
done by monitoring the thermostats and extracting the amount of energy
added/subtracted from the respective regions. While the direct method is
intuitive and straightforward to implement in simulations, it suffers from
serious practical difficulties, as was discussed in Refs. Müller-Plathe, 1997
and Schelling, Phillpot, and Keblinski, 2002: $\bm{J}$ is a quantity with
large fluctuations, such that large temperature gradients must be used in
order to distinguish trend from noise. These temperature gradients are also
usually out of the range of experimental measurements, such that the results
might not be directly comparable to experiments. Additionally, a large $\nabla
T$ leads to non-linearities near the thermostated regions, which violates the
linear response assumption governing Eq. 1.

2) $Impose$ $\bm{J}$, $measure$ $\nabla T$
An alternative to the direct method was proposed by Müller-Plathe: Müller-
Plathe (1997); Zhang  _et al._ (2005) since $\bm{J}$ is the problematic
quantity, one can, by contrast with the direct method, impose the flux, which
is therefore known exactly. Further, the temperature gradient $\nabla T$
becomes the response quantity. Since $\nabla T$ is measured over a reasonable
number of particles, it results in much better statistical averages. The heat
flux is obtained by periodically swapping molecular center-of-mass velocities
of “hot” and “cold” molecules between “cold” and “hot” regions, in order to
make the “cold” region colder, and the “hot” region hotter. This is known as
the $reverse$ non-equilibrium molecular dynamics ($R$NEMD), or Müller-Plathe,
technique, and is the method employed in this paper. The $R$NEMD method has
been validated against both the direct method and the more computationally
extensive and explicit Green-Kubo method. Müller-Plathe (1997); Schelling,
Phillpot, and Keblinski (2002); Izvekov, Chung, and Rice (2011)

Figure 1: $R$NEMD scheme. Here, the sample consists of 12$\times$3$\times$3
unit cells oriented in the $<$100$>$ direction. Each layer consist of
1$\times$3$\times$3 unit cells and PBC are applied on all directions, such
that the cold layer is repeated in the figure. Velocity swaps occur between
the “cold” and “hot” regions, producing a heat flux $\bm{J}$ away from the hot
region and a resulting temperature gradient $\nabla T$ that is extracted from
the temperature profile. In this particular case ($T$=300 K, $P$=0 GPa), we
obtained $T(z)=390.4-2.13*|z|$ and $T(z)=390.6-2.14*|z|$ with the origin in
the middle of the hot “region”.

In detail, the $R$NEMD implementation is as follow (see also Fig. 1):

  1. 1.

Equilibrate the system at the target average temperature, using for instance a
Langevin thermostat (NVT simulations).

  2. 2.

Define hot and cold regions and intermediate layers, and identify which
molecules belong to each region.

  3. 3.

Remove the thermostat and evolve the whole system in the constant energy
ensemble (NVE) for a time $t$.

  4. 4.

Identify the molecule in the hot region with the smallest kinetic energy, and
the molecule in the cold region with the largest kinetic energy.

  5. 5.

Exchange the atomic velocities of the two molecules.

  6. 6.

Repeat N times from step 3.

The MD simulations (steps 1 and 3) were performed with the LAMMPS code;
Plimpton (1995) while the identification and swapping procedures, as well as
the processing of inputs/outputs were performed with in-house tools. At each
swap of velocities (step 5), the exact amount of heat exchanged is known, and
the flux along the direction $\alpha$ is given by:

$\bm{J}_{\alpha}=\frac{\sum_{\textrm{transfers}}\frac{m}{2}(v_{h}^{2}-v_{c}^{2})}{2\times
A\times t},$ (2)

where $m$ is the mass of the molecule, $v_{h}$ and $v_{c}$ the center-of-mass
velocities of the hot and cold molecules, respectively, $A$ the cross-
sectional area of the sample, and t the total time of the simulation. The
factor of 2 accounts for the two directions of the heat transfer (left and
right, see Fig. 1). Thermal equilibrium is reached after a number of cycles,
see Fig. 2. The temperature gradient is then calculated by performing a linear
fit of the temperature as a function of the distance from the hot region (Fig.
1). The average temperature in each layer is obtained according to
$T=\nicefrac{{2E_{\textrm{kin}}}}{{3Nk_{B}}}$, with $E_{\textrm{kin}}$ the
total kinetic energy in the layer, $N$ the number of atoms in the layer, and
$k_{B}$ the Boltzmann constant.

From Fig. 2, we notice that the total energy exhibits a nearly constant but
very small drift during the course of the simulation. The origin of the drift
is unknown. Simple NVE simulations with an identical system showed similar
behavior, such that the swapping procedure is not at fault; the effect also
remains when a smaller timestep (0.25 vs. 0.4 fs) is used. In the supplemental
material to Ref. Mathew  _et al._ , 2018, a study using a SB-type FF
parameterized for TATB (2,4,6-triamino-1,3,5-trinitrobenzene), Mathew et al.
mention a discontinuity of the potential energy function when the three-center
angles approach 180∘. We wonder if a similar effect might be at play for RDX.
The effect being very small ($<2\times 10^{-9}$ eV/molecule/timestep), its
impact could also have been overlooked, and in any case the issues was not
investigated further and was deemed as irrelevant in the current study.

Armed with the heat flux determined from Eq. 2, and the thermal gradient
extracted from the temperature profiles (Fig. 1), the thermal conductivity of
the sample is then directly determined from Eq. 1. The last step is to take
into account finite size effects that lead to scattering, as was mentioned in
Ref. Schelling, Phillpot, and Keblinski, 2002. Therein, the authors discuss
the specific task of simulating thermal conductivity in a solid under periodic
boundary conditions, and show that the size effects lead to a linear
dependence of the resistivity (1/$\kappa$) with respect to the inverse of the
sample length $L_{\alpha}$:

$\frac{1}{\kappa}=\frac{1}{\kappa_{\infty}}-\frac{\lambda}{L_{\alpha}},$ (3)

where $\kappa_{\infty}$ is the thermal conductivity extrapolated to an
infinite sample (the value we report in this paper), and $\lambda$ a
coefficient. For each case (orientation, T, P), we run simulations with
samples of four different lengths, and the thermal conductivity is obtained by
fitting these four data points with Eq. 3. Each simulation consists of 250000
timesteps of equilibration (100 ps), followed by 20000 swaps (5 ns), of which
the last 10000 are used for production.

Figure 2: Profiles from the $R$NEMD simulation for $\alpha$-RDX at 300 K, 0
GPa, in the $<$100$>$ orientation (12$\times$3$\times$3 unit cells). Top:
temperature profiles for the hot region (red), cold region (blue), and sample
(black). Middle: total energy profile. Bottom: kinetic energy exchanged at
each swap. The total simulation time of 5 ns corresponds to 20000 swaps.

###  II.3 MD sample preparation

$\alpha$-RDX is an orthorhombic crystal with space group $Pbca$ and 8
molecules (168 atoms) per unit cell McCrone (1950); Choi and Prince (1972). In
order to determine the equilibrium sample size for given P/T initial
conditions, we first perform simulations with a system of 3$\times$3$\times$3
unit cells in the NPT (isothermal–isobaric) ensemble. Pressure and temperature
are imposed via Nose-Hoover barostat/thermostat, with damping parameters of
100 and 1000 timesteps, respectively; the orthorhombic symmetry is imposed by
fixing the three angles of the supercell at 90 degrees, but the three
dimensions are allowed to vary independently. A timestep of 0.5 fs is used for
the NPT simulations, and the lattice constants are averaged over the last 1 ns
of a 2 ns run. In order to further validate the SB-FF, we also determined the
lattice constants of $\gamma$-RDX (orthorhombic $Pca2_{1}$, 8 molecules/unit
cell) Davidson  _et al._ (2008) as a function of pressure, using a
3$\times$3$\times$3 supercell.

For the thermal conductivity calculations, we built samples of 8, 12, 16, and
24 unit cells in length and a cross section set at 3$\times$3 unit cells. In
order to investigate the anisotropy of the thermal conductivity, the
directions along the three perpendicular lattice vectors of the crystal are
considered ($<$100$>$, $<$010$>$, and $<$001$>$). These specific crystal
orientations were chosen since they produce an orthorhombic configuration that
allows for PBC simulations without the use of specific techniques to build
infinite samples in arbitrary directions (see Ref. Kroonblawd  _et al._ ,
2016). The resulting samples contain between 12096 and 36288 atoms, with
dimensions at 300 K and 0 GPa between 10.7 and 32.3 nm for samples oriented
along $<$100$>$ (3.5$\times$3.2 nm2 cross-section), 9.2 and 27.7 nm
(4.0$\times$3.2 nm2 cross-section) for $<$010$>$, and 8.4 and 25.3 nm for
$<$001$>$ (4.0$\times$3.5 nm2 cross-section). The size of the smallest
longitudinal dimension is still near twice the estimated phonon mean-free path
in $\beta$-HMX, Chitsazi  _et al._ (2020) which we use as rough estimate for
RDX as well.

Before running the $R$NEMD simulations, the samples with the appropriate
lattice constant were equilibrated at the target temperature for 100 ps with a
Langevin thermostat. The NVE runs between velocity swaps had a duration of 250
fs, in accordance with previous studies; Kroonblawd and Sewell (2016a) we also
independently checked on shorter runs (5000 swaps of equilibration followed by
5000 swaps for production) that the NVE run duration did not impact the
results, along with the cross-section size, see Table 1. Apart from the
smallest cross-section, thermal conductivities lie within $\sim$10% of each
other. The 3$\times$3 cross-section / 250 fs swap frequency combination was
thus used in this work. A timestep of 0.4 fs is used for the $R$NEMD
simulations.

Table 1: Sensitivity of the thermal conductivity $\kappa$ to the swap frequency and cross-section of the sample. System: $<$100$>$ $\alpha$-RDX, T=300 K, P=0 GPa, length=8 unit cells. cross-section | swap frequency (fs) |  $\kappa$ (W.m-1.K-1)
---|---|---
2$\times$2 | 250 | 0.409$\pm$0.010
3$\times$3 | 250 | 0.237$\pm$0.003
4$\times$4 | 250 | 0.261$\pm$0.004
3$\times$3 | 500 | 0.248$\pm$0.003
3$\times$3 | 750 | 0.232$\pm$0.002

###  II.4 Experimental

Impulsively stimulated thermal scattering (ISTS) was performed on $<$100$>$,
$<$010$>$, and $<$001$>$ oriented RDX crystals. RDX crystals were grown from
acetone solution after removal of contaminants and imperfections from the
starting material by recrystallization. RDX crystals were then indexed, cut,
and polished for optical clarity for use in transmission in ISTS. The
resulting slab geometry RDX crystals of approximately 7$\times$7$\times$0.7-1
mm3 size were affixed to a rotation stage and oriented perpendicularly to the
probe beam. A set of translation stages was used to position each crystal at
the optimal probe volume for ISTS.

ISTS is a non-contact optical technique for measuring thermal conductivity and
sound speed in-situ. Abramson, Brown, and Slutsky (1999); Dennett and Short
(2018); Johnson  _et al._ (2012); Käding  _et al._ (1995); Maznev, Nelson, and
Rogers (1998); Rogers, Maznev, and Nelson (2012); Rogers  _et al._ (2000);
Rogers, Yang, and Nelson (1994); Tokmakoff, Banholzer, and Fayer (1993) A pair
of 28 ps pump pulses (500 Hz repetition rate) at 1064 nm and continuous wave
532 probe beams were generated using the $\pm$1 order of a diffractive optical
element. A spatial mask removed extraneous orders for both beams. The pump
beams and probe were collimated and then focused into the sample. At the
sample the pump and probe beams were approximately 100 $\mu$m and 50 $\mu$m in
diameter respectively. Pump beams were crossed at $\sim$4.8∘ in the RDX
sample. A small amount of the pump beam energy was absorbed by the crystal
generating a thermally induced grating. The thermal excitation launched
counter propagating pressure waves that travel at the sample acoustic
velocity. The magnitude of the thermal grating decreases via conduction to
ambient conditions. The probe beam was diffracted by the generated grating.
The diffracted beam was sent to a DC-600 MHz silicon amplified detector
(Thorlabs FPD610-FC-VIS) to measure the acoustic velocity and thermal grating
decay. A 1 GHz oscilloscope was used to record the probed signal. Since RDX
has weak absorption at the pump wavelength, the diffracted signal was small
enough that a heterodyne technique was employed, wherein a spatially
coincident local oscillator (LO) was mixed with the diffracted signal to
simultaneously increase weak diffraction signals and reduce parasitic
oscillator contributions. Maznev, Nelson, and Rogers (1998); Rogers, Maznev,
and Nelson (2012); Rogers  _et al._ (2000) A calcite phase plate was used to
control the relative phase between the signal and LO beams. In phase and out
of phase were defined to be 0∘ and 180∘, respectively, relative between the LO
and signal beam phases. 512 shots were averaged per phase point during data
collection. Thermal conductivity and sound speed data were collected at 5∘
radial increments through a full revolution. Shown in Fig. 3 is a schematic of
the experimental apparatus. The LO and diffracted signal were spatially masked
to reduce contributions from light scattering in the sample. Preliminary
discussion of ISTS in $\alpha$-RDX was provided in Ref. Lazarz  _et al._ ,
2020, although all results presented here were obtained from new measurements.

Figure 3: Schematic of the experimental apparatus. A box-CARS geometry was
used to spatially separate the signal and pump beams after transmission
through the sample. A 1.0 neutral density filter was placed in the local
oscillator path. Due to the Bragg condition the local oscillator and
diffracted probe are coincident after the sample.

##  III Results

###  III.1 Lattice constant and thermal expansion coefficient of $\alpha$-RDX
as a function of pressure and temperature

Table 2: Comparison of the lattice constants and volume of the $\alpha$ and $\gamma$-RDX unit cells obtained in this work with the SB-FF to experimental and computational values from the literature. Error percentages with respect to experimental values are in parentheses. $\alpha$-RDX at 300 K and 0 GPa |  $\gamma$-RDX at 300 K and 5.2 GPa
---|---
parameter | this work | DFT-D [Sorescu and Rice, 2010] | DFT-D [Hunter  _et al._ , 2013] | Exp [Olinger, Roof, and Cady, 1978] |  | this work | SB-FF [Josyula, Rahul, and De, 2014] | DFT-D [Sorescu and Rice, 2010] | Exp [Davidson  _et al._ , 2008]
a(Å) | 13.486 (+2.3) | 13.237 (+0.4) | 13.282 (+0.8) | 13.182 |  | 12.718 (+1.2) | 12.71 (+1.2) | 12.699 (+1.1) | 12.565
b(Å) | 11.546 (-0.2) | 11.391 (-1.6) | 11.419 (-1.3) | 11.574 |  | 11.057 (+1.2) | 11.05 (+1.1) | 10.918 (-0.1) | 10.930
c(Å) | 10.552 (-1.5) | 10.770 (+0.6) | 10.736 (+0.3) | 10.709 |  | 9.657 (+1.9) | 9.64 (+1.7) | 9.503 (+0.3) | 9.477
V(Å3) | 1642.93 (+0.6) | 1623.94 (-0.6) | 1628.27 (-0.3) | 1633.86 |  | 1357.91 (+4.3) | 1354.9 (+4.1) | 1317.66 (+1.2) | 1301.5
Figure 4: a: Lattice constants of $\alpha$ and $\gamma$-RDX as a function of
pressure at 300 K, compared to experimental results from Refs. Olinger, Roof,
and Cady, 1978; Davidson  _et al._ , 2008 and zero temperature DFT+D results
from Ref. Sorescu and Rice, 2010. b: Lattice constants of $\alpha$-RDX as a
function of temperature at 0 GPa from NPT-MD simulations, compared to
experimental results from Bolotina and Pinkerton Bolotina and Pinkerton (2015)
and Sun et al.. Sun  _et al._ (2011) c: Lattice constants of $\alpha$-RDX as a
function of temperature and pressure from NPT-MD simulations. The cusps
observed at P$\geq$3 GPa correspond to transition to the $\gamma$-phase during
the MD simulations. d: Coefficient of thermal expansion in $\alpha$-RDX as a
function of temperature, at 0 GPa. Results from this work are compared to
experimental results from Cady Cady (1972) and Bolotina & Pinkerton. Bolotina
and Pinkerton (2015)

We first compare the equilibrium lattice constant and unit cell volume of
$\alpha$ and $\gamma$-RDX (at 5.2 GPa for the latter), and compare it to
experimental and $ab\ initio$ results, see Table 2. Included are results from
DFT+D calculations (DFT with a correction for dispersion forces to reproduce
crystalline equilibrium lattice constant). It can be seen that the SB-FF
reproduces the lattice constants and unit cell volume with good accuracy, with
an error compared to experiments of 0.6% for the volume of the $\alpha$-RDX
unit cell, and a maximum error on the lattice constant of 2.3% (1.6% for
DFT+D). For $\gamma$-RDX, the error is below 5% for the volume of the unit
cell (1.2% for DFT+D), and 1.9% at most for the lattice constants. Our results
also agree well with previous calculations performed with the SB-FF for
$\gamma$-RDX.

Figure 4-a shows the lattice constants of $\alpha$ and $\gamma$-RDX as a
function of pressure, at 300 K. Results obtained with the SB-FF (this work)
are compared to experiments Olinger, Roof, and Cady (1978); Davidson  _et al._
(2008) and DFT+D calculations. Sorescu and Rice (2010) We again notice the
good performance of the SB-FF, with errors on the order of the DFT+D errors.
In addition, the qualitative changes in $a,b,c$ lattice constants between the
$\alpha$ and $\gamma$ phases are reproduced correctly: $a$ shows almost no
change, $b$ increases, $c$ decreases. Figure 4-b reports the temperature
dependence of the lattice constants and volume in $\alpha$-RDX, at zero
pressure, and compared to experimental results using X-ray diffraction data.
Sun  _et al._ (2011); Bolotina and Pinkerton (2015) The SB-FF reproduces
qualitatively and quantitatively the experimental temperature dependence of
the lattice constants and volume, with errors below 3% for the lattice
constant and excellent agreement on the volume, albeit through a compensation
of errors on the $a$ and $c$ lattice constants. The increase between 50 and
500 K is 1.5% for $a$ (1.1% from experiments), 2.8% for $b$ (3.8%), and 3.7%
for $c$ (3.4%). For the volume, our work predicts a total increase over the
temperature range of 8.2%, while Bolotina & Pinkerton report 8.4%. Bolotina
and Pinkerton (2015)

Lastly for the lattice constants, Fig. 4-c reports the variation of $a$, $b$,
and $c$ as a function of pressure between 250 and 500 K, in $\alpha$-RDX. We
see that the pressure hinders the lattice expansion, with a marked effect on
$c$ especially: the variation between the 250 and 500 K lattice constants goes
from 1.0, 1.6, 2.5% at 0 GPa to 0.4, 1.0, and 0.7% at 2GPa, for $a$, $b$, and
$c$ respectively. Additionally, we note a cusp for $b$ and $c$ above 325 K,
starting at 3 GPa. This marks the phase transition to the $\gamma$ phase,
with, accordingly, an increase in $b$ and a decrease in $c$, while $a$ remains
about the same. While the experimental phase transition is reported closer to
4 GPa, Davidson  _et al._ (2008) the procedure employed here does not
guarantee that 3 GPa is the “true” transition pressure with the SB-FF. Small
systems used in MD simulations lead to artificial constraint through periodic
boundary conditions, and the simulation timescale is orders of magnitude
smaller than in experiments. What can be said with confidence is that
simulations at 3 GPa and above should be closely monitored for a possible
phase transition.

Table 3: Analytic expressions for the lattice constants and thermal expansion
coefficients in $\alpha$-RDX as a function of pressure and temperature,
compared to available experimental results. Lattice constants are given by
L(T)=$a_{0}+a_{1}\times 10^{-4}\times T+a_{2}\times 10^{-7}\times T^{2}$,
volume by V(T)=$a_{0}+a_{1}\times 10^{-1}\times T+a_{2}\times 10^{-4}\times
T^{2}$, thermal expansion coefficients by $\alpha$(T)=$a_{0}\times
10^{-5}+a_{1}\times 10^{-7}\times T$. Last column denotes the temperature
range over which the data was fitted.

Lattice constants and volume
---
| 100 | 010 | 001 | V | T
| $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{0}$ | $a_{1}$ | $a_{2}$ |
Ref. Bolotina and Pinkerton,2015 | 13.1290 | 1.316 | 3.400 | 11.3557 | 7.236 | 4.163 | 10.5640 | 2.200 | 10.69 | 1575.14 | 1.463 | 2.745 | 90–300 K
P=0 GPa | 13.3853 | 2.126 | 4.137 | 11.3440 | 6.363 | 1.146 | 10.3758 | 3.100 | 9.373 | 1575.75 | 1.563 | 2.280 | 50–500 K
P=1 GPa | 13.1679 | 1.837 | 1.049 | 11.1796 | 5.147 | 0.302 | 10.1956 | 2.896 | 2.279 | 1501.08 | 1.308 | 0.569 | 250–500 K
P=2 GPa | 13.0075 | 1.738 | 0.520 | 11.0792 | 3.490 | 1.315 | 10.0586 | 3.030 | -0.246 | 1449.63 | 1.083 | 0.240 | 250–500 K
P=3 GPa | 12.8803 | 1.046 | 0.200 | 10.9791 | 3.745 | 1.352 | 9.9506 | 3.296 | -1.600 | 1407.16 | 1.061 | -0.035 | 250–325 K
P=4 GPa | 12.7500 | 1.894 | 0.400 | 10.9255 | 1.774 | 2.000 | 9.8428 | 4.597 | -3.760 | 1371.09 | 1.0681 | -0.021 | 250–325 K
P=5 GPa | 12.6580 | 1.490 | 0.514 | 10.8760 | 0.544 | 2.819 | 9.7805 | 3.317 | -1.291 | 1346.50 | 0.678 | 0.025 | 250–400 K

Thermal expansion coefficients
---
| 100 |  | 010 |  | 001 |  | V |  | T
| $a_{0}$ | $a_{1}$ |  | $a_{0}$ | $a_{1}$ |  | $a_{0}$ | $a_{1}$ |  | $a_{0}$ | $a_{1}$ |  |
Ref. Bolotina and Pinkerton,2015 | 1.007 | 0.512 |  | 6.399 | 0.663 |  | 2.134 | 1.97 |  | 9.539 | 3.14 |  | 90–300 K
Ref. Sun  _et al._ ,2011 | 3.07 | – |  | 8.28 | – |  | 9.19 | – |  | 20.7 | – |  | 303–443 K111Fitted as constant
Ref. Cady,1972 | 2.439 | 0.842 |  | 8.531 | 0.673 |  | 7.391 | 2.073 |  | 18.33 | 3.625 |  | 113–408 K222Terms in $T^{2}$ and beyond ignored, $<$100$>$ fitted up to 423 K.
P=0 GPa | 1.604 | 0.603 |  | 5.618 | 0.163 |  | 3.089 | 1.717 |  | 10.30 | 2.478 |  | 50–500 K
P=1 GPa | 1.400 | 0.154 |  | 4.608 | 0.031 |  | 2.848 | 0.421 |  | 8.854 | 0.606 |  | 250–500 K
P=2 GPa | 1.339 | 0.082 |  | 3.167 | 0.225 |  | 3.009 | -0.056 |  | 7.514 | 0.250 |  | 250–500 K
P=3 GPa | 0.818 | 0.307 |  | 3.408 | -0.082 |  | 3.302 | -0.325 |  | 7.530 | -0.101 |  | 250–325 K
P=4 GPa | 1.487 | 0.060 |  | 1.634 | -0.357 |  | 4.641 | -0.762 |  | 7.761 | -0.346 |  | 250–325 K
P=5 GPa | 1.179 | 0.079 |  | 0.513 | 0.511 |  | 3.380 | -0.268 |  | 5.067 | 0.323 |  | 250–400 K

From these results, we extracted analytic functions to predict the lattice
constants and volume as a function of pressure and temperature, which are
reported in Table 3. The quadratic expression also leads to a linear relation
for the coefficient of linear thermal expansion in the direction $L$,
$\alpha_{L}$(T):

$\alpha_{L}(T)=\frac{1}{L}\frac{\textrm{d}L(T)}{\textrm{d}T},$ (4)

and similarly for the coefficient of volume thermal expansion $\alpha_{V}$(T).
The thermal expansion coefficients are also reported in Table 3, and compare
well with recent experimental results from Ref. Bolotina and Pinkerton, 2015.
The thermal expansion coefficients are plotted in Fig. 4-d, along with
experimental results from Refs. Cady, 1972; Bolotina and Pinkerton, 2015.

###  III.2 Thermal conductivity in $\alpha$-RDX as a function of temperature

Figure 5: Thermal resistivity (1/$\kappa$) as a function of the inverse sample
length (1/$L$) and temperature for $\alpha$-RDX in the $<$100$>$, $<$010$>$,
and $<$001$>$ orientations, at 0, 1, and 2 GPa. The solid lines are fit to the
data using Eq. 3.

In Fig. 5, we show the dependence of the thermal resistivity (1/$\kappa$)
versus the inverse of the sample length (1/$L$), which is predicted to be
linear according to Eq. 3. For each point, the thermal conductivity
$\kappa(L,T,P)$ was obtained via Eq. 1. Heat exchanges were summed over 1000
swaps, and the sums averaged over 10000 swaps, which gives us the average heat
flux $<$$\bm{J}$$>$ and the associated error $\delta(\bm{J})=\sigma_{\bm{J}}$
($\sigma_{\bm{J}}$ is the standard deviation). The temperature gradient
$\nabla(T)$ is obtained by fitting the temperature profiles, which produces
two values since there are two directions of propagation (see Fig. 1); the
standard deviation between the two values is used as error on $\nabla(T)$:
$\delta(\nabla T)=\sigma_{\nabla T}$ (the standard error on the fit itself is
negligible). Using standard error propagation rules Bevington (1969) and
neglecting the covariance between $\nabla T$ and $\bm{J}$, the error
$\delta(\kappa)$ on the thermal conductivity $\kappa$ is obtained as:

$\delta(\kappa)=\sqrt{\kappa^{2}\left(\frac{\delta(\bm{J})^{2}}{<\bm{J}>^{2}}+\frac{\delta(\nabla
T)^{2}}{<\nabla T>^{2}}\right)},$ (5)

and the error on the thermal resistivity $1/\kappa$ is then simply
$\delta(1/\kappa)=\delta(\kappa)/\kappa^{2}$. All fits were performed with the
R software package, R Core Team (2017) and each point was weighted by a factor
equal to the inverse of the standard deviation of the measurement, when
available, to account for the statistical uncertainty during fitting.

The thermal resistivity data in Fig. 5 follows mostly a linear trend, which
validates the behavior predicted by Eq. 3, and allows us to extract the
extrapolated size-independent thermal conductivity, $\kappa_{\infty}$, which
is the quantity most relevant to higher scale models and to compare to
experiments. The results are reported in Figs. 6-a–c for the temperature
dependence.

The first important information that can be obtained from Figs. 6-a–c is that
$\alpha$-RDX exhibits anisotropic thermal conductivity, with
$\kappa_{001}>\kappa_{100}>\kappa_{010}$, and a difference of about 20%
between $\kappa_{001}$ and $\kappa_{010}$ at 300 K. Previous results obtained
with the SB-FF and a slightly different $R$NEMD method by Izvekov et al.
predict the same ordering for $\kappa_{001}$, $\kappa_{010}$, and
$\kappa_{001}$, and a difference of 24.6% between the two extremes at 300 K
(more discussion in section IV).

It is also clear from Figs. 6-a–c that the thermal conductivity of
$\alpha$-RDX follows $\kappa\propto 1/T$ in the 200–400 K regime (simulations
performed above 400 K occasionally led to a phase transition and were thus
discarded). This is in agreement with theory that predicts a $1/T$ dependence
for a solid without defects at “high” temperature, i.e. $T>>T_{\Theta}$, where
$T_{\Theta}$ is the Debye temperature of the solid. Klemens (1969) Using the
isotropic continuum approximation and the sound speed extracted from the SB-FF
and from experiments for $\alpha$-RDX, Ref. Izvekov, Chung, and Rice, 2011
reports $T_{\Theta}$=100 and 125 K, respectively, while Rey-Lafon & Bonjour
estimate $T_{\Theta}\simeq$60 K from the acoustic phonons at low temperature.
Izvekov, Chung, and Rice (2011); Rey-Lafon and Bonjour (1973) Standard
theories show that the thermal conductivity varies as $\kappa\propto T^{3}$ at
low temperature, however this is solely due to the temperature dependence of
the heat capacity. Kittel (1966) With heat capacity being a constant in
classical MD, there is no contribution to reduce the thermal conductivity at
low temperature, and therefore a monotonic behavior can be expected.
Additional calculations performed in the 50–200 K range for $\alpha$-RDX
$<$100$>$ confirm a monotonous behavior, however the data departs from the
linear behavior starting around 100–150 K (see Fig. 6-a), suggesting that
extrapolation of the trends extracted in the 200–400 K range to lower
temperature should be taken with extreme caution. Interestingly, Izvekov et
al. Izvekov, Chung, and Rice (2011) reported a maximum of the thermal
conductivity for single crystal RDX when using their NEMD method, at around
275 K, for all three directions. While that maximum seems to correlate with
experimental observations compiled by Miller, Miller (1995) this is most
likely fortuitous: Miller reports limited data points with large error bars,
and used polycrystalline samples with density $\rho$=1.716 g/cm3, vs.
$\sim$1.8 g/cm3 for ideal RDX crystals. In addition, a rigorous calculation
with the Green-Kubo method, also performed in Ref. Izvekov, Chung, and Rice,
2011, did not reproduce this result, and the authors acknowledged the absence
of a clear reason why a non-monotonic thermal conductivity was observed in a
perfect crystal within the classical approximation.

Figure 6: Thermal conductivity of $\alpha$-RDX in the $<$100$>$, $<$010$>$,
and $<$001$>$ orientations, as a function of temperature at 0 GPa (a), 1 GPa
(b), and 2 GPa (c), and as a function of pressure at 200 K (d), 300 K (e), and
400 K (f). Lines serve as guide to the eye (the fit to the data is discussed
in Section III.4). (a), insert: $<$100$>$ results down to 50 K, showing the
breakdown of the linear behavior below 100–150 K. Note the log scale on the
$x$ axis. (e), insert: $<$010$>$ and $<$001$>$ results up to 4 GPa, showing
that the linear behavior holds at higher pressure.

###  III.3 Thermal conductivity in $\alpha$-RDX as a function of pressure

Figures 6-d–f reports the thermal conductivity of $\alpha$-RDX between 0 and 2
GPa for each orientation, at $T$=200, 300, and 400 K (additional results from
our calculations performed at 250 and 350 K are omitted for clarity).
Preliminary calculations at 300 K showed that more than one simulation
exhibited signs of a phase transition at 3 and 4 GPa, which prevented us from
performing a reliable linear fit according to Eq. 3. This occurred
specifically for samples in the $<$100$>$ orientation, but also for one case
in the $<$001$>$ orientation, at 3 GPa (in that case the fit was still
performed, based on 3 sample sizes). The onset of the phase transition is
marked by changes in pressure and potential energy, however, since the volume
of the sample is fixed at the beginning of the simulation and constrained
throughout, the suspected transformation to $\gamma$ cannot be completed, and
an unphysical structure is obtained. While the experimental
$\alpha\rightarrow\gamma$ pressure phase transition is closer to 4 GPa, the
SB-FF was shown to lead to phase transition at different pressures during
uniaxial compression (depending on the axis), while transition under
hydrostatic compression was not observed up to 9 GPa. Munday  _et al._ (2011)
In addition to the accuracy of the force field itself, the size of the system
and timescale of the simulation play a role here, since PBC introduce
artificial constraints that can affect phase transition behavior. As a result,
calculations were limited to the 0–2 GPa range, but the 3 and 4 GPa results at
300 K, are still presented in Fig. 6-e to show that the trends still hold.

Our results suggest that $\kappa$ is a linear function of pressure.
Specifically, $\kappa_{010}$ and $\kappa_{001}$ follow similar quantitative
trends, while $\kappa_{100}$ is less affected by pressure. As a result, the
difference between $\kappa_{010}$ and $\kappa_{001}$ remains about the same as
a function of pressure; $\kappa_{100}$ increases the most with respect to the
other two directions, even crossing with $\kappa_{010}$ at 200 K, $\sim$1.5
GPa (the crossover pressure increases beyond 2 GPa at higher temperature,
though the trends are similar). Experimentally, the linear trends should be
measurable beyond uncertainties, with an increase of nearly 30% between 0 and
2 GPa for $\kappa_{001}$ at 200 K, for instance.

While some theoretical models suggest a turnover of $\kappa(P)$ in molecular
crystals at some pressure, Long  _et al._ (2012); Fan  _et al._ (2017) this
behavior is due to phonon softening, which usually occurs during a phase
transition. Nakanishi, Nagasawa, and Murakami (1982) Since our calculations
were limited, by design, to a single phase of RDX, a monotonic relation
between the thermal conductivity and pressure is to be expected.

###  III.4 Analytical expression for $\kappa(T,P)$

Based on the results presented in sections III.2 and III.3, we propose a
simple analytical expression that represents the $T$ and $P$-dependence of the
thermal conductivity,

$\kappa(T,P)=\kappa_{0}+\frac{\alpha}{T}+\beta P$ (6)

where $a$,$b$, and $c$ are coefficients extracted, for each orientation, from
the data obtained for $T$=200–400 K and $P$=0–2 GPa. The parameters are
presented in Table 4, and yield excellent agreement with the data: average
errors between the predicted values and the data are less than 3% for all
orientations. A better agreement ($\leq$1.1% on average for each orientation)
can be obtained by using a more complex form, i.e. a second order function
$\kappa(T,P)=a+b/T+c~{}P+d~{}P/T+f/T^{2}+g~{}P^{2}$; however, the underlying
theory does not support a quadratic dependence on the inverse temperature,
which may lead to drastic deviations when extrapolation is used. A different,
empirical, form was recently employed in Ref. Leiding  _et al._ , 2021 to fit
the average thermal conductivity in $\beta$-HMX, which incidentally seems to
better match previous MD data for liquid HMX outside of the fitting range:
$\kappa=\kappa_{0}(\textrm{exp}[-\alpha T]+\beta P)$. Using this function did
not lead to a significant improvement of the fit in our case and, as a result,
the linear form of Eq. 6 is preferred.

Finally, we used the parameters presented in Table 4 and Eq. 6 to map the
thermal conductivity as a function of pressure and temperature, extrapolated
to $T$=500 K and $P$=5 GPa, i.e. enclosing and exceeding the experimental
domain of stability of $\alpha$-RDX. Dreger (2012) For each orientation, we
observe that $\kappa(T,P)$ decreases monotonically from the high-$P$/low-$T$
region toward the low-$P$/high-$T$ region.

Table 4: Parameters obtained from fitting the $\kappa(T,P)$ data for each orientation, according to Eq. 6. Fits performed to the data at T=200–400 K, P=0–2 GPa. Units are $\kappa_{0}$: W.m-1.K-1; $\alpha$: W.m-1; $\beta$: 10-9m2.s-1.K-1 (=W.m-1.K-1.GPa-1). Standard error on the parameters provided in parentheses. orientation | $\kappa_{0}$ | $\alpha$ | $\beta$
---|---|---|---
$<$100$>$ | 0.228 (0.007) | 29.85 (1.99) | 0.0667 (0.002)
$<$010$>$ | 0.156 (0.013) | 36.75 (3.31) | 0.0906 (0.004)
$<$001$>$ | 0.151 (0.013) | 63.26 (3.98) | 0.0841 (0.003)
Figure 7: Predicted thermal conductivity of $\alpha$-RDX in the $<$100$>$,
$<$010$>$, and $<$001$>$ orientations, as a function of pressure and
temperature up to 5 GPa and 500 K, using Eq. 6 and the parameters presented in
Table 4. Contours are plotted at 0.1 W.m-1.K-1 intervals.

###  III.5 Experimental Results

Figure 8: Left: Raw RDX $<$010$>$ ISTS measurements. The sound speed can be
calculated from the frequency in the insert. The decay of the signal can be
used to calculated the thermal conductivity. Right: Average thermal
conductivities for $<$100$>$, $<$010$>$, and $<$001$>$ RDX crystal
orientations through ISTS. Markers correspond to indexed crystal orientation,
not the probed orientation. Probed directions are labelled at peaks and
troughs.

An RDX $<$010$>$ ISTS measurement is shown in Fig. 8, left. Two time domain
signals are present in the recorded measurement: the 10’s nanosecond fast
damped oscillatory portion and the 20 microsecond slow decay. The damped
oscillatory portion can be used to calculate the sound speed based on the
frequency of the damped oscillations and the wavelength of the probe beam
Rogers, Maznev, and Nelson (2012). Rather than use a power spectrum fit
Rogers, Maznev, and Nelson (2012); Rogers, Yang, and Nelson (1994); Sun  _et
al._ (2008) to find the frequency of the oscillatory portion, a nonlinear fit
based on an analysis of variance was used,

$I=\textrm{Acos}(f_{osc}t-B)e^{-Ct}+D$ (7)

where I is the measured signal, $A-D$ are fitting constants, and $f_{osc}$ is
the frequency of the damped oscillations, and is directly related to the sound
speed. The thermal conductivity can be calculated by the exponential decay of
the signal after the acoustic portion. A moving average of 100 points was
applied to both the fast and slow signals to decrease fitting times and
further constrain the nonlinear fit to the measured signal.

Presented in Fig. 8, right, are averaged experimental in phase and out of
phase thermal conductivity results for RDX $<$100$>$, $<$010$>$, and $<$001$>$
oriented crystals. The reported thermal conductivities in Fig. 8, right, are
for the indexed crystal orientation. The probed thermal conductivity was
orthogonal to the reported orientation i.e. $<$100$>$ reported will result in
probing of the $<$010$>$ and $<$001$>$ orientations. Probed orientations were
inferred by comparing sound speed measurements performed in these experiments
to velocity measurements from ISTS Sun  _et al._ (2008) and Brillouin
scattering. Bolme and Ramos (2014) RDX has known anisotropy in sound speed,
Sun  _et al._ (2008) and was expected to have anisotropy in thermal
conductivity as well. The thermal conductivity oscillated between 0.28-0.31
W.m-1.K-1 for the $<$100$>$ RDX probed direction, 0.29-0.31 W.m-1.K-1 for
$<$010$>$ RDX probed direction, and 0.275-0.30 W.m-1.K-1 for the $<$001$>$ RDX
probed direction. The corresponding RDX orientations to probe direction are
labelled in Fig. 8, right. The average thermal conductivities are
approximately 0.308, 0.298, and 0.283 W.m-1.K-1 along the $<$001$>$,
$<$100$>$, and $<$010$>$ orientations, respectively.

The aperiodic noise present in the thermal conductivity measurements is mostly
attributed to signal to detector alignment. Minor defects on the surface of
crystals resulted in small pointing changes for the diffracted signal and LO
into the photodiode. RDX crystal habits when formed from acetone lead to
difficulty in cutting certain crystal facets, like the $<$001$>$ RDX
orientation. Connick and May (1969a); Halfpenny, Roberts, and Sherwood (1984)
The small pointing variations in the LO and signal resulted in smaller
measured signals decreasing the signal to noise ratio. Additionally, the
pointing variance increased deviation from the overall sinusoidal trends for
the thermal conductivity as seen in Fig. 8, left for the RDX $<$100$>$ probe
direction. LO pointing was optimized on a shot to shot basis to increase the
signal to noise ratio.

##  IV Discussion

Figure 9: Comparison of the results from this work ($R$NEMD and ISTS) with $R$NEMD results from Izvekov et al., Izvekov, Chung, and Rice (2011) and Sakano et al. Sakano  _et al._ (2018), and experimental results from McGuire & Tarver, McGuire and Tarver (1981) Fedoroff & Clift, Fedoroff and Clift (1960) Rogers, Rogers (1975) Loftus & Gross, Loftus and Gross (1959) Dobratz & Crawford, Dobratz and Crawford (1985) Zinn & Malder, Zinn and Mader (1960) Miller, Miller (1997) Faubion, Faubion (1976) Shoemaker et al., Shoemaker, Stark, and Taylor (1985), Hanson-Parr & Parr, Hanson-Parr and Parr (1999) and Lawless et al.. Lawless  _et al._ (2020) For the latter, the thermal conductivity is actually extracted from finite element simulations to reproduce experimental temperature profiles. The shaded region is bounded by extreme cases for $R$NEMD, determined from parameters and errors in Table 4. Table 5: Thermal conductivity of $\alpha$-RDX at near-ambient temperature and pressure conditions. Source |  |  $\kappa$ (W.m-1.K-1) | T(K) |  $\rho$ (g.cm3)
---|---|---|---|---
This work $R$NEMD |  $<$100$>$ | 0.316–0.344 | 293 | 1.80
|  $<$010$>$ | 0.257–0.306 | 293 | 1.80
|  $<$001$>$ | 0.341–0.394 | 293 | 1.80
This work ISTS |  $<$100$>$ | 0.294–0.301 | 293 | 1.8
|  $<$010$>$ | 0.280-0.287 | 293 | 1.8
|  $<$001$>$ | 0.308-0.309 | 293 | 1.8
Izvekov et al. Izvekov, Chung, and Rice (2011) ($R$NEMD) |  $<$100$>$ | 0.342 | 300 | 1.805
|  $<$010$>$ | 0.317 | 300 | 1.805
|  $<$001$>$ | 0.406 | 300 | 1.805
Sakano et al. Sakano  _et al._ (2018) ($R$NEMD) |  $<$100$>$ | 0.594 | 300 | 1.805
McGuire & Tarver McGuire and Tarver (1981) | - | 0.260 | 293 | 1.8
Fedoroff & Clift Fedoroff and Clift (1960) | - | 0.292 | 29.3 | 1.533
Rogers Rogers (1975) | - | 0.106 | 293 | 1.806
Loftus & Gross Loftus and Gross (1959) | - | 0.205 | 293 | 1.650
Dobratz & Crawford Dobratz and Crawford (1985) | - | 0.106 | 295 | 1.810
Zinn & Malder Zinn and Mader (1960) | - | 0.293 | 293 | 1.8
Miller Miller (1997) | - | 0.193 | 293 | 1.716
Faubion Faubion (1976) | - | 0.073 | 293 | 1.660
Hanson-Parr & Parr Hanson-Parr and Parr (1999) | - | 0.213–0.305333Authors applied an empirical correction to $\kappa$ in order to account for sample porosity. | 293 | 1.60–1.64
Lawless et al. Lawless  _et al._ (2020) | - | 0.09–0.22 | 300 | 1.105–1.703444$\kappa$ is obtained via finite element analysis simulations to match experimental temperature profiles.

The results above provide $\kappa(T,P)$ along the three primitive lattice
vectors of $\alpha$-RDX, in a relevant regime of temperatures and pressures.
There are other crystalline orientations that would be relevant in the
crystal; for instance, $(210)$ and $(111)$ surfaces are observed in
experimental samples, Connick and May (1969b) and thermal conductivity through
the direction perpendicular to those would be useful. However, $\alpha$-RDX
has an orthorhombic unit cell (space group $Pbca$), where a two-rank tensor
such as the thermal conductivity has only three non-zero components, all
independent. Kelly and Groves (1970) As a result, the results presented in
this work allow for the full determination of the thermal conductivity, in any
arbitrary direction.

We compare the thermal conductivity obtained from this work with results from
the literature in Fig. 9 and Table 5. We use the parameters and associated
errors provided in Table 4 to obtain limiting values for the thermal
conductivity, which is bounded at $P$=0 by
($\kappa_{0}-\delta(\kappa_{0}))+\frac{\alpha-\delta(\alpha)}{T}$ and
$(\kappa_{0}+\delta(\kappa_{0}))+\frac{\alpha+\delta(\alpha)}{T}$, see Fig. 9.
Firstly, we note the good agreement between our $R$NEMD results and the ISTS
points obtained in this work (results from ISTS can resolve the anisotropy in
thermal conductivity as the probe is diffracted only in the thermal grating
direction): the trend for orientation dependence is the same for the two
methods, and the values agree between 1 and 20%, depending on the orientation,
as it appears that the $R$NEMD overestimates $\kappa_{100}$ and
$\kappa_{001}$. Taking into account the limiting cases shown in Fig. 9, we
find that $\kappa_{010}$ from ISTS is within the predicted range, while the
values for $\kappa_{100}$ and $\kappa_{001}$ are now within 6 and 11% of the
lower predictions, respectively. Considering possible experimental
uncertainty, the accuracy of the SB-FF, and additional errors arising from the
simulations and fit of the data, this level of agreement is very encouraging.
It also suggests that the results obtained with $R$NEMD can indeed be used to
model the thermal conductivity of HE crystals for which experimental
measurements are not available and non-trivial to obtain (pressure dependence,
in particular, would require the precise alignment of crystals within diamond
anvil cells and a simultaneous control of the temperature). Simulations can
also guide experiments to focus on specific results, like the cross-over
behavior observed in Fig. 6 between $\kappa_{100}$ and $\kappa_{010}$ as a
function of pressure.

As mentioned earlier, our $R$NEMD results also compare well with the work of
Izvekov et al. using a similar method, Izvekov, Chung, and Rice (2011) with
$\kappa_{001}>\kappa_{100}>\kappa_{010}$ in the 225–300 K temperature range at
zero pressure, and the average thermal conductivities
$\bar{\kappa}=(\kappa_{100}+\kappa_{010}+\kappa_{001})/3$ at 300 K within 3%
of each other. Our results do not reproduce the maximum in the thermal
conductivity observed by Izvekov et al. at 275 K, however that is most likely
a fortuitous result, as discussed in Sec. III.2 and acknowledged in Ref.
Izvekov, Chung, and Rice, 2011. An additional result obtained via $R$NEMD by
Sakano et al. Sakano  _et al._ (2018) for $<$001$>$ oriented RDX is also
shown, higher than the ISTS results by over 90%, and both our ($\sim$60%) and
Izvekov’s ($\sim$45%) $R$NEMD results. Sakano et al. used a similar setup as
ours, but different sample sizes and parameters, notably a swapping frequency
of 3 ps vs. 250 fs in our case. While Sakano et al. state that smaller swap
frequencies induce “nonconstant heat flux and too large thermal gradient”, we
have not observed such issues. Similarly, while Sakano et al. fitted the
$1/\kappa$ vs. $1/L$ relation for $L$ between 32 and 59 nm due to a breakdown
of the linearity for smaller samples, ours and Izvekov’s results show the
appropriate linear behavior for samples between 8 and 25 nm (as seen in Fig.
5), and 6 to 55 nm, respectively. These results, in addition to previous
considerations on sufficient sample size in a similar molecular crystal
Chitsazi  _et al._ (2020), and the good agreement with ISTS results, give us
confidence that the present results are representative of $\alpha$-RDX.

Finally, Fig. 9 shows the experimental results of Miller Miller (1997) and
Shoemaker et al., Shoemaker, Stark, and Taylor (1985) which were fitted to a
non-monotonic relation by Miller, although at this scale the error bars make
it difficult to assess the true form of the curve. Additional experimental
results are shown, most of it for polycrystalline or powder samples with
density lower than that of single crystal $\alpha$-RDX under similar
conditions, with results averaged over many crystal orientations which lowered
the thermal conductivity values. As such, direct comparison with our
predictions for single crystals is difficult. However, our results agree well
with the experimental results presented by Hanson-Parr & Parr Hanson-Parr and
Parr (1999) (see Fig. 9), who also used powders but explicitly corrected for
porosity. Thermal conductivity from different sources at or near ambient
conditions are also summarized in Table 5. Additional ISTS experiments are
currently underway on oriented single crystals to probe more temperatures and
test the $R$NEMD trends.

##  V Conclusion

We have used $R$NEMD simulations with the Smith-Bharadwaj non-reactive force
field to determine the thermal conductivity, $\kappa$, of $\alpha$-RDX in the
$<$100$>$, $<$010$>$, and $<$001$>$ directions, as a function of temperature
and pressure for $T$=200-400 K and $P$=0–2 GPa. The SB-FF is shown to
reproduce the lattice constants and thermal expansion coefficients of
$\alpha$-RDX in good agreement with experiments. We find that $\kappa$ varies
linearly with the inverse temperature, in agreement with theory. In addition,
we also observe a linear dependence of $\kappa$ with pressure. As a result, we
suggest a simple linear form for $\kappa(T,P)$ for each orientation, and
extract the corresponding coefficients which are then used to extrapolate
$\kappa(T,P)$ up to 500 K and 5 GPa.

Critically, ISTS measurements on RDX single crystals were performed to provide
the first direct validation between MD and experimental results for
orientation-dependent thermal conductivity in HE crystals. ISTS results
validate the $R$NEMD anisotropy trend, and the numbers agree between 1 and 20%
depending on the orientation, for an average agreement around 10%. These
results can be used to parameterized mesoscale models that require
anisotropic, temperature and pressure dependent, thermal properties for
$\alpha$-RDX.

##  VI Acknowledgements

Work presented in this article was supported by the Laboratory Directed
Research and Development program of Los Alamos National Laboratory under
project number 20180100DR. This research used resources provided by the Los
Alamos National Laboratory Institutional Computing Program. Los Alamos
National Laboratory is operated by Triad National Security, LLC, for the
National Nuclear Security Administration of U.S. Department of Energy
(Contract No. 89233218CNA000001).

## References

  * Bowden and Yoffe (2009) F. P. Bowden and Y. D. Yoffe,  _Initiation and Growth of Explosion in Liquids and Solids_ (Cambridge University Press, 2009).
  * Storm, Stine, and Kramer (1990) C. B. Storm, J. R. Stine, and J. F. Kramer, “Sensitivity relationships in energetic materials,” in  _Chemistry and Physics of Energetic Materials_ , NATO Advanced Science Institutes Series, Series C, Mathematical and Physical Sciences, Vol. 309, edited by S. N. Bulusu (Kluwer Academic Publ., 1990) pp. 605–639.
  * McGuire and Tarver (1981) R. R. McGuire and C. M. Tarver, “Chemical decomposition models for the explosion of confined HMX, TATB, RDX, and TNT explosives,” in  _Proceedings of the 7th International Symposium on Detonation_ (Office of Naval Research, Arlington, VA, 1981) pp. 56–64.
  * Fedoroff and Clift (1960) B. T. Fedoroff and G. D. Clift,  _Encyclopedia of Explosives and Related Items_ (Picatinny Arsenal, Dover, NJ, 1960).
  * Rogers (1975) R. N. Rogers, “Thermochemistry of explosives,” Thermochim. Acta 11, 131–139 (1975).
  * Loftus and Gross (1959) G. Loftus and G. Gross,  _Thermal and Self Ignition Properties of Six Explosives, NBS Report 6548_ (National Bureau of Standards, Gaithersburg, MD, 1959).
  * Dobratz and Crawford (1985) B. M. Dobratz and P. C. Crawford,  _LLNL Explosives Handbook, Properties of Chemical Explosives and Explosive Simulants, Report UCRL 52997, Rev. 2_ (LLNL, Livermore, CA, 1985).
  * Zinn and Mader (1960) J. Zinn and C. L. Mader, “Thermal initiation of explosives,” J. Appl. Phys. 31, 323–328 (1960).
  * Miller (1997) M. S. Miller,  _Thermophysical Properties of RDX, ARL-TR-1319_ (Army Research Lab, Aberdeen Proving Ground, MD, 1997).
  * Faubion (1976) B. D. Faubion,  _Thermal Conductivity of RDX_ (Mason and Hanger-Silas Mason Co., Inc., Amarillo, TX, 1976).
  * Shoemaker, Stark, and Taylor (1985) R. L. Shoemaker, J. A. Stark, and R. E. Taylor, “Thermophysical properties of propellants,” HTHP 17, 429–435 (1985).
  * Hanson-Parr and Parr (1999) D. M. Hanson-Parr and T. P. Parr, “Thermal properties measurements of solid rocket propellant oxidizers and binder materials as a function of temperature,” J. Energ. Mater. 17, 1–48 (1999).
  * Lawless  _et al._ (2020) Z. D. Lawless, M. L. Hobbs, M. J. Kaneshige, Z. D. Lawless, M. L. Hobbs, and M. J. Kaneshige, “Thermal conductivity of energetic materials,” J. Energ. Mater. 38, 214–239 (2020).
  * Kittel (1966) C. Kittel, “Thermal poperties of insulators,” in  _Introduction to Solid State Physics, Third Edition_ (John Wiley & Sons, 1966) Chap. 6.
  * Kroonblawd and Sewell (2016a) M. P. Kroonblawd and T. D. Sewell, “Theoretical determination of anisotropic thermal conductivity for crystalline 1,3,5-triamino-2,4,6-trinitrobenzene (TATB),” J. Chem. Phys. 139, 074503 (2016a).
  * Kroonblawd and Sewell (2016b) M. P. Kroonblawd and T. D. Sewell, “Theoretical determination of anisotropic thermal conductivity for initially defect-free and defective TATB single crystals,” J. Chem. Phys. 141, 184501 (2016b).
  * Bedrov, Smith, and Sewell (2000) D. Bedrov, G. D. Smith, and T. D. Sewell, “Thermal conductivity of liquid octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX) from molecular dynamics simulations,” Chem. Phys. Lett. , 64–68 (2000).
  * Chitsazi  _et al._ (2020) R. Chitsazi, M. P. Kroonblawd, A. Pereverzev, and T. Sewell, “A molecular dynamics simulation study of thermal conductivity anisotropy in $\beta$-octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine ($\beta$-HMX) ,” Model. Simul. Mater. Sci. Eng. 28, 025008 (2020).
  * Müller-Plathe (1997) F. Müller-Plathe, “A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity,” J. Chem. Phys. 106, 6082 (1997).
  * Zhang  _et al._ (2005) M. Zhang, E. Lussetti, L. E. S. de Souza, and F. Müller-Plathe, “Thermal conductivities of molecular liquids by reverse nonequilibrium molecular dynamics,” J. Phys. Chem. B 109, 15060–15067 (2005).
  * Munday  _et al._ (2011) L. B. Munday, P. W. Chung, B. M. Rice, and S. D. Solares, “Simulations of high-pressure phases in RDX,” J. Phys. Chem. B 115, 4378–4386 (2011).
  * Josyula, Rahul, and De (2014) K. Josyula, Rahul, and S. De, “Thermomechanical properties and equation of state for the gamma-polymorph of hexahydro-1,3,5-trinitro-1,3,5-triazine,” RSC Adv. 4, 41491–41499 (2014).
  * Hooks  _et al._ (2015) D. E. Hooks, K. J. Ramos, C. A. Bolme, and M. J. Cawkwell, “Elasticity of crystalline molecular explosives,” Propellants, Explos. Pyrotech. 40, 333–350 (2015).
  * Olinger, Roof, and Cady (1978) B. Olinger, B. Roof, and H. Cady, “The linear and volume compression of $\beta$-HMX and RDX to 9 GPa,” in  _Proc. Symposium (Intern.) on High Dynamic Pressures (C.E.A., Paris, France)_ (1978) pp. 3–8.
  * Davidson  _et al._ (2008) 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, “Explosives under pressure–the crystal structure of $\gamma$-RDX as determined by high-pressure X-ray and neutron diffraction,” Cryst. Eng. Comm. 10, 162–165 (2008).
  * Dreger (2012) Z. A. Dreger, “Energetic materials under high pressures and temperatures: stability, polymorphism and decomposition of RDX,” J. Phys. Conf. Ser. 377, 012047 (2012).
  * Cawkwell  _et al._ (2016) M. J. Cawkwell, D. J. Luscher, F. L. Addessio, and K. J. Ramos, “Equations of state for the $\alpha$ and $\gamma$ polymorphs of cyclotrimethylene trinitramine,” J. Appl. Phys. 119 (2016).
  * Smith and Bharadwaj (1999) G. D. Smith and R. K. Bharadwaj, “Quantum chemistry based force field for simulations of HMX,” J. Phys. Chem. B 103, 3570–3575 (1999).
  * Bedrov  _et al._ (2001) D. Bedrov, C. Ayyagari, G. D. Smith, T. D. Sewell, R. Menikoff, and J. M. Zaug, “Molecular dynamics simulations of HMX crystal polymorphs using a flexible molecule force field,” J. Comput. Mater. Des. 8, 77–85 (2001).
  * Weingarten and Sausa (2015) N. S. Weingarten and R. C. Sausa, “Nanomechanics of RDX single crystals by force displacement measurements and molecular dynamics simulations,” J. Phys. Chem. A 119, 9338–9351 (2015).
  * Josyula, Rahul, and De (2019) K. Josyula, Rahul, and S. De, “In silico study of $\alpha$-$\gamma$ phase transformation in hexahydro-1,3,5-trinitro-1,3,5-triazine,” Comput. Mater. Sci. 170, 109180 (2019).
  * Cawkwell  _et al._ (2010) M. J. Cawkwell, K. J. Ramos, D. E. Hooks, and T. D. Sewell, “Homogeneous dislocation nucleation in cyclotrimethylene trinitramine under shock loading,” J. Appl. Phys. 107 (2010).
  * Ramos  _et al._ (2010) K. J. Ramos, D. E. Hooks, T. D. Sewell, and M. J. Cawkwell, “Anomalous hardening under shock compression in (021)-oriented cyclotrimethylene trinitramine single crystals,” J. Appl. Phys. 108 (2010).
  * Bidault and Pineau (2018) X. Bidault and N. Pineau, “Granularity impact on hotspot formation and local chemistry in shocked nanostructured RDX,” J. Chem. Phys. 149, 224703 (2018).
  * Schelling, Phillpot, and Keblinski (2002) P. K. Schelling, S. R. Phillpot, and P. Keblinski, “Comparison of atomic-level simulation methods for computing thermal conductivity,” Phys. Rev. B 65, 144306 (2002).
  * Izvekov, Chung, and Rice (2011) S. Izvekov, P. W. Chung, and B. M. Rice, “Non-equilibrium molecular dynamics simulation study of heat transport in hexahydro-1,3,5-trinitro-s-triazine (RDX),” Int. J. Heat Mass Transf. 54, 5623–5632 (2011).
  * Plimpton (1995) S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” J. Comput. Phys. 117, 1–19 (1995).
  * Mathew  _et al._ (2018) N. Mathew, M. P. Kroonblawd, T. Sewell, and D. L. Thompson, “Predicted melt curve and liquid-state transport properties of TATB from molecular dynamics simulations,” Mol. Simulat. 44, 613–622 (2018).
  * McCrone (1950) W. C. McCrone, “Crystallographic data. 32. RDX (cyclotrimethylenetrinitramine),” Anal. Chem. 22, 954–955 (1950).
  * Choi and Prince (1972) C. S. Choi and E. Prince, “The crystal structure of cyclotrimethylenetrinitramine,” Acta Crystallogr. B 28, 2857–2862 (1972).
  * Kroonblawd  _et al._ (2016) M. P. Kroonblawd, N. Mathew, S. Jiang, and T. D. Sewell, “A generalized crystal-cutting method for modeling arbitrarily oriented crystals in 3D periodic simulation cells with applications to crystal–crystal interfaces,” Comput. Phys. Commun. 207, 232–242 (2016).
  * Abramson, Brown, and Slutsky (1999) E. H. Abramson, J. M. Brown, and L. J. Slutsky, “Applications of impulsive stimulated scattering in the Earth and planetary sciences,” Annu. Rev. Phys. Chem. 50, 279–313 (1999).
  * Dennett and Short (2018) C. A. Dennett and M. P. Short, “Thermal diffusivity determination using heterodyne phase insensitive transient grating spectroscopy,” J. Appl. Phys. 123, 215109 (2018).
  * Johnson  _et al._ (2012) J. A. Johnson, A. A. Maznev, M. T. Bulsara, E. A. Fitzgerald, T. C. Harman, S. Calawa, C. J. Vineis, G. Turner, and K. A. Nelson, “Phase-controlled, heterodyne laser-induced transient grating measurements of thermal transport properties in opaque material,” J. Appl. Phys. 111, 023503 (2012).
  * Käding  _et al._ (1995) O. W. Käding, H. Skurk, A. A. Maznev, and E. Matthias, “Transient thermal gratings at surfaces for thermal characterization of bulk materials and thin films,” Appl. Phys. A 261, 253–261 (1995).
  * Maznev, Nelson, and Rogers (1998) A. A. Maznev, K. A. Nelson, and J. A. Rogers, “Optical heterodyne detection of laser-induced gratings,” Opt. Lett. 23, 1319–1321 (1998).
  * Rogers, Maznev, and Nelson (2012) J. A. Rogers, A. Maznev, and K. A. Nelson, “Impulsive stimulated thermal scattering,” in  _Characterization of Materials_ (American Cancer Society, 2012) pp. 1–17.
  * Rogers  _et al._ (2000) J. A. Rogers, A. A. Maznev, M. J. Banet, and K. A. Nelson, “Optical generation and characterization of acoustic waves in thin films: Fundamentals and applications,” Annu. Rev. Mater. Sci. 30, 117–157 (2000).
  * Rogers, Yang, and Nelson (1994) J. Rogers, Y. Yang, and K. Nelson, “Elastic modulus and in-plane thermal diffusivity measurements in thin polyimide films using symmetry-selective real-time impulsive stimulated thermal scattering,” Appl. Phys. 58, 523–534 (1994).
  * Tokmakoff, Banholzer, and Fayer (1993) A. Tokmakoff, W. F. Banholzer, and M. D. Fayer, “Thermal diffusivity measurements of natural and isotopically enriched diamond by picosecond infrared transient grating experiments,” Appl. Phys. A 56, 87–90 (1993).
  * Lazarz  _et al._ (2020) J. D. Lazarz, S. D. McGrane, R. Perriot, C. Bolme, M. J. Cawkwell, and K. J. Ramos, “Anisotropic thermal conductivity and elasticity of RDX using impulsive stimulated thermal scattering,” AIP Conf. Proc. 2272, 060023 (2020).
  * Sorescu and Rice (2010) D. C. Sorescu and B. M. Rice, “Theoretical predictions of energetic molecular crystals at ambient and hydrostatic compression conditions using dispersion corrections to conventional density functionals,” J. Phys. Chem. C 114, 6734–6748 (2010).
  * Hunter  _et al._ (2013) S. Hunter, T. Sutinen, S. F. Parker, C. A. Morrison, D. M. Williamson, S. Thompson, P. J. Gould, and C. R. Pulham, “Experimental and DFT-D studies of the molecular organic energetic material RDX,” J. Phys. Chem. C 117, 8062–8071 (2013).
  * Bolotina and Pinkerton (2015) N. B. Bolotina and A. A. Pinkerton, “Temperature dependence of thermal expansion tensors of energetic materials,” J. Appl. Crystallogr. 48, 1364–1380 (2015).
  * Sun  _et al._ (2011) J. Sun, X. Shu, Y. Liu, H. Zhang, X. Liu, Y. Jiang, B. Kang, C. Xue, and G. Song, “Investigation on the thermal expansion and theoretical density of 1,3,5-trinitro-1,3,5-triazacyclohexane,” Propellants, Explos. Pyrotech. 36, 341–346 (2011).
  * Cady (1972) H. H. Cady, “Coefficient of thermal expansion of pentaerythritol tetranitrate and hexahydro-1,3,5-trinitro-1,3,5-triazine (RDX),” J. Chem. Eng. Data 17, 369–371 (1972).
  * Bevington (1969) P. R. Bevington, in  _Data Reduction and Error Analysis for the Physical Sciences_ (McGraw-Hill Book Company, 1969) Chap. 4.
  * R Core Team (2017) R Core Team,  _R: A Language and Environment for Statistical Computing_, R Foundation for Statistical Computing, Vienna, Austria (2017).
  * Klemens (1969) P. G. Klemens, “Thermal conductivity of solids,” in  _Thermal Conductivity_ , edited by R. P. Tye (Academic Press, 1969) Chap. 1.
  * Rey-Lafon and Bonjour (1973) M. Rey-Lafon and E. Bonjour, “Étude de la chaleur specifique de la trinitro-1,3,5 hexahydro-s-triazine cristallisée détermination experimentale et calcul à partir des fréquences de vibration infrarouges et raman,” Mol. Cryst. Liq. Cryst. 24, 191–199 (1973).
  * Miller (1995) M. S. Miller, “Thermophysical properties of cyclotrimethylenetrinitramine,” J. Thermophys. Heat Transfer 8, 803–805 (1995).
  * Long  _et al._ (2012) Y. Long, Y. G. Liu, F.-D. Nie, and J. Chen, “A method to calculate the thermal conductivity of HMX under high pressure,” Philos. Mag. 92, 1023–1045 (2012).
  * Fan  _et al._ (2017) H. Fan, Y. Long, L. Ding, J. Chen, and F.-D. Nie, “A theoretical study of elastic anisotropy and thermal conductivity for TATB under pressure,” Comput. Mater. Sci. 131, 321–332 (2017).
  * Nakanishi, Nagasawa, and Murakami (1982) N. Nakanishi, A. Nagasawa, and Murakami, “Lattice stability and soft modes,” J. Phys. Colloq. 43 (C4), C4–35–C4–55 (1982).
  * Leiding  _et al._ (2021) J. A. Leiding, K. A. Velizhanin, R. Perriot, M. J. Cawkwell, T. D. Aslam, and S. A. Andrews, “A new parameterization of a 1-step thermal decomposition model of PBX-9501,” to be sumitted (2021).
  * Sun  _et al._ (2008) B. Sun, J. M. Winey, N. Hemmi, Z. A. Dreger, K. A. Zimmerman, Y. M. Gupta, D. H. Torchinsky, and K. A. Nelson, “Second-order elastic constants of pentaerythritol tetranitrate and cyclotrimethylene trinitramine using impulsive stimulated thermal scattering,” J. Appl. Phys. 104, 073517 (2008).
  * Bolme and Ramos (2014) C. A. Bolme and K. J. Ramos, “The elastic tensor of single crystal rdx determined by brillouin spectroscopy,” J. Appl. Phys. 116, 183503 (2014).
  * Connick and May (1969a) W. Connick and F. May, “Dislocation etching of cyclotrimethylene trinitramine crystals,” J. Cryst. Growth 5, 65–69 (1969a).
  * Halfpenny, Roberts, and Sherwood (1984) P. Halfpenny, K. Roberts, and J. Sherwood, “Dislocations in energetic materials: IV. The crystal growth and perfection of cyclotrimethylene trinitramine (RDX),” J. Cryst. Growth 69, 73–81 (1984).
  * Sakano  _et al._ (2018) M. Sakano, B. Hamilton, M. M. Islam, and A. Strachan, “Role of molecular disorder on the reactivity of RDX,” J. Phys. Chem. C 122, 27032–27043 (2018).
  * Connick and May (1969b) W. Connick and F. G. J. May, “Dislocation etching of cyclotrimethylene trinitramine crystals,” J. Cryst. Growth 5, 65–69 (1969b).
  * Kelly and Groves (1970) A. Kelly and G. W. Groves, in  _Crystallography and Crystal Defects_ (Longman, 1970) Chap. 4.