## 1 Introduction  
Molecular dynamics (MD) simulations are commonly used to generate trajectories at atomistic resolution for molecular systems or materials. Typically, they take advantage of thermostats and barostats so as to control the temperature and pressure and to generate conformations representative of well-defined statistical ensembles [1]. Pressure is often controlled in an isotropic manner [2], so that the shape of the simulation cell is fixed and its volume is modified by means of a uniform scaling of the lattice vectors. This is necessary, for instance, when simulating liquid systems or solvated molecules, which do not offer any resistance to shear. However, when simulating solids, it is convenient to let the cell shape change so as to allow the simulated system to relax to its optimal geometry, to explicitly simulate phase transitions, and to enable the calculation of stress-strain curves. The idea of performing MD simulations with flexible cells was pioneered by [PERSON] and [PERSON] [3, 4], who introduced an extended Lagrangian formalism to propagate the nine variables associated with the generating vectors of the periodic lattice. The cell is, thus, provided with an inertia, and its acceleration is controlled by the difference between the internal pressure and the external stress, resulting in a second-order differential equation for the components of the lattice vectors. Monte Carlo algorithms can also be used to sample cell deformations [5, 6], although their application has been limited. A number of variants of the original Parrinello-Rahman method has been published [7, 8, 9, 10, 11, 12], all of them based on a second-order differential equation. The possibility to add friction and noise terms [13, 14, 15] or to combine MD with Monte Carlo [16] has been discussed. These formulations, however, retained the second-order behavior of the algorithms discussed above. An important drawback of second-order equations is that, if parameters are not tuned properly, the volume and shape of the simulation cell might be subject to fluctuationsthat are difficult to dampen. A notable exception is the first-order Berendsen barostat [17] in its anisotropic form, which is usually considered effective in the equilibration phase, but does not sample any known ensemble. A common approach is, thus, to first equilibrate a system with the Berendsen barostat and then simulate it with a second-order barostat.  
In a recent paper, we introduced a first-order algorithm named stochastic cell rescaling (SCR), which resembles the isotropic Berendsen barostat but additionally includes a stochastic contribution, resulting in the correct isobaric distribution [18]. A semi-isotropic variant, suitable for simulating constant tension systems such as membranes and interfaces, was also developed and tested in Ref. [18]. Here, we develop the formalism for the completely anisotropic version, where all the cell components are allowed to change and external stress can be controlled. The isotropic and semi-isotropic variants can be seen as constrained cases of the anisotropic scheme introduced here. The algorithm is implemented in several MD codes (the educational MD software SimpleMD, GROMACS [19], and LAMMPS [20]) and tested on a number of systems (Lennard-Jones solids, ice, gypsum, and gold), where it is also compared with commonly employed alternatives, namely, the Parrinello-Rahman (PR) [3] and the Martyna-Tuckerman-Tobias-Klein (MTTK) [8; 9] algorithms.