Pressure control in interfacial systems: Atomistic simulations of vapor nucleation

A large number of phenomena of scientiﬁc and technological interest involve multiple phases and occur at constant pressure of one of the two phases, e.g., the liquid phase in vapor nucleation. It is there-fore of great interest to be able to reproduce such conditions in atomistic simulations. Here we study how popular barostats, originally devised for homogeneous systems, behave when applied straightfor-wardly to heterogeneous systems. We focus on vapor nucleation from a super-heated Lennard-Jones liquid, studied via hybrid restrained Monte Carlo simulations. The results show a departure from the trends predicted for the case of constant liquid pressure, i.e., from the conditions of classical nucleation theory. Artifacts deriving from standard ( global ) barostats are shown to depend on the size of the simulation box. In particular, for Lennard-Jones liquid systems of 7000 and 13 500 atoms, at conditions typically found in the literature, we have estimated an error of 10–15 k B T on the free-energy barrier, corresponding to an error of 10 4 –10 6 s (cid:12) 1 (cid:27) (cid:12) 3 on the nucleation rate. A mechanical ( local ) barostat is proposed which heals the artifacts for the considered case of vapor nucleation. ' 2018 Author(s).


I. INTRODUCTION
Atomistic simulations are routinely used to investigate a variety of multiphase nanoscale systems, such as bubbles, drops, solid walls in contact with fluids, and solutions. In order to reproduce experimentally relevant conditions in small simulation samples far from the thermodynamic limit, barostats are needed to control the pressure.
The principle inspiring many barostats used in molecular dynamics (MD) is to generate the correct equilibrium distribution for the isothermal-isobaric or isoenthalpic-isobaric ensembles evolving an extended system of equations for the generalized degrees of freedom connected to the particles and simulation box. The force driving the expansion or compression of the system is the imbalance between the current instantaneous pressure, which depends on the positions and momenta of all the particles, and the target pressure. Also the dynamics of the particles are affected by the imbalance between the present and target pressure via the coupling with the degrees of freedom of the simulation box. [1][2][3][4] Because the instantaneous pressure depends on all particles, in the following we will refer to this class of barostats as global barostats. Global barostats are also used in Monte Carlo simulations. In this case, one typically alternates particles and volume moves. 5 The volume move is accepted or rejected depending on the instantaneous enthalpy (H + PV, where H is the Hamiltonian, P is the target a) simone.meloni@uniroma1.it b) alberto.giacomello@uniroma1.it pressure, and V is the volume of the sample) of the system before and after the move.
Pressure control is relevant also for the simulation of a variety of multi-phase systems, which is beyond the original scope of global barostats. What sets these systems apart is that different subdomains can have different pressures. A broad range of phenomena falls into this class, including homogeneous and heterogeneous vapor nucleation, [6][7][8][9][10][11][12] nucleation of polymorphic crystals, [13][14][15][16] dissolution of bubbles and droplets, and condensation or evaporation. In this work we show that in such cases, in which the relative amount of the two phases changes along the process, the pressure of the preexisting bulk metastable phase might change during the process when one uses global barostats, which is different from the condition at which experiments are carried out.
Here, in order to appraise these effects, we consider the case of vapor nucleation from a homogeneous metastable liquid. We present a simple macroscopic theory based on the sharp-interface model explaining the behavior of global barostats and their effects on nucleation. Atomistic simulations are performed for a Lennard-Jones (LJ) liquid in the same nominal thermodynamic conditions as those available in the literature 7,8 (both references use global barostats). A hybrid restrained Monte Carlo (hRMC) scheme 17,18 is adopted in order to cope with the problem of rare events 19 typical of nucleation and in order to compute the related free energy profile; the volume of the largest bubble is used as the order parameter. 20 The good agreement between macro-and microscopic results suggests that the intuitive argument of domains at different pressures is, indeed, at the origin of the artifacts associated with global barostats. A solution to these artifacts consists in using a local barostat that imposes the (local) force balance between a piston and the contacting liquid. Simulations are run using the local barostat showing that, at variance with global barostats, this approach is able to maintain the liquid pressure constant at the target value all along the process.
The manuscript is organized as follows. A macroscopic, sharp-interface model is introduced in Sec. II. In the same section, a microscopic formulation of the problem is presented. It is shown that within the sharp-interface limit the two representations are consistent. In Sec. III A, the simulation campaign is described in detail, while in Sec. III B we validate the local barostat for homogeneous systems. In Secs. III C and III D, the results are discussed. Section IV is left for conclusions.

