## 4 Discussion  
In this work, we developed the fully anisotropic version of the stochastic cell rescaling (SCR) barostat that was introduced in Ref. [18]. The algorithm was then implemented in three different MD software packages and tested on a number of crystalline systems. When possible, results were compared with alternative algorithms already available in the tested codes. The isotropic and semi-isotropic versions introduced in Ref. [18] could be obtained as constrained cases of the equations of motion developed in this work. Additionally, this method could be shown to be equivalent to a Parrinello-Rahman barostat augmented with friction and noise terms, provided that the high-friction and zero-mass limit is taken. The anisotropic SCR barostat is based on a first-order stochastic differential equation that describes the relaxation of the lattice vectors to their equilibrium values. The deterministic part of the equation resembles the anisotropic Berendsen barostat [17]. A crucial advantage of the anisotropic SCR barostat with respect to the Berendsen barostat is that, thanks to the noise term, the correct cell fluctuations are sampled. A single input parameter should be provided, \(\tau_{P}\), that allows to enforce an approximate value for the relaxation time of the volume. For this approximation to hold, one should also provide an a priori estimation of the isothermal compressibility of the system \(\beta_{T}\). In practice, only the value of \(\frac{\beta_{T}}{\tau_{P}}\) is used in the resulting equations of motion. Hence, an incorrect estimate of the isothermal compressibility would result in an observed relaxation time different from \(\tau_{P}\). This issue is common with the Berendsen barostat [17]. Nevertheless, the analyses carried out in this work showed that the method is robust against variations of \(\tau_{p}\) within 2-3 orders of magnitude, suggesting no significant effects on the sampled \(NST\) distribution if the input \(\beta_{T}\) is estimated with the correct order of magnitude. Two integrators were developed and tested. One of them is simpler and could be straightforwardly integrated in existing MD codes. The other one is more complicated as it requires propagating position and velocities simultaneously to lattice vectors, and was only implemented in the SimpleMD code. The second integrator has the advantage of being time-reversible and, thus, allowing detailed balance violations to be quantified. Figure 10: Results from the simulations of the gold crystal systems in LAMMPS, employing the MTTK and SCR anisotropic barostats in presence of a linearly increasing shear stress. The plot shows the trajectory for the modulus of the cell vector \(\mathbf{b}\), which is a good variable to describe the crystal breaking. Selected structures are shown with the corresponding sampling times along the horizontal axis. In our algorithm, velocities could be optionally scaled when the lattice vectors were modified. We propose to use a scaling matrix for velocities that is the inverse transpose of the matrix used to scale positions. This is different from the standard procedure used in Ref. [8; 9; 16], where it can be shown that both positions and velocities were rescaled with exponential matrices in the limit of small \(\Delta t\), and these matrices were the inverse of each other. To the best of our knowledge, this idea is new, and allows rotations to be handled consistently on position and velocities, which is expected to be more effective when simulating systems with holonomic constraints. This advantage was present when global rotations were removed by keeping the cell matrix triangular, since, in this case, the transformation matrix was not symmetric. It is, instead, irrelevant when implementing the MTTK algorithm constraining the cell momentum matrix to be symmetric [8], since, in this case, the transformation matrix was also symmetric. Indeed, this approach has been argued to be more effective in simulating systems with constraints [1].