*2.1. Theoretical Background and System Preparation*

The elastic stiffness tensor *Cij* (written here in Voigt notation) was calculated using the direct method, which is based on the following equation

$$\mathbb{C}\_{ij} = \frac{\partial \sigma\_i}{\partial \varepsilon\_j}. \tag{1}$$

Here, *ε<sup>j</sup>* denotes one of six linear strain components and *σ<sup>i</sup>* is one of six Cauchy stress components. For monoclinic crystals such as *β*-HMX, the stiffness tensor contains 13 independent non-zero elastic coefficients. We chose linear strain rather than the often-used Lagrange strain because it is work-conjugate to the Cauchy stress (see [13], Section 2.2.4) and it is the Cauchy stress that is typically provided as output from MD simulations. For the small strains used in this work we can expect any differences arising due to the use of linear rather than Lagrange strains to be negligible. The partial derivatives in Equation (1) were calculated numerically using finite differences. To achieve this, we adopted the following four-step procedure.

First, a crystal cell of *β*-HMX (P21/n space group setting) was equilibrated for 1 ns at the temperature and pressure of interest using the isobaric-isothermal (NPT) ensemble. Next, 100 ps of isochoric-isothermal (NVT) simulation was performed using the time-averaged crystal cell parameters from the NPT trajectory. In the third step, small linear strains were imposed on the NVT-equilibrated crystal cell and NVT simulations were performed on the strained cells. Finally, the elastic coefficients *Cij* were obtained by numerical differentiation of the stresses with respect to strains.

We used a 3D-periodic crystal cell consisting of 6 × 4 × 5 unit-cell replications along the *a*, *b*, and *c* crystallographic axes, respectively. This crystal cell contains 240 HMX molecules (6720 atoms). Similar cell sizes for both *β*-HMX and other molecular-crystal high explosives are thought to be large enough to produce size-converged elastic properties of bulk crystals [6,14]. Temperatures *T* = 300, 500, 700, 900, and 1100 K and pressures *P* = 10<sup>−</sup>4, 5, 10, 20, and 30 GPa were considered. These pressures were chosen to match the interval considered in Ref. [15]. At atmospheric pressure (10−<sup>4</sup> GPa), *β*-HMX is known to undergo a phase transition to the *α* polymorph (orthorhombic) at approximately 420 K and the *δ* polymorph (hexagonal) when heated above approximately 450 K [16,17]; *δ*-HMX melts at about 550 K. However, with the force field that we use in this work the crystal remains stable (in the sense that melting does not occur) at 700 K and 1 atm for simulations lasting several nanoseconds. The crystal does melt at 1 atm somewhere between 700 and 900 K. Therefore, for the pressure *P* = 1 atm, we only considered temperatures up to 700 K. However, we observed that *β*-HMX remains crystalline all the way to 1100 K at 5 GPa and higher pressures. This is consistent with the MD-based melt curve of *β*-HMX reported recently by Kroonblawd and Austin [18]. They predicted that at 5 GPa *β*-HMX melts at 1320 K and the melting temperatures get higher for pressures above 5 GPa. Thus, overall, we studied 23 temperature and pressure combinations. We emphasize that the results below for the higher temperatures should be interpreted with caution because HMX may decompose on relevant time scales at those temperatures, at least at 1 atm [16].

When choosing the strain size, one has to consider the following arguments. On the one hand, one wants the strain to be sufficiently small to yield a reliable finite-difference approximation for the numerical evaluation of derivatives. On the other hand, if the strain is too small, the resulting change in stress is also small and, because there is a significant thermal noise in stress calculations, long simulation times are required to obtain converged results. Considering these factors, strains *ε<sup>i</sup>* were chosen to be ±0.004 when *i* = 1, 2 or 3 and ±0.008 when *i* = 4, 5 or 6. The six strains required to obtain the full set of elastic coefficients were applied by distorting the simulation box in the following way. If the edge vectors of an unperturbed triclinic cell are given by vectors **a**, **b**, and **c**, then the geometry of the simulation box can, with no loss of generality, be described by the upper triangular matrix

$$h\_0 = \begin{pmatrix} a\_x & b\_x & c\_x \\ 0 & b\_y & c\_y \\ 0 & 0 & c\_z \end{pmatrix} \tag{2}$$

by aligning **a** with the positive *x* axis and requiring **b** to lie in the *xy* plane [19]. Then, the upper triangular matrix *h* for the the strained cell is given by

$$h = (I \pm \epsilon\_i) h\_{0\prime} \tag{3}$$

where *I* is the 3 × 3 identity matrix and  *<sup>i</sup>* is one of the following six upper triangular matrices

$$
\varepsilon\_1 = \begin{pmatrix} \varepsilon\_1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad \varepsilon\_2 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \varepsilon\_2 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad \varepsilon\_3 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \varepsilon\_3 \end{pmatrix}, \tag{4}
$$

$$
\varepsilon\_4 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & \varepsilon\_4 \\ 0 & 0 & 0 \end{pmatrix}, \qquad \epsilon\_5 = \begin{pmatrix} 0 & 0 & \varepsilon\_5 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad \epsilon\_6 = \begin{pmatrix} 0 & \varepsilon\_6 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \tag{5}
$$

Each of the strained simulation cells was equilibrated in the NVT ensemble for 100 ps after which the average stress tensor over a 4 ns NVT production trajectory was accumulated. For the pressure of 1 atm and temperatures of 300, 500, and 700 K the changes in stress tensor were small because the crystal is less rigid at the low pressure. Therefore, for those three cases, the trajectory was run for 12 ns to achieve better convergence of the elastic coefficients. Once the stress tensor components were obtained, the numerical derivatives in Equation (1) were evaluated as follows

$$\mathbf{C}\_{ij} \approx \frac{\sigma\_i(\boldsymbol{\varepsilon}\_j) - \sigma\_i(-\boldsymbol{\varepsilon}\_j)}{2\boldsymbol{\varepsilon}\_j},\tag{6}$$

where *σi*(*εj*) denotes the time-averaged strain tensor component for the equilibrated crystal cell with strain *ε<sup>j</sup>* applied. Because of the elastic tensor symmetry, we expect *Cij* = *Cji* when *i* = *j*. This was indeed the case, with differences between *Cij* and *Cji* typically less than 0.2 GPa. For these off-diagonal elements of the elastic tensor, we report the average of *Cij* and *Cji*. As part of our analysis we also obtained the eight elastic coefficients that are expected to be zero in a monoclinic crystal. These coefficients were, indeed, very small, never exceeding 0.08 GPa in magnitude. Three examples of the full 6 × 6 *Cij* matrices obtained in this work are given in the Supplementary Materials. In addition to the elastic coefficients, the Voigt, Reuss, and Voigt–Reuss–Hill average bulk and shear moduli were calculated using standard expressions [8,20].