II. THEORETICAL ANALYSIS OF CONTINUUM AND ATOMISTIC MODELS OF A TWO-PHASE LIQUID/VAPOR SYSTEM
We focus on the homogeneous nucleation of a vapor bubble in a metastable liquid. This deceptively simple case allows us to analyze the shortcomings of standard barostats in dealing with multiphase systems domains at different pressures. The same arguments should also apply to a variety of other multiphase systems, including heterogeneous vapor nucleation and condensation.
We start by introducing a simple continuum model of vapor nucleation-the sharp-interface-and the associated classical nucleation theory, CNT. 21 This model is based on a number of approximations, including the fact that the interface is ideally sharp, that are sometimes violated in actual systems. Nevertheless, within these approximations, it allows us to obtain an explicit dependence of the liquid pressure and of the energetics of the process on the volume of the vapor bubble, which helps understanding the shortcomings of standard (global) barostats. In the results, Sec. III C, we will illustrate that, even when the sharp-interface model approximations are violated, e.g., when the system is relatively close to the critical point, this theory captures the qualitative trend of the data.
In the sharp-interface model, it is assumed that the bulk properties of the fluids are valid up to the interface, where a sharp change in these properties occurs. The liquid and vapor domains are assumed to be uniform and isotropic. In particular, the diagonal terms of the stress tensor are all equal and the off-diagonal terms are zero (this hypothesis is consistent with the empirical observation of simulation data, Fig.  SM1 of the supplementary material). At the (infinitesimal) interface, these conditions are no longer met and the tangential and normal components of the stress tensor to the surface are different. 22 Within the sharp-interface model, this imbalance is translated into a surface tension γ acting at the dividing surface, which has an indirect influence on the liquid and vapor pressures via the (extended) Laplace equation. 10 In such a system, the average pressure of the whole sample reads where, consistently with the sharp-interface model, we assumed that the pressure field is of the form P(x) = P L θ L (x) functions of the liquid and vapor domains, respectively. 23 The interfacial terms do not contribute directly to the average pressure because the interface is sharp, i.e., it has an infinitesimal volume. V L and V V are the volumes of the liquid and vapor phases, respectively, V = V L + V V is the total volume, and χ V = V V /V and χ L = V L /V are the vapor and liquid volume fractions. An atomistic justification of Eq. (1) is given below. Equation (1) can be used to quantify the variation of the liquid pressure during an isothermal and isobaric bubble nucleation event. A closed set of equations for evaluating the liquid pressure can be obtained adding the extended Laplace law introduced in Ref. 10 or, if one is only interested in the liquid pressure at the critical nucleus, its conventional form valid for extremal points of the free energy. Here we use a simpler empirical approach: we assume that P V is constant and equal to the vapor tension at the simulated temperature; this approximation is then validated by atomistic simulations.
Conventional barostats used in atomistic simulations, 1,2,4 which have been designed for homogeneous systems, control the average pressure of the sample, P. Thus, within the sharp-interface model, the pressure of the liquid in a sample containing one vapor bubble of volume V V is where the dependence of the various terms on the volume of the bubble is made explicit. Since vapor nucleation occurs when P V > P L , Eq. (2) shows that the actual liquid pressure decreases along nucleation and that the driving force of the process, ∆P = P V P L = (P V P)/(1 χ V ), grows along it instead of remaining constant as it happens in actual experiments. Equation (2) can be used in conjunction with CNT to quantify the effect of conventional barostats on the free-energy profile of the process in a finite-size system. In CNT, where it is assumed that the pressure of the liquid is constant along the process (P 0 L ), the free energy difference between the liquid containing a bubble of volume V V and the reference bulk liquid reads 21 where N v is the number of vapor atoms in the bubble, µ V (P 0 L ) and µ L (P 0 L ) are the chemical potentials of the vapor and liquid phases at P 0 L , respectively, γ is the surface tension, and A is the area of the liquid/vapor interface. The second equality in Eq. (3) follows from a first order expansion of chemical potentials around the vapor tension P V . γ is assumed to be the planar surface tension of the two phases at coexistence.
Assuming that the liquid is incompressible and that, as said above, P V is constant and equal to the vapor tension, the free energy profile at variable liquid pressure is We remark that, owing to the many assumptions of CNT, Eq. (4) does not necessarily describe in quantitative terms atomistic results, but it is certainly useful to explain what are the potential artifacts connected with the use of conventional barostats on the free energy profile. In Sec. III A, atomistic simulations implementing various methods for controlling the pressure will be used to quantify these effects on the free-energy profile and nucleation barrier.
(2)], the effect of conventional barostats is that of reducing the barrier [Eq. (4)] as compared to the case of constant liquid pressure [Eq. (3)]. In Fig. 1, we report both the free-energy profile ∆G 0 (V V ) according to Eq. (3) (black line) and the free-energy profile ∆G(V V ) according to Eq. (4) (red and blue lines for systems of 7000 and 13 500 particles, respectively).
In Eqs. (3) and (4), the free energy is computed setting P 0 L = 0.026, P V = 0.046, and γ = 0.098 for the reference liquid pressure, vapor tension, and surface tension at T = 0.855. 7 [Lennard-Jones units are used throughout the article: temperature, pressure, length, and time are reported in reduced units, /k B , /σ 3 , σ, and σ(m/ ) 1/2 , respectively]. The liquid volume is assumed to be constant during nucleation and consistent with the bulk density of atomistic systems of N = 7000 and 13 500 particles: V L = N/ρ L , where ρ L = 0.58 is the metastable liquid density at the current pressure and temperature of simulations. Given the difference between the liquid and vapor densities (ρ V = 0.08), this approximation has a minor effect on free energy. Global barostat free energy profiles are shown in Fig. 1 and compared with CNT results. This comparison shows that the free energy profiles with 7000 and 13 500 are below the CNT one; in particular, the nucleation barrier ∆G † , i.e., the difference between the maximum and initial free energy, follows the trend ∆G † 7000 < ∆G † 13500 < ∆G † CNT . Indeed, this is consistent with the observation that the driving force ∆P grows along nucleation for global barostats, and its growth is more marked for the smaller sample.
In order to extend these results to more general systems, it is worth estimating the error affecting the free energy barrier as a function of the size of the system and thermodynamic conditions due to the global barostat. Figure 2(a) reports the size of the sample corresponding to an error on the barrier of 10 k B T as a function of ∆P = P 0 L − P V . As expected, the closer the system is to two-phase coexistence (∆P = 0), the larger the critical nucleus is, and the larger must be the sample to keep the error under the prescribed threshold. Our model suggests that in the physical conditions studied in previous work 7,8 and in the present work, the atomistic system should contain at least ≥10 4 particles in order to have an error on the free-energy barrier ≤10 k B T. Panel b of the same figure presents the percent error on the free-energy barrier as a function of the ratio between the total volume and the volume of the critical bubble, V/V † V . The continuum sharp-interface model shows that, independently of the thermodynamic conditions, simulation boxes 15 times larger than the critical bubble are necessary to have errors on the barrier ≤10%.
The microscopic expression for the pressure of an isotropic system consisting of n particles interacting via a pair potential is where p i and m i are the momentum and the mass of the ith particle, f ij is the force between ith and jth particles, and r ij is their (vector) distance. If we consider a two-phase system containing n L bulk liquid, n V bulk vapor particles, and n int interface particles, we can rewrite Eq. (5) as the sum of three terms, associated with the liquid, vapor, and interface domains, When the interface is vanishingly small, the contribution of the corresponding term is negligible and the pressure of the sample is expressed as the sum of the first and second term. When the liquid and vapor domains are large enough, these terms can be interpreted as the liquid and vapor pressures, 24 and Thus, consistently with the macroscopic sharp-interface model in Eq. (1), when the interface thickness is negligible, Eq. (6) reduces to the volume-weighted average of the liquid and vapor

A. Simulation details
We considered a system composed of particles interacting via the truncated and force shifted (TFS) Lennard-Jones (LJ) potential, analogous to those considered in Refs. 7 and 8, where with a cut-off radius r c = 2.5. In the TFS-LJ potential, the pair particle forces go to zero smoothly as r goes to r c . The liquid vapor phase diagram of the TFS-LJ system has been reported in Refs. 7 and 25. We compute the vapor nucleation free-energy barrier as a function of the largest vapor bubble in the system, V V , estimated using the M-method. 20 The method consists of several steps. (i) Particles are labeled as liquid-like if they have more than five particles closer than 1.6 σ, and vapor-like otherwise.
(ii) The simulation box is partitioned into cells. The size of the cells is chosen such that they can contain at most one particle. A cell is labeled liquid or vapor if it contains a liquid-like or vapor-like particle. Empty cells are classified analysing both the first and second neighbors cells: If the number of nearest neighbor face-sharing empty/vapor cells is 7 or more, also the number of second nearest neighbor face-sharing empty/vapor is evaluated; if also the number of these cells is 7 or more, the original empty cell is labeled as vapor. (iii) Finally, a cluster analysis is performed on the vapor cells and the size of the largest bubble is established as the total volume of largest cluster of interconnected cells, i.e., cells sharing a face or a corner (Fig. 3).
To study vapor nucleation, we employ the hybrid Restrained Monte Carlo (hRMC) approach, 17,18,26 which is well suited for non-analytical collective variables (CV), such as the size of the largest vapor bubble used here. hRMC allows to sample the conditional probability density function at the current value of the volume of the vapor bubble, and to compute conditional averages. Thus, one can estimate the mean force by the conditional average of the observable −k( 17,27 which can be numerically integrated to obtain the free energy profile along the nucleation process. An in-depth explanation of the hRMC method is given in the Appendix.

hRMC with a global barostat
A typical MC method for sampling constant pressure ensembles consists in alternating particles and volume moves. Particle moves are accepted or rejected according to the Metropolis criterion, which will be detailed below for the case of hRMC. In volume moves, a random, isotropic expansion/compression is generated and particle positions are rescaled accordingly. The move is accepted or rejected on the basis of the energy and PV values before and after the move.
In the first step, a short NVE MD simulation is integrated, starting from the previous configuration and with the momenta extracted from a Maxwell-Boltzmann distribution at the relevant temperature. The acceptance probability is where H and H are the extended Hamiltonian of the system before and after the move, respectively. The extended Hamiltonian is the sum of kinetic, K(p), and (physical) potential, U(r), energies plus a biasing potential energy term which forces the system to visit configurations in which V V fluctuates around the target value V * V : k is the coupling constant determining the degree of fluctuations allowed to the volume of the bubble (see the Appendix for more details). The second MC step consists in a change of the volume of the system.
where H and H and V and V are the extended Hamiltonians and volumes of the system before and after the move, respectively, P is the target pressure, and N is the number of particles.

hRMC with the local barostat
To overcome the artifacts due to global barostats, we also adopt a local barostat, which consists in enclosing the system between two moving walls of particles to which a constant additional force f is applied [ Fig. 4(a)]. The wall particles interact with the fluid via a suitable potential (here LJ) and, at stationarity, the total force F = fn Wall exerted on the liquid by the n Wall particles is equal and opposite to that exerted by the fluid particles on the walls, i.e., when the external pressure F/A, with A the area of the walls, is equal to the liquid one, P L . Thus, with the present barostat, stationarity is determined by the (local) balance between the forces of the piston and of the liquid in contact with it rather than on the average pressure of the sample, including vapor domains.
In the present work, each wall is made of two layers of TFS-LJ atoms (50 times heavier than the fluid ones) in the fcc lattice configuration. The LJ parameters are WW = 10 WF and FF = WF (W = wall, F = fluid). In Fig. 4(b), the calibration curve P L vs f is reported for a bulk TFS-LJ liquid; this graph shows that the macroscopic prediction P L = fn W /A is fulfilled, confirming the mechanical balance mechanism by which the local barostat controls the liquid pressure. Figure 4(b) reports data obtained with different values of WF and σ WF indicating that the local barostat does not sensitively depend on the chosen solid-liquid interaction potential. In other words, the local barostat is rather robust and does not require fine tuning of the solid-liquid interaction.
Other local barostats can also be adopted, e.g., that are based on a non-interacting particles gas, 28 but we found the moving walls one to be simpler to use in the presence of a gas phase.
The hRMC simulation protocol used to implement the local barostat is the following. A short MD NVE trajectory of both fluid and solid particles is integrated, initializing particles momenta from a Maxwell-Boltzmann distribution. The boundary conditions are free in the direction orthogonal to the walls and periodic in the other directions. The acceptance probability reads In this case, the extended Hamiltonian isH(p, r, where the sum runs over the 2 × n W particles of the moving walls and z i is their position in the direction orthogonal to the walls.
Before closing this section, it is worth mentioning that the use of the local barostat is not limited to simple atomic fluids. For example, one can use the local barostat also with molecular fluids such as water. We show this by simulating a small box of TIP4P/Ew 29 water with two pistons, in which the wall particles interact with the oxygen atoms of water molecules

B. Validation of the local barostat
We validated the local barostat by comparing results against those obtained with a global one for bulk systems. In particular, we focused on the distribution of instantaneous pressures and on the phase diagram (Fig. 6). One notices that  the instantaneous pressure distribution obtained with the local barostat is, within the error bars, the same as that obtained with the global one. Also the liquid and vapor branches of the TFS-LJ binodal obtained with the local barostat match very well with literature data. 25 We also considered the case of more complex molecular fluids, by comparing the liquid branch of the binodal of TIP4P/Ew water obtained by the local and global barostats and literature data 31 [Fig. 6(c)]. Also in this case, there is a very good matching of local barostat results with reference data.

C. Vapor bubble nucleation
Simulations of vapor nucleation are performed at T = 0.885 and P = 0.026, i.e., the same conditions used in the literature. 7,8 We considered two computational samples containing 7000 and 13 500 particles. These samples are relatively large, in particular, the second system is larger than those used in the literature. [6][7][8] For each sample, we computed the freeenergy profile vs bubble volume with both the global and the local barostats. The mean forces are estimated at a set of 20 values of the bubble volume V V of the largest vapor bubble in the sample (see the Appendix and Refs. 17 and 27).
As a first remark, we notice that results obtained with the local barostat for the two samples of different size are in good agreement between them (Fig. 7) and with the CNT predictions (Fig. 1). The barrier and critical size are slightly smaller in the atomistic case; this effect is well known (see, e.g., Ref. 6) and is associated with the limits of the continuum model, namely, to the idealized sharp-interface. With the global barostat, the system shows a significant dependence of the free-energy profile on the sample size. In particular, the barriers are (22 ± 1) k B T and (30 ± 1) k B T for the small and large samples, respectively, both significantly smaller than the value measured with the local barostat, (40 ± 1) k B T and (39 ± 1) k B T for the small and large samples, respectively. These results confirm that, in order to have an accurate prediction of the nucleation barrier, free of finite size effects arising from the pressure control, one has either to simulate very large samples or to resort to a local barostat.
The errors in the free-energy barriers are reflected with exponential sensitivity on the nucleation rates, which are one of the final goals of the simulations of nucleation. Assuming that the nucleation rate follows a CNT-like relation, k = k 0 exp(∆G † /k B T ), and assuming that the kinetic prefactor k 0 is not affected by how pressure is controlled, one estimates differences of 4-6 orders of magnitude between the local and global barostat rates, depending on the size of the sample. Even larger errors are expected in the case of fewer particles often used in the older literature.

D. Effect of the barostat on the properties of the liquid and vapor domains
The sharp-interface interpretation of the effect of the global barostat on the free energy, discussed in Sec. II, is that the liquid pressure decreases along vapor nucleation. Here we investigate the variation of the pressure of the liquid domain and other properties of the system as a function of the vapor bubble size with the global and local barostats to validate the theoretical predictions.

Density
We start by analyzing the dependence of the (conditional ensemble averaged) radial density field, ρ(r; V V ) (r is the distance from the center of the bubble), on the type of barostat. ρ(r; V V ) has been computed for both the 7000 and 13 500 particles samples and with both barostats in a radial range encompassing the bubble, interface, and liquid domains. We considered samples containing bubbles of several sizes, from very small to supercritical ones. Very small bubbles, V V ≤ 700, do not present well defined vapor domains. For bubbles larger than this threshold [ Figs. 8(a) and 8(b), V V = 1500, to be compared with a critical nucleus of V † V ∼ 2500], the radial density presents the expected profile with bulk vapor and liquid domains separated by an interface. The first observation is that with both barostats and for both samples the interface, the region in which the density changes from low (vapor) to high (liquid) values, is rather thick, ∼8. This large value is not surprising considering that simulations are performed at pressure and temperature conditions relatively close to the critical point.
A second observation is that there are important differences between the radial densities obtained with the two barostats. With the local barostat, the density field of both samples shows two plateaus at small and large r [see the insets of Fig. 8(a)], corresponding to the vapor and liquid domains, respectively. The density in the bubble is very close to the value corresponding to the vapor tension, which confirms the reliability of the approximation on the value of P V used in Sec. II. At the other end of the radial range, the density in the bulk liquid reaches the expected value. With the global barostat, on the contrary, in the smaller system, the radial density does not reach the vapor and liquid plateaus. In particular, the value of the radial density at the last point is 3.5% lower than the liquid bulk value at the target pressure and temperature. In the large sample, the radial density reaches the target liquid density value but the curve presents a significant slope in this domain, which suggests that it does not correspond to the bulk liquid. This is confirmed by independent NVT simulations performed at the average density of the last four points of ρ(r; V V ), in which we measured the total vapor fraction, χ tot V , i.e., the vapor fraction due to all bubbles present in the liquid [ Fig. 8(c)]. Our results show that the system presents two regimes: for densities close to the bulk value, the one measured in the liquid domain of samples containing a small nucleating bubble, χ tot V is small and constant; for densities corresponding to samples containing larger nucleating bubble, χ tot V is large and grows with V V , i.e., with decreasing ρ. 32 This confirms that, with the global barostat and in the presence of critical bubbles, the liquid does not behave as a bulk liquid. Concerning ρ(r; V V ) at small r, in the vapor region the radial density is slightly above the target value.
We believe that the remarkable effect of the global barostat on the density has two main reasons: (i) the relatively large compressibility of the LJ liquid and (ii) the thick interface at the present thermodynamic conditions. We expect that for less compressible liquids, e.g., water, and at thermodynamic conditions further from the critical point, the effect of the global barostat on the density would be smaller. This does not mean that in these cases the barostat-related artifacts on the energetics of nucleation would be smaller, simply it might be more difficult to identify that simulations are performed with an inappropriate setup.

Liquid pressure
It is important to evaluate the pressure of the liquid domain in order to validate the assumptions behind the effect of global barostat. In Fig. 9, we report the pressure of a liquid control volume far from the vapor bubble and from the solid walls computed via Eq. (7) with the prescriptions of Irving and Kirkwood. 33 These results show the expected decreasing trend of P L with the bubble size. However, since the sub-domains are small, the large statistical error of the estimated pressure makes it difficult to draw reliable conclusions. This is especially critical in samples containing larger bubbles in which the limited bulk liquid domain imposes to use very small control volumes. Thus, we also follow a different approach, which consists in first determining the mean density in the liquid domain, which converges with the number of hRMC steps faster than the local pressure, and then computing the pressure via an independent NVT simulation of a bulk liquid with 3000 particles at the measured density. The density of the bulk sample is set to the average density of the last four points of the radial profile for J. Chem. Phys. 148, 064706 (2018) FIG. 9. Comparison between the pressure computed via Eq. (7) and via an NVT simulation at the density of the bulk liquid domain in a system containing a bubble of volume V v . The solid black line represents the target value of P L . These results show that the two approaches are equivalent but the former has a much larger statistical error associated with it. selected values of the bubble volume [ Fig. 8(c)]. Results show that pressures estimated with both methods are consistent, with lower errors connected with the second one (Fig. 9).
Despite the improved statistical accuracy, due to the relatively large scattering of the density values (see the insets of Fig. 8), estimated with the second approach is also limited the overall accuracy of pressure estimated with the second approach is also limited. Thus, one should focus on the FIG. 10. Liquid pressure as a function of the bubble size. The red and blue symbols represent the pressure controlled by a global barostat for the samples of 7000 and 13 500 particles, respectively; purple and green points refer to the pressure for samples controlled by the local barostat. The red and blue lines are the continuum predictions for the liquid pressure [Eq. (2)]. The black line represents the target liquid pressure. In the figure, we also report the colormaps of the density field of two snapshots of the samples with 7000 (upper panel) and 13 500 (lower panel) particles at V V ∼ 2500. These snapshots show that the departure of the pressure from the target value is due to the interaction of the thick interfaces with their periodic images. This problem for bubbles close to the critical size has already been put forward by Meadley and Escobedo 8 for their simulations on a sample of 10 000 particles at the same thermodynamic conditions. When a bubble interacts with its periodic image, the radial density in the liquid domain (Fig. 8) used to compute the pressure is reduced and the pressure decreases. qualitative effects of barostats on the P L vs V V curves. With the local barostat, the liquid pressure is almost constant all along the process and very close to the target value, typically within the statistical error from the reference pressure (Fig. 10). For samples containing larger bubbles, one observes a small reduction of liquid pressure, which is related to the overlap of the bubble with its periodic images that lowers the "liquid" density. On the contrary, in the case of global barostat, the pressure significantly decreases with the bubble size. This occurs with both samples but the phenomenon is enhanced in the case of 7000 particles. With the large sample, the liquid pressure is initially close to the target value and then deviates for V V ≥ 1000.
The dependence of the pressure on the bubble volume and, for a given V V , on the number of particles in the sample is consistent with the analysis of Sec. II. However, atomistic simulations show a larger deviation from the target pressure than that predicted by the sharp-interface model. We believe that this is due to two reasons: (i) the limited accuracy in the estimate of the pressure via the density of the liquid domain 34 and (ii) the presence of a very thick interface, which is not taken into account in the sharp-interface model, i.e., in Eq. (1) one discards both (a) the continuous change of the normal pressure in going from the liquid to the vapor domain and (b) the tangential contribution, which differs from the normal one. Nevertheless, it is remarkable that even in conditions very far from those of Sec. II the theoretical predictions are in qualitative (pressure) and quantitative (nucleation barrier) agreement with atomistic results.

IV. CONCLUSIONS
In this work, we have addressed the issue of controlling pressure in vapor nucleation from a metastable liquid. Our theoretical analysis and numerical simulations show that global barostats result in an underestimation of the liquid pressure, which is particularly severe far from two-phase coexistence. In turn, this can bring artifacts on the driving force and, ultimately, on the free energy of the process.
According to our analysis based on the sharp-interface model, in order to have an error on the nucleation barrier ≤10% in a simulation in which the global pressure is set equal to a target value, the simulation box volume should be around 15 times larger than the critical bubble volume.
To confirm the theoretical predictions, we have performed hRMC simulations aimed at computing the free-energy profile along the nucleation pathway. Atomistic data show a qualitative agreement with the theoretical predictions.
Our results suggest that simulations using standard barostats, [6][7][8]35 if not performed on a reasonably large system size, might be affected by errors of the order of 10-15 k B T on the barrier height, corresponding to an error of 10 4 -10 6 s 1 σ 3 on the rate.
Finite-size effects associated with pressure control can be eliminated by replacing conventional, global barostats, developed for single-phase systems, with a local barostat, which controls the pressure of the liquid phase. This can be achieved by adding moving solid walls interacting with the liquid via, e.g., a Lennard-Jones potential. The walls, to which an external force is applied, act as pistons compressing the liquid at the desired pressure all along the nucleation process and make it possible to perform constant liquid pressure simulations even with small simulation boxes (e.g., 7000 particles in the present system).
To conclude, to have an accurate prediction of the nucleation barrier, free of finite-size effects arising from the pressure control, one has either to simulate significantly larger samples or to resort to a local barostat.

SUPPLEMENTARY MATERIAL
See supplementary material for the figure with the components of the stress matrix in 30 control volumes in a bulk (LJ) liquid using the global and local barostats.

APPENDIX: RESTRAINED MONTE CARLO
In our simulations, we estimate the free energy barriers using the hRMC approach presented in Refs. 17 and 18. In this method, the atoms are subject to the extended potential U(r) + U k (r). Here U(r) is the TFS-LJ interaction, U k (r) = k/2(V V (r) − V * V ) 2 is the biasing term where V V (r) is the current volume of the largest bubble in the system estimated with the M-method 20 and V * V is the target value of the bubble volume.
Following Ref. 19, we show how the free energy can be reconstructed from restrained simulations. Consider the average where Z k (V * V ) ≡ ∫ dr exp − βU k (r, V * V ) and Z = ∫ dr exp − βU(r) is the canonical partition function. Since Z is independent of V * V , it was introduced in the second equality in (A1) in order to interpret f k (V * V ) as the derivative of Recalling that the Landau free energy of a variable is defined as F(V * V ) = − β −1 ln P V V (V * V ), we find that in the proper limit Eq. (A1) is an estimate of the derivative of the free energy, ∇ V * V F(V * V ). The mean force (A1) can be estimated using hRMC and the relative free energy via integration. The conventional approach of MC, in which a single particle is subjected to a random displacement, makes simulations inefficient since the order parameter, which in this case is computed through the expensive procedure described in the main text, must be evaluated at each step. Therefore, in order to increase the efficiency, we use the hybrid Monte Carlo (hMC) approach in which at each time step the displacement of a single atom is replaced by a collective displacement according to a short MD trajectory. MD is started from the current particle configuration, while momenta are extracted from a Maxwell-Boltzmann distribution at the given temperature. Differently from standard MC, the acceptance criterion takes into account also the initial and final kinetic energy. In hMC, the Hamiltonian which generates the MD dynamics could be different from the one adopted in the acceptance test. 36,37 Since the dependence on r in V V (r) is non-analytical, here we choose to use the physical potential U(r) to generate the dynamics and the restrained potential U k (r) for the MC acceptance.