## 2 Materials and Methods
### Integrating the Equations of Motion
#### 2.3.2 Time-Reversible Integrator  
The time-reversible algorithm was implemented as follows:  
1. Propagate volume for \(\Delta t/2\). 2. Propagate cell shape at constant volume. 3. Propagate volume for \(\Delta t/2\). The first step was performed by propagating \(\sqrt{V}\) with the Euler algorithm [18] for a time \(\Delta t/2\):  
\[V^{\prime}=\left(\sqrt{V}-\frac{\beta_{T}\sqrt{V}}{2\tau_{p}}\bigg{(}P_{0}- \frac{\text{Tr}\mathbf{P}_{\text{int}}}{3}+\frac{\text{Tr}(\mathbf{h}\Sigma \mathbf{h}^{T})}{3V}-\frac{k_{B}T}{2V}\bigg{)}\frac{\Delta t}{2}+\sqrt{\frac {k_{B}T\beta_{T}}{2\tau_{p}}}\sqrt{\frac{\Delta t}{2}}\mathcal{R}\right)^{2}, \tag{15}\]  
where \(\mathcal{R}\) is a standard Gaussian number. The second step was performed by applying a deformation matrix computed as:  
\[\widetilde{\mathbf{\mu}}=\exp\left(\widetilde{\mathbf{A}}\Delta t+\sqrt{\frac{2 \beta_{T}k_{B}T\Delta t}{3V^{\prime}\tau_{p}}}\widetilde{\mathbf{\mathcal{R}}} \right), \tag{16}\]where \(\mathbf{\widetilde{A}}=-\frac{\beta_{T}}{ST_{p}}\Big{[}\Big{(}\mathbf{P}_{\text{ int}}-\frac{\text{Tr}\mathbf{P}_{\text{int}}}{3}\Big{)}+\frac{\mathbf{h} \mathbf{\Sigma}\mathbf{h}^{T}}{V}-\frac{\text{Tr}(\mathbf{h}\mathbf{\Sigma} \mathbf{h}^{T})}{3V}\Big{]}\) and \(\widetilde{\mathbf{\mathcal{R}}}=\mathbf{\mathcal{R}}-\frac{\text{Tr} \mathbf{\mathcal{R}}}{3}\). The matrix exponential was computed using the Pade approximation [29]. Since the argument of the exponential was constructed as a traceless matrix, \(\widetilde{\boldsymbol{\mu}}\) has exactly unit determinant, even when computed within this approximation. Note that Equation (16) uses the internal pressure and volumes evaluated _before_ the volume update of Equation (15), _except_ for the prefactor of the noise term. Choosing a noise prefactor proportional to \(\frac{1}{V^{\prime}}\) made sure that in the limit of small \(\Delta t\) or large \(\tau_{P}\), that is when the stochastic term dominates Equation (16), the probability to generate a given deformation matrix \(\widetilde{\boldsymbol{\mu}}\) or its inverse would be equal; thus, minimizing detailed balance violations. Finally, the square root of the volume was evolved again for a time \(\Delta t/2\):  
\[V^{\prime\prime}=\left(\sqrt{V^{\prime}}-\frac{\beta_{T}\sqrt{V}}{2\tau_{P}} \bigg{(}P_{0}-\frac{\text{Tr}\mathbf{P}_{\text{int}}}{3}+\frac{\text{Tr}( \mathbf{h}\mathbf{\Sigma}\mathbf{h}^{T})}{3V}-\frac{k_{B}T}{2V}\bigg{)}\frac{ \Delta t}{2}+\sqrt{\frac{k_{B}T\beta_{T}}{2\tau_{P}}}\sqrt{\frac{\Delta t}{2}} \mathcal{R}^{\prime}\right)^{2}, \tag{17}\]  
where \(\mathcal{R}^{\prime}\) is a standard Gaussian number independently drawn from \(\mathcal{R}\). Note that the combination of Equations (15) and (17) is identical to applying one of the two equations using a twice-as-large \(\Delta t\). Since the volume is unaffected by Equation (16), its dynamics is identical to that resulting from the reversible integrator proposed in Ref. [18]. Overall, the lattice vectors were, thus, multiplied by a scaling matrix \(\boldsymbol{\mu}=\widetilde{\boldsymbol{\mu}}\sqrt[3]{V^{\prime\prime}/V}\). Both the Euler and the time-reversible algorithms were formulated here and implemented in the codes by directly propagating the deformation matrix instead of the cell matrix, since the two formulations were equivalent but the first one minimized the number of matrix operations. Importantly, the equations of motion (Equation (10)) only depend on the instantaneous cell matrix and not on the deformation with respect to the initial cell matrix. Hence, different choices of the initial cell matrix are expected to lead to the same distributions once the system is equilibrated. Following Ref. [18], positions and momenta were evolved, concurrently, with lattice vectors and updated at half step.