## Appendix B Effective Energy Drift Calculation  
The effective energy drift [24; 30] for the time-reversible integrator developed in Section 2.3.2 could be computed by introducing an auxiliary momentum-like variable \(\boldsymbol{\alpha}\) independent of the state of the system and with a Gaussian equilibrium distribution \(\mathcal{P}(\boldsymbol{\alpha})\propto e^{-\boldsymbol{\alpha}^{2}/2}\), and then using the high-friction limit of the integrator introduced in Ref. [30]. This allowed to write the shape update \(\mathbf{h}_{i+1}=\widetilde{\boldsymbol{\mu}}\,\mathbf{h}_{i}\), with \(\widetilde{\boldsymbol{\mu}}\) defined in Equation (16), as:  
\[\boldsymbol{\alpha}_{i} =\widetilde{\mathcal{R}}_{i}, \tag{12a}\] \[\boldsymbol{\alpha}_{i+1/2} =\boldsymbol{\alpha}_{i}+\widetilde{\mathbf{A}}_{i}\frac{\sqrt{ \Delta t}}{b}\,,\] (12b) \[\mathbf{h}_{i+1} =\exp\left(b\sqrt{\Delta t}\,\boldsymbol{\alpha}_{i+1/2}\right) \mathbf{h}_{i}\,,\] (12c) \[\boldsymbol{\alpha}_{i+1} =\widetilde{\mathcal{R}}_{i+1}\,, \tag{12d}\]where \(b=\sqrt{\frac{2\beta\gamma k_{B}T}{3\gamma^{\prime}\tau_{p}}}\) and \(i+1/2\) is an intermediate time between \(i\) and \(i+1\). \(V^{\prime}\) is the volume obtained from Equation (15) after the first half-step isotropic propagation. Note that \(\widetilde{\mathbf{A}}_{i}\) instead was computed using the non-propagated volume \(V\). The contribution to the effective energy drift given by Equation (A3) was defined as:  
\[\Delta\widetilde{H}_{i\to i+1}=-k_{B}T\log\frac{\mathcal{P}\big{(}- \boldsymbol{\kappa}^{i+1},\mathbf{h}^{i+1}\big{)}\,\Pi\big{(}(-\boldsymbol{ \kappa}^{i+1},\mathbf{h}^{i+1})\to(-\boldsymbol{\kappa}^{i},\mathbf{h}^{i}) \big{)}}{\mathcal{P}\big{(}\boldsymbol{\kappa}^{i},\mathbf{h}^{i}\big{)}\,\Pi \big{(}(\boldsymbol{\kappa}^{i},\mathbf{h}^{i})\to(\boldsymbol{\kappa}^{i+1}, \mathbf{h}^{i+1})\big{)}}\,, \tag{14}\]  
where the factors \(\Pi\) are the transition probabilities of the forward and backward moves. Note that the backward move was obtained by changing the sign of the momentum-like variables \(\boldsymbol{\kappa}\), in a generalized detailed-balance fashion. The contributions to Equation (14) could be computed by writing the joint probability distributions of the elements of the matrix \(\widetilde{\boldsymbol{\mathcal{R}}}\) as functions of \(\boldsymbol{\kappa}\) and \(\mathbf{h}\), keeping in account that \(b\) was untouched in the change of the shape. It was also convenient to define \(\Delta\boldsymbol{\epsilon}=\widetilde{\mathbf{A}}_{i}\,\Delta t+b\sqrt{ \Delta t}\,\widetilde{\boldsymbol{\mathcal{R}}}_{i}\) as the argument of the matrix exponential of Equation (16). As a final result, the effective energy drift produced by the change of shape was computed as:  
\[\Delta\widetilde{H}=\Delta K+\Delta U+\Delta U_{\text{el}}+k_{B}T\sum_{ \alpha\beta}\!\left[\frac{\Delta t}{2b^{2}}\Delta(\widetilde{A}_{\alpha\beta })^{2}+\frac{1}{b^{2}}\Delta\epsilon_{\alpha\beta}\Big{(}\widetilde{A}_{ \alpha\beta}(t)+\widetilde{A}_{\alpha\beta}(t+\Delta t)\Big{)}\right]. \tag{15}\]  
By substituting \(b\) in Equation (15), one obtains:  
\[\Delta\widetilde{H}=\Delta K+\Delta U+\Delta U_{\text{el}}+\frac{3V^{\prime} \tau_{p}}{2\beta\tau_{p}}\sum_{\alpha\beta}\!\left[\Delta(\widetilde{A}_{ \alpha\beta})^{2}\frac{\Delta t}{2}+\Delta\epsilon_{\alpha\beta}\Big{(} \widetilde{A}_{\alpha\beta}(t)+\widetilde{A}_{\alpha\beta}(t+\Delta t)\Big{)} \right]. \tag{16}\]  
Note that the term proportional to \(\Delta\epsilon_{\alpha\beta}\) tends to cancel the increment in the kinetic and potential energies in the limit of small \(\Delta t\). The term proportional to \(\Delta(\widetilde{A}_{\alpha\beta})^{2}\) has a prefactor \(V^{\prime}\) that changes at each step, and, thus, leads to an accumulation along the simulated trajectory. This is different from the corresponding term in Equation (14) of Ref. [30], which only depended on the initial and final state.