{'fraction_non_alphanumeric': 0.033948420303345564, 'fraction_numerical': 0.008979904725401083, 'mean_word_length': 4.3771529515677905, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Isothermal-isobaric molecular dynamics using stochastic velocity rescaling

The authors present a new molecular dynamics algorithm for sampling the isothermal-isobaric ensemble. In this approach the velocities of all particles and volume degrees of freedom are rescaled by a properly chosen random factor. The technical aspects concerning the derivation of the integration scheme and the conservation laws are discussed in detail. The efficiency of the barostat is examined in Lennard-Jones solid and liquid near the triple point and compared with the deterministic Nos\'{e}-Hoover and the stochastic Langevin methods. In particular, the dependence of the sampling efficiency on the choice of the thermostat and barostat relaxation times is systematically analyzed.


I. INTRODUCTION
Starting from the breakthrough article of Andersen 1 many schemes have been proposed to modify the Newton equations of motion so as to perform molecular dynamics (MD) in the isothermal-isobaric ensemble, where the number of particles, the external pressure and the external temperature (N P T ) are fixed. Most of these schemes are based on an extended-Lagrangian formulation. Before proceeding to overview the details of the extended system method, it is better to separate two issues: pressure and temperature control. Within Andersen's framework the control of pressure was achieved by adding to the Hamiltonian auxiliary dynamical variables, which represented the volume V and its conjugate momentum. The resulting configurations were sampled from the N P H ensemble, where the enthalpy H was constant. This approach has been further extended to allow anisotropic variation of the cell shape and size. 2,3 Several almost-equivalent schemes have been proposed, with similar performance and similar parameters, often in combination with different thermostats (see e.g. Refs. 4,5).
In the Andersen work, the temperature was regulated by additional stochastic collisions, leading to a sampling of the N P T ensemble. The idea underlying the Andersen's thermostat was thus similar to the use of Monte Carlo, and inherited from Monte Carlo the characteristic weak points such as discontinuity of the trajectories, lack of a conserved quantity, at least in the original formulation, and artifacts in the computation of dynamical properties. Another intrinsic disadvantage of the algorithm due to its stochastic nature was the not-reproducibility of the trajectory, unless the same random number generator was employed.
In subsequent refinements, the Nosé-Hoover deterministic thermostat 4,6,7 replaced the stochastic one. For this purpose, a new auxiliary variable was introduced, playing the role of a friction on the particles and/or on the cell dynamics and controlling the temperature. This method allowed for the definition of a conserved quantity, crucial to check the integration accuracy. However, a weak point was that the time-evolution of the kinetic energy was described by a second-order differential equation, resulting in an algorithm which is inefficient for equilibration. Even worst, the lack of ergodicity provided a poor control on the temperature in difficult cases such as harmonic solids and, in general, systems with normal modes. 8 Thermostat chains suggested in Ref. 9 partially solved these problems, however at the price of a more complex formulation and a larger number of parameters.
A possible alternative to the deterministic thermostat is stochastic Langevin dynamics. 10 Its combination with extended-system control of the pressure, the so-called Langevin piston scheme, 11 was proven to be efficient and ergodic even in difficult cases. The lack of a conserved quantity has been a traditional drawback of stochastic MD, which was however recently solved. 12 However, it is worth noting that the choice of the friction drastically influences the dynamics and, indirectly, the sampling efficiency. On the one hand, a too small friction coefficient may cause a very long equilibration time. On the other hand, a too large friction coefficient may induce strong perturbations of the system dynamics hindering the exploration of the phase space.
A further possibility for temperature and pressure control is the deterministic scheme suggested by Berendsen et al. 13 This scheme is not based on an extended Lagrangian, and the volume of the simulation box is driven by the difference between the internal pressure and the target one. Similarly, the kinetic energy is driven by its deviation from the target one. The algorithm is intuitive, easy to implement and efficient in the equilibration phase. Unfortunately it cannot be used when the correct distribution is required (e.g. when fluctuations have to be calculated), since it does not sample the N P T ensemble. Moreover, in difficult cases the Berendsen scheme systematically transfers energy to the slowest degrees of freedom and can lead to wrong results. 11,14 Hence the "common wisdom" suggests that the initial equilibration should be carried out with the Berendsen scheme, while the ensemble averages should be calculated from the subsequent runs performed with an extended-system method.
In a recent paper 15 we proposed a thermostat which can be regarded as a stochastic version of the Berendsen scheme. Similarly to the Nosé-Hoover scheme, 4,6,7 it is a global thermostat, in the sense that it acts only on the total kinetic energy of the system, and produces configurations in the canonical ensemble. In addition, this scheme allows defining a conserved quantity, does not suffer of ergodicity problems in solids 15 and has been used successfully in practical applications for equilibration purposes 16 or to perform ensemble averages. 17,18 Finally, it has a minimal impact on the system dynamics, thus being an optimal choice for the computation of dynamical properties. 19 The object of the present article is to combine the global stochastic thermostat introduced in Ref. 15 with a standard barostat 1,4 to sample the isothermal-isobaric ensemble. The properties of the resulting algorithm are discussed in details and illustrated on Lennard-Jones (LJ) liquid and solid systems. A special attention is dedicated to the parameterization of the thermostat and barostat relaxation times. In particular, we provide a framework for the choice of the barostat relaxation time which is based on the conservation of the total effective enthalpy, thus on the accuracy of the resulting sampling. For the thermostat relaxation time, we discuss its impact on the sampling efficiency by measuring the autocorrelation times of several relevant physical quantities. In this sense, we compare the new scheme with the Nosé-Hoover 4,6,7 and with the Langevin piston barostats. 11 This paper is organized as follows: In Section II we formulate the most important aspects of the theoretical framework linking statistical theory to equations of motion. In particular, we demonstrate that the stochastic rescaling algorithm samples the isothermic-isobaric distribution. Section III contains several applications of our scheme to LJ fluid and solid systems. Some comments on the practical implementation of the barostat and a reliable discretization scheme are collected in the Appendix.

II. THEORY
A. Isothermal-isobaric ensemble We consider a system of N particles, described by coordinates q i , momenta p i , masses m i , contained in a box of variable volume V subject to an external pressure P ext . With q and p we indicate the set of coordinates q i and p i , respectively. In the isotropic N P T ensemble the probability distribution is usually defined as 20,21 where K(p) = i |p i | 2 /(2m i ) is the kinetic energy, U (q, V ) is the potential energy and β = (k B T ) −1 , with k B the Boltzmann constant and T the temperature. The expression for P N P T has been addressed in the literature and is not unique. For instance, Attard 22 proposed to use the logarithm of the volume as the integration variable, or, equivalently, an extra V −1 term in the probability distribution. Nevertheless we will adopt the more standard definition in Eq. (1).
When MD simulations are performed in absence of external forces (i.e. j ∂U/∂q j = 0) the momentum of the center of mass of the system is conserved. Thus, it is convenient to adopt as dynamical variables the centerof-mass coordinate q cm and momentum p cm , and the coordinates r and momenta π relative to the center of mass After the change of variables, the N P T ensemble reads where K(π) = i |π i | 2 /(2m i ) is the kinetic energy relative to the center of mass and K cm (p cm ) = |p cm | 2 /(2 i m i ) is the kinetic energy of the center of mass. The delta-functions in Eq. (3) are related to the conservation of total momenta and inertia in Eqs. (2) and to the fact that, in the center-of-mass framework, j m j r j = 0 and j π j = 0. As a further simplification, the redundant variables q cm and p cm and can be integrated out, resulting in the distribution Here the additional V term comes from the fact that the size of the domain on which q cm is integrated is proportional to V . The contribution of the extra V factor is almost negligible for systems with more than a few tens of particles, though we will include it to have a formally correct ensemble.
The distribution in Eq. (4) provides ensemble averages which are equivalent to those from Eq. (1), and should be used as a target ensemble for any MD simulation in the N P T ensemble when external forces are not present. In the next two Subsections we provide an algorithm to sample this distribution.

B. Pressure control
We first consider the problem of fixing the pressure, adopting a solution similar to that of Refs. 1,4,5. Since we work in the reference frame of the center of mass, we choose an initial configuration satisfying i π i = 0 and i m i r i = 0. Then we evolve it according to the following equations of motioṅ Here f i = −∂U/∂r i are the forces, η is proportional to the relative change rate of the volume and W is the inertia of the piston, which determines the typical time scale for volume changes. The instantaneous internal pressure P int is given by where the derivative d/dV is performed at fixed scaled coordinates, as discussed in more detail in the Appendix. Equations (5) are similar to those introduced by Hoover, 4 but for an additional term in Eq. (5c). The original Hoover formulation, as observed by the author himself, gives a slightly wrong ensemble. With our correction, when the dynamics is combined with a thermostat as shown in the next Subsection, the exact N P T distribution is sampled. An alternative solution to this problem have been proposed by Martyna et al. 5 . It can be easily shown that the following quantities are conserved during the dynamics: The first two conservation laws come from the absence of external forces. The last quantity H is close to the enthalpy of the system (K +U +P ext V ), and its initial value H 0 depends on the initial configuration and velocities. When an ensemble of systems is evolved according to Eqs. (5), their distribution evolves according to a corresponding Liouville-like equation. Defining the (6N + 2)dimensional vector x = (r, π, V, η), Eqs. (5) can be written in compact form asẋ = G(x), and the Liouville-like equation on the density asρ = −∂ x (Gρ) = −Lρ, wherê L is the Liouville operator. Written explicitly in terms of the phase-space variables r, π, V and η, the Liouville operator readŝ The long time distribution for the ensemble can be found solving the equationLP = 0. As it can be verified by direct substitution, the following class of functions solveŝ LP N P H = 0 The delta functions here come from the conservation laws in Eqs. (7). Assuming that the system is ergodic an MD simulation based on Eqs. (5) will produce configurations according to this distribution. An alternative derivation of Eq. (9) can be done introducing the phase-space compressibility as it is done for instance in Ref. 23. Rigorously speaking, Eq. (9) is not the standard N P H distribution. Indeed, here the constant of motion is H, which, as discussed above, slightly deviates from the enthalpy. Moreover, there is an additional V −1 prefactor. However, since we are interested in N P T sampling, we ignore these discrepancies and proceed in the combination of the N P H scheme with a thermostat.

C. Temperature control
We now combine the presented constant-enthalpy scheme with a thermostat. We first extend the N P T ensemble of Eq. (4) so as to include η as a dynamical variable. Since η is a velocity, we associate to it a kinetic energy and define the extended N P T ensemble as where K * (π, η) = K(π) + W η 2 2 is the total kinetic energy, which includes also the barostat kinetic energy. When the η variable is integrated out, this distribution is identical to the target one in Eq. (4). Thus, we will use this distribution as a target ensemble for the N P T algorithm. It is easy to verify that this distribution is stationary with respect to the Liouville-like operator in Eq. (8), i.e.LP N P T (π, r, V, η) = 0. The only reason why a time average performed with Eqs. (5) will not give results consistent with the N P T ensemble is that Eqs. (5) have an extra constant of motion (namely H). To sample the N P T ensemble one has to introduce a thermostat, which changes at every step the value of H in such a way that the N P T ensemble is stationary. Since P N P T is the product of a term depending on the potential energy and a term depending on kinetic energy, the thermostat can be designed to act only on the latter one, whose distribution in the Gibbs ensemble is known a priori and is Here N * f = 3N − 2 is the number of degrees of freedom, including the volume (+1) and excluding the center of mass (−3). In the following, we describe how this task is accomplished by the stochastic velocity rescaling. 15 In the stochastic velocity rescaling, the application of the thermostat consists of a rescaling of the velocities of the system. We here rescale not only the momenta π but also the barostat velocity η, by a factor α: This change affects the total kinetic energy K * = K + W η 2 2 , which includes also the barostat kinetic energy. Following Ref. 15, we calculate at each step the rescaling factor by propagating the total kinetic energy according to the equation The propagation is done for a half timestep ∆t/2 as described in the Appendix, so that α 2 = K * (t + ∆t/2)/K * (t). Here τ T is the relaxation time of the thermostat,K * = N * f /(2β) is the average kinetic energy, and dW is a Wiener noise. 24 At variance with the scheme in Ref. 15, here the cell degree of freedom needs to be counted, and the rescaling procedure affects both the velocities of the particle and the velocity of the piston η.
As it has been shown previously, 15 the stationary distribution associated with the stochastic dynamics in Eq. (12)  Our algorithm consists of an alternation of constantenthalpy steps [Eq. (5)] and a thermostat step [Eq. (12)], as discussed in details in the Appendix. It is worthwhile noting that other schemes can be implemented in the same fashion, simply by changing the thermostat step. In particular, the Langevin piston algorithm 11 is easily implemented using the thermostat step as discussed in Refs. 12,25, applying the friction and noise to the particles and the cell velocities, and removing the velocity of the center of mass afterwards so as to eliminate the total force on the center of mass and obtain the same ensemble. Similarly, a standard Nosé-Hoover scheme may be coupled with the constant-enthalpy integrator.

D. Control of sampling errors
In the practical implementation of MD one integrates Eqs. (5) using a finite timestep. This invariably introduces small systematic errors in the sampled distribution.
The traditional way to control this error is to perform a microcanonical simulation and to verify the conservation of the total energy. This approach is easily extended to variable cell simulations in the N P H ensemble, where the conservation of the enthalpy is used. An equivalent conserved quantity can also be defined for simulations performed with the Nosé-Hoover scheme in the N P T ensemble. In a recent paper we have introduced a formalism which allows to define a conserved quantity, called effective enthalpy, also for stochastic dynamics. This has been done for our stochastic velocity rescaling, 15 and has been later extended to Langevin dynamics. 12 Here we discuss how this concept is generalized to an effective enthalpy in variable cell simulations. To this aim, we first repeat some of the ideas introduced in the previous papers, then we show how the compressibility is taken into account for our variable cell algorithm.
We consider the (6N + 2)-dimensional vector x = (p, r, η, V ) in the phase space. In principle MD equations would produce a continuous trajectory x(t). However, in the practical implementation the time is discretized so that only a discrete sequence of snapshots of the system is available, x 0 , x 1 , . . . , which is then used to evaluate ensemble averages. In our scheme, each point x i is obtained from the previous one x i−1 applying alternatively Eq. (5) and the thermostat. We now consider the move from x i to x i+1 as a proposal for a Monte Carlo move. Even if in an MD simulation all the moves are accepted, it is instructive to calculate the acceptance rate for that move. To enforce detailed balance relative to the target distribution P N P T (x), the acceptance rate in the Metropolis scheme should be set equal to The transition matrix M (x ′ ← x) is the probability of generating x ′ as the next point in the sequence given that the present point is x. The star denotes time reversal of velocities, i.e. if x = (p, r, η, V ) then x * = (−p, r, −η, V ). Clearly, P N P T (x) = P N P T (x * ). The detailed balance enforced by Eq. (13) is not the standard one but its generalization for odd and even variables discussed in Ref. 24. We then associate an effective enthalpyH i to each point in the sequence in such a way that the acceptance rate in Eq. (13) is equal to min 1, e −β(Hi+1−Hi) . By simple substitution, the change of effective enthalpy between subsequent points is defined as We now proceed to an explicit evaluation of the increment of effective enthalpy at each MD step. Our scheme is a sequence of two different operations: the thermostat and the time evolution of Eq. (5). When the thermostat is applied, either stochastic rescaling or Langevin, its propagation is done analytically, so that detailed balance is exactly satisfied. This gives an acceptance equal to 1 and a vanishing change in the effective enthalpy. 15 On the other hand, when Eqs. (5) are propagated, the change in effective enthalpy is not trivial and has to be evaluated explicitly. The change in effective enthalpy is decomposed in two terms as .
(15) The first term, substituting the ensemble probability in Eq. (10) and the definition of enthalpy in Eq. (7c), is where V i is the volume at step i. The second term of Eq. (15) is calculated reminding that the rate between the backward transition matrix and the forward transition matrix is equal to the Jacobian of the transform which, as shown in the appendix, is (17) Combining Eqs. (15), (16) and (17) one obtains: so that the increment in effective enthalpy is exactly equal to the increment in the enthalpy.
In conclusion, the calculation of the effective enthalpy is very similar to the calculation of the effective energy in Ref. 15. The effective enthalpy is obtained as a cumulated sum of all the increments of H in the constant-enthalpy dynamics. Alternatively, one can obtain the effective enthalpy as the total enthalpy H minus the sum of all the contributions to the kinetic energy given to the thermostat. In the limit of zero time step, the effective enthalpy is exactly conserved. In real applications, the effective enthalpy is expected to exhibit random fluctuations and a small systematic drift. In principle, if exact sampling is required, the change in effective enthalpy can be used to perform hybrid Monte Carlo simulations, 26 or to reweight properly the obtained configurations. 27 In practice, one can just monitor its systematic drift and use it to quantify the sampling errors. If the drift is too large, one should decrease the integration time step.

III. EXAMPLES
We now test the discussed algorithm on a couple of model systems. We consider a system of 256 particles interacting through LJ potential in a cubic box with periodic boundary conditions. Simulations are carried out both in the solid fcc phase and in the liquid one at temperature β −1 = 0.692 and pressure P ext = 0, close to the triple point. We use everywhere LJ reduced units. The interactions are truncated at 2.5, and the force is smoothed between 2.25 and 2.5. Long-range corrections which take into account the difference between the true LJ potential and our truncated approximation are added to the total energy of the system. To this aim, a constant radial-distribution-function is assumed for distances larger than 2.25, 20 and the corrections to the internal pressure are evaluated as derivatives of the energy correction with respect to the cell volume. The dynamics is computed using an in-house code.
For convenience, the piston mass W is defined as N f β −1 τ 2 P , where τ P is a time describing the time-scale of the volume fluctuations. 7 For the Nosé-Hoover scheme discussed in Section III, the thermostat mass is chosen as Q = N f β −1 τ 2 T , where the thermostat relaxation time τ T measures the thermostat strength.
In the following, three kinds of thermostats are considered: (a) our scaling procedure (stochastic rescaling), which is stochastic and global; (b) Langevin piston, which is stochastic and local; (c) Nosé-Hoover, which is deterministic and global. This choice allows us to investigate separately which is the effect of a global thermostat versus a local one and how the stochastic nature of the algorithm affects its performance. We investigate the role played by the relevant simulation parameters, namely the timestep ∆t, thermostat relaxation time τ T , the barostat relaxation time τ P . For the Langevin scheme, we choose a friction γ = (2τ T ) −1 , for the reasons discussed in Ref. 19. The presented averages are obtained from 10 7 steps simulations for each possible choice of the simulation parameters.

A. Timestep, barostat mass and effective-enthalpy drift
We consider here the choice of the timestep ∆t and of the barostat relaxation time τ P . As it will become clear later, the rationale behind the choice of these two parameters is the same, thus we discuss them at the same time.
In Fig. 1 we show the dependence of the effectiveenthalpy drift dH dt on ∆t, fixing the other parameters at τ P = 0.5 and τ T = 0.2. The customary procedure in MD simulations is to choose a ∆t which is as large as possible, with the constraint that the drift of the total energy is smaller than a given threshold. In our case, an equivalent check can be done on the effective enthalpy introduced in Subsection II D or, for the Nosé-Hoover, on its conserved quantity. The exact value of the threshold is rather arbitrary, but still it gives a quantitative estimate of the sampling errors given by the finite timestep, and still can be used to compare the relative accuracy of two different simulations. The two global schemes (Nosé-Hoover and stochastic rescaling) exhibit a similar drift, which grows with ∆t. Also for the Langevin scheme the drift increases with the timestep, but faster than in the global schemes. For further calculations, we choose a  We now consider the choice of τ P . We recall that there is a relationship between the mass of the particles (and of the barostat) and the time-step. As an example, a fourfold decrease of all the masses produces exactly the same trajectory as a two-fold increase of the time-step, thus giving a two times faster exploration. In the same way, for a fixed value of the particles mass, a lower barostat mass will provide a faster volume evolution at the price of a lower accuracy. Thus, the choice of the barostat mass, W = N f β −1 τ 2 P , is based on a compromise between a small mass and a reasonable effective-enthalpy drift, and is very similar to the choice of ∆t. In Fig. 2 we show the dependence of the effective-enthalpy drift on τ P . Since the drift induced by a too small τ P is due to sampling errors of a single degree of freedom (the volume), it is partially hidden by the larger drift due to the inexact integration of the particles degrees of freedom. However, it is clear that for τ P ≥ 0.5 the drift due to the volume is negligible. We notice that in none of the presented simulations we have observed an appreciable decoupling between the barostat and the internal degrees of freedom in the case of small τ P . We also computed the autocorrelation time of a few relevant quantities, and we noticed that the sampling efficiency always increases as τ P is decreased. Thus, the final rule turns out to be a τ P which is as small as possible, but with a reasonable effective-enthalpy conservation. We thus choose τ P = 0.5 to continue our tests. We now proceed to a systematic study of the effects of the other relevant simulation parameter, namely the thermostat relaxation time τ T . To this aim, we fix a timestep ∆t = 0.005 and a barostat relaxation time τ P = 0.5, and perform simulations for a wide range of τ T . As it is seen in Fig. 3, the effect of τ T on the effective-enthalpy drift is rather weak, at least for the two global schemes.
Indeed, for these choices of the parameters, the origin of the drift is not in the thermostat but in the constantenthalpy step. The drift in Langevin is larger than that of the other schemes. For τ T we cannot use the same criterion that we used to optimize τ P . Indeed the two stochastic schemes (rescaling and Langevin) are stable for any possible choice of τ T . This is due to the fact that, in these two cases, the thermostat equations are linear and can be integrated in an analytic fashion (see Appendix). This is not true for the NH scheme, which turns out to be be unstable for τ T ≈ 0.02, at least using the present integration scheme. Moreover, it is not obvious that a smaller τ T should lead to faster sampling, as it will be seen in the following.
We thus proceed in an analysis of the statistical efficiency as a function of the thermostat parameter, similar to what has been done for the N V T ensemble. 19 A quantitative measure of the time needed for performing ensemble averages of a given quantity X is given by its autocorrelation time τ X (see e.g. Ref. 28). Since a common goal of N P T simulation is the calculation of averages of quantities like the enthalpy and the volume of the system, we consider here the cases where X = H and X = V . The autocorrelation time of the total enthalpy also provides a rough estimate of the time needed for the equilibration procedure. However, it must be noticed that to achieve a fast equilibration also the autocorrelation time of the fluctuations of the enthalpy should be short. Thus, we also calculate the correlation time for the fluctuations of these quantities, obtained with X = ∆H 2 ≡ (H − H ) 2 and X = ∆V 2 ≡ (V − V ) 2 respectively. The first one is related to the bulk modulus and the second to the heat capacity at constant pressure. Thus, these autocorrelation times provide the statistical efficiency in the calculation of these important physical observables.
The simulation results for a number of runs are summarized in Fig. 4 and Fig. 5 for the liquid and the solid respectively. In all cases, the performance of stochastic rescaling is equal or better than that of the local Langevin thermostat. For large τ T , the two stochastic schemes have the same behavior. For small τ T , the performance of the local Langevin scheme is worse. This happens for both the solid and liquid systems, and indicates that tuning the performance of Langevin thermostat is very difficult since it is efficient only for a limited friction range. These problems in Langevin thermostat are related to the fact that it is local, thus it disturbs considerably the trajectory, hindering diffusion (in the liquid) or thermalization of the slow modes (in the solid). This is in agreement to what has been observed in Ref. 19 for a liquid in the N V T ensemble.
On the other hand, the Nosé-Hoover thermostat is global and, similarly to the velocity rescaling, has a very small impact on the trajectory. In the present calculations, performed at the triple point, we don't expect the solid to be harmonic enough to introduce ergodicity problems. Its performances are excellent when looking  at the autocorrelation time of enthalpy (τ H ) and volume (τ V ), for both the solid and the liquid. For the solid, a resonance effect leads to a peak in these autocorrelation times for τ T ≈ 2. However, fluctuations of the enthalpy (τ ∆H 2 ) are well thermalized only for a limited range of choices for τ T . This is related to the fact that the thermostat is optimally coupled only to the modes with a given frequency, and is not able to thermalize the others. The fact that the average enthalpy has a short autocorre-lation time is due to cancellations in the autocorrelation function, and is a signature of the enthalpy oscillations around the correct value. The damping of these oscillations can be rather slow if τ T is not chosen properly.
The stochastic rescaling scheme combines the stability of the stochastic Langevin with the efficiency of the global Nosé-Hoover. The coupling parameter τ T in the stochastic rescaling scheme has the meaning of thermostatstrength, and a smaller τ T always guarantees a better efficiency. On the contrary, in the Nosé-Hoover scheme τ T also determines which modes are optimally coupled with the thermostat, and it is more difficult to tune. Finally, the stochastic rescaling is comparable with the Nosé-Hoover in computation of averages for small τ T , surpassing it in computation of fluctuations.

C. Dynamical properties
It is also interesting to check if the rescaling scheme preserves the dynamics. To investigate this point, we consider the calculation of dynamical properties. For the LJ liquid, we calculate the self-diffusion coefficient D from the mean square displacement. In order to overcome the practical problem of computing the mean square displacement in a variable cell, first we calculate it using the scaled coordinates V −1/3 q i , then we scale the obtained diffusion coefficient of a factor V 2/3 to recover the correct value. The diffusion coefficient as a function of the thermostat relaxation time τ T is shown on Fig. 6. It is evident that both global schemes (velocity rescaling and Nose-Hoover) do not affect the diffusion. In particular, the calculated value (D = 0.03) is in agreement with that obtained in microcanonical simulations in similar conditions. On the other hand, when the local Langevin scheme is used with a too high friction the dif-fusion is hindered. Remarkably, the choices of τ T which lead to a slower diffusion correspond to the choices that lead to a worse sampling (see the previous Subsection). These time can also be compared with the typical collision time, defined as the velocity autocorrelation time τ v . It is easy to show that τ v is related to D by τ v = βmD. In the limit of large τ T (small friction), it turns out to be τ v = 0.04. When the collisions with the thermal bath have the same frequency that the natural interparticle collisions, the thermostat turns out to be the bottleneck for diffusion. In the solid, where the diffusion coefficient is vanishing, we opt for the calculation of the full velocity-velocity autocorrelation function. The results are shown in Fig. 7, where the autocorrelation function is plotted for all the considered thermostats, using a reasonable thermostat relaxation time τ T = 0.1, and for a reference microcanonical simulation at similar conditions. From the plot it is clear that the global schemes do not affect significantly the dynamical properties. Even if these simulations are performed in the N P T ensemble, they can be used to calculate correct dynamical properties. This is not true for the local Langevin scheme, where the autocorrelation function is strongly affected by the thermostat. This is an additional confirmation that the impact on the dynamical properties is not due to the stochastic nature but to the locality of the Langevin thermostat.

IV. CONCLUSIONS
A new method has been introduced to perform MD simulations at constant temperature and pressure. The method is stochastic but does not affect the dynamical properties. The concept of effective energy, which provides a constant of motion for stochastic algorithms and which was already introduced in previous papers, 12,15 has been extended in the present scheme to an effective enthalpy. A systematic analysis of the role played by the control parameter has been shown, focused on the choice of the thermostat and barostat relaxation times. Moreover, a comparison of the new scheme with standard Nosé-Hover and Langevin barostats has been shown. The new scheme appear to be simpler to use and easier to control. Further work is required to apply the presented algorithm to the case of nonequilibrium MD, where the Gibbs' ensemble does not hold. 29

APPENDIX A: INTEGRATION SCHEME
At variance with Ref. 15, we follow here a timereversible scheme which is consistent with the Trotter decomposition scheme. 30,31 The Trotter decomposition is usually adopted to decompose the Liouville (or Fokker-Planck) operator which describes the dynamics of an ensemble of realizations.
For instance, if the exact solution of the Liouville equation for a finite time increment ∆t is ρ(t + ∆t) = e −(L1+L2+L3)∆t ρ(t), its approximated solution is ρ(t + ∆t) ≈ e −L1 ∆t 2 e −L2 ∆t 2 e −L3∆t e −L2 ∆t 2 e −L1 ∆t 2 ρ(t). Equivalently, one can state that if the equations of motion are in the formẋ = G 1 (x) + G 2 (x) + G 3 (x) an approximate propagator can be obtained applying in the proper order the analytical solutions of the equationsẋ = G 1 (x), x = G 2 (x) andẋ = G 3 (x). A necessary condition is that the three decoupled equations can be solved exactly.
We here decompose the equations of motion as a sum of three parts: the first is the propagation of the thermostat, and the others are two equations of motion whose sum is equivalent to the N P H dynamics. The flow of the integration scheme is as follows: 1. Propagate the thermostat for a time ∆t/2.
2. Propagate velocities for a time ∆t/2, according tȯ 3. Propagate positions and velocities for a time ∆t, according toṙ 4. Propagate again velocities for a time ∆t/2, as in step 2.
5. Propagate again the thermostat for a time ∆t/2, as in step 1.
The sum of Eqs. (A1) and (A2) is equivalent to the original equations of motion (5). This choice for the splitting is arbitrary, but it is motivated by the fact that Eqs. (A1) and (A2) can be solved analytically also for finite time.
We now give explicit expressions for the propagation of the single stages. Stages 2 and 4 are propagated solving Eqs. (A1): The instantaneous pressure here is calculated as where f ij (t) is the force between particles i and j and r ij is their distance, evaluated by taking into account the periodic boundary conditions within the minimumimage convention. 20,32 The last term ∂U/∂V takes into account any explicit dependence of the potential energy on the volume, such as the long range corrections. 5,20 Stage 3 is propagated solving Eqs. (A2): r i (t + ∆t) = sinh[η(t)∆t] η(t) p i (t) m i + e η(t)∆t r i (t) (A5a) V (t + ∆t) = e 3η(t)∆t V (t) (A5b) p i (t + ∆t) = e −η(t)∆t p i (t).
It is also useful to evaluate the compressibility associated to this propagators, defined as the Jacobian of the corresponding transforms. The transforms in Eqs. (A3) and (A5) are linear, with Jacobian equal to respectively 1 and e 3η(t)∆t = V (t + ∆t)/V (t). Thus, when a volume element in the phase space is transformed by Eqs. (A3) and (A5) its measure is changed by a a factor V (t + ∆t)/V (t). This contribution is crucial for a proper evaluation of detailed-balance violations as it is done in Subsection II D.
We now describe the propagation of the thermostat stages (1 and 5). The thermostat can be either velocity rescaling, Langevin or Nosé-Hoover. We here simply