## 2 Materials and Methods
### Scaled Coordinates and Internal Pressure  
Whenever the lattice vectors were modified, it was convenient to act on positions and, optionally, momenta, so as to keep them consistent with the modified cell. The internal pressure \(\mathbf{P}_{\text{int}}\) defined in Equation (9) depends on the change in the target distribution \(\mathcal{P}_{\text{NST}}(\{\mathbf{q},\mathbf{p}\},\mathbf{h})\) when the lattice vectors \(\mathbf{h}\) are modified, and its calculation depended on how positions and momenta were affected by the scaling. Concerning positions, the standard recipe was to modify them using the same transformation matrix that was used to modify the lattice vectors, that is \(\mathbf{q}_{i}^{\prime}=\mathbf{h}^{\prime}\mathbf{h}^{-1}\mathbf{q}_{i}\), where \(\mathbf{h}^{\prime}\) denotes the cell matrix after the barostat had been propagated. This is equivalent to consider the change of the cell to be performed at constant scaled positions, \(\mathbf{s}_{i}=\mathbf{h}^{-1}\mathbf{q}_{i}\), and ensured that interparticle distances were modified in a way that did not depend on the periodic images in which the particles were located. Concerning momenta, in Ref. [18], we considered two possible strategies: (a) scaling momenta with an inverse factor than that used for positions and (b) leaving momenta unchanged. The first formulation should have been properly generalized to handle the anisotropic scaling used in this paper, where lattice vectors were multiplied by a deformation matrix rather than by a scalar. Here, we proposed to scale the momenta as \(\mathbf{p}_{i}^{\prime}=\left(\mathbf{h}(\mathbf{h}^{\prime})^{-1}\right)^{T} \mathbf{p}_{i}\), which corresponds to performing the cell deformation at constant scaled momenta defined as \(\boldsymbol{\pi}_{i}=\mathbf{h}^{T}\mathbf{p}_{i}\). Notably, momenta were scaled using the transpose of the inverse matrix employed for scaling the positions. If the matrix acting on the positions, \(\mathbf{h}^{\prime}\mathbf{h}^{-1}\), was a rotational matrix, it could be easily shown that positions and momenta would be identically rotated; thus, preserving the consistency between interatomic distances and interatomic velocities. More generally, it could be easily seen that scalar products between positions (or distances) and velocities (or relative velocities) were preserved with these scaling operations. Thus, for any transformation, if the relative velocity of two atoms connected by a constrained bond was orthogonal to the constrain before the transformation, it would also be so after the transformation. This might lead to better stability in presence of bond constraints. The definition of the internal pressure in this formulation is \(\mathbf{P}_{\text{int}}=\frac{\gamma}{\gamma}(\mathbf{K}-\Xi)\), where \(\mathbf{K}\) and \(\Xi\) are called kinetic energy tensor and _virial tensor_, respectively, and are defined as \(K_{\alpha\beta}=\sum_{i=1}^{N}\frac{p_{i}^{\prime}p_{i}^{\beta}}{2m_{i}}\) and \(\Xi_{\alpha\beta}=\frac{\gamma}{2}\sum_{i}h_{\beta i}\frac{\partial U}{ \partial h_{i}}\). The last derivative should be taken at constant scaled coordinates and might include contributions from long-range corrections [1]. In the second formulation, where momenta were not affected by the scaling, the internal pressure was defined as \(\mathbf{P}_{\text{int}}=\frac{\gamma}{\gamma}\left(\frac{N_{\alpha\beta} \mathbf{k}_{\beta}\mathbf{T}}{2}\mathbf{I}-\Xi\right)\). The idea of computing the internal pressure from the average kinetic energy rather than from the instantaneous one was proposed first in Ref. [26]. This formulation might have advantages related to the fact that it is less dependent on how the kinetic energy is computed [26; 27]. However, it is expected to be less effective when coupled with constraints or rigid bodies, since the consistency between positions and momenta would not be maintained. In addition, the evaluation of the internal pressure in presence of constraints is more difficult [18].