Monte Carlo simulations in the unconstrained ensemble

The unconstrained ensemble describes completely open systems whose control parameters are chemical potential, pressure, and temperature. For macroscopic systems with short-range interactions, thermodynamics prevents the simultaneous use of these intensive variables as control parameters, because they are not independent and cannot account for the system size. When the range of the interactions is comparable with the size of the system, however, these variables are not truly intensive and may become independent, so equilibrium states defined by the values of these parameters may exist. Here, we derive a Monte Carlo algorithm for the unconstrained ensemble and show that simulations can be performed using chemical potential, pressure, and temperature as control parameters. We illustrate the algorithm by applying it to physical systems where either the system has long-range interactions or is confined by external conditions. The method opens up a new avenue for the simulation of completely open systems exchanging heat, work, and matter with the environment.

The Metropolis Monte Carlo (MC) method [1] vastly contributed to the understanding of many physical phenomena [2,3]. Different versions of the method have been devised, applied to various statistical ensembles as, e.g., microcanonical [4], canonical [1], grand canonical [5], semi-grand canonical [6], isothermal-isobaric [7,8], isostress-isostrain [9], and distinct variants of the Gibbs ensemble [10,11]. All these ensembles include at least one extensive variable as a control parameter, such as energy, volume or number of particles: little attention has been paid to MC methods in which the control parameters are the chemical potential µ, pressure P , and temperature T . Indeed, thermodynamics tells us that the latter intensive quantities are not independent and cannot account for the size of a macroscopic system with short-range interactions [2,12]. When applied to these systems, the µP T ensemble requires the addition of the equation of state linking µ, P , and T or, conversely, it can be used to infer such a link [13]. Reported MC methods [14,15] have taken advantage of this fact by considering a constrained µP T ensemble in which the examined systems are partially closed in either volume or number of particles, thus finding the underlying relation µ(P, T ) of the system under scrutiny.
In contrast, an unconstrained ensemble with µ, P , and T as independent control parameters can be properly defined for small systems [16], confined systems [17,18], and for long-range interacting systems [19,20], which are intrinsically nonadditive [21,22] and have an additional degree of freedom that may render µ, P , and T independent. This makes it possible to study completely open systems exchanging heat, work, and matter with the environment. In this Letter we derive an elementary MC scheme for simulations in the unconstrained ensemble and show that it consistently combines the MC algorithms of the grand canonical and isothermal-isobaric ensembles. By testing the method for simple physical systems that can be analytically evaluated, we identify the role of interactions in equilibrium states of completely open systems.
Since dealing with completely open systems is rather unconventional, we first recall some concepts about their thermostatistics [16,19]. Considering for simplicity a one-component system where the combination of the first and second laws of thermodynamics is expressed by dE = T dS −P dV +µdN , where E, S, V , and N are the energy, entropy, volume, and number of particles, respectively, one defines a statistical ensemble by taking a collection of N independent replicas of the system, with total energy, entropy, volume, and number of particles E t = N E, S t = N S, V t = N V , and N t = N N , respectively. The energy balance for the ensemble becomes [16] where the last term accounts for the variation of E t when N varies at constant S t , V t , and N t ; the quantity E is called subdivision potential [16] or replica energy [22,23].
No assumption has been made on the nature of the system, so Eq. (1) showing that the natural variables of the replica energy are µ, P , and T . One realizes that conventional thermodynamics [12] focuses on systems in which E = 0 by imposition, but this is not the general situation. Systems with E = 0 are nonadditive [22], because their entropy is not a linear homogeneous function of E, V , and N . As discussed below, E can be derived [16] from a partition function Υ(µ, P, T ) such that E (µ, P, To establish the basis for the MC algorithm in the unconstrained ensemble, we closely follow Ref. [2], extending the standard arguments used for other ensembles. For simplicity we consider a system in a cubic box of side L = V 1/3 (the extension to a rectangular box is immediate) in which particles have coordinates r i , i = 1, . . . , N , and define scaled coordinates by r i = V 1/3 s i . The canonical partition function of the system takes the form (2) where s N ≡ (s 1 , . . . , s N ), β = 1/k B T , Λ is the de Broglie thermal wavelength, and U(s N ; V ) is the potential energy. We assume that the system is coupled to two independent reservoirs of ideal gas particles at temperature T : reservoir a that exchanges particles with the system and reservoir b that exchanges volume (see Fig. 1). Reservoir a has N a − N particles and volume V a , reservoir b has N b particles and volume V b − V , and their canonical partition functions Q a (N a − N, V a , T ) and Q b (N b , V b − V, T ) are easily obtained from Eq. (2) by setting U = 0. The partition function of the total system including the reservoirs is Q tot = Q a QQ b . Thus, the probability density of observing the system with N particles and volume V is 0 dV Q tot , so the most probable values of N and V are those that minimize the total free energy F tot = −k B T ln Q tot . Taking the limit of infinite reservoirs, the quantities that survive, besides T , are the chemical potential µ of reservoir a and the pressure P of reservoir b (see the Supplemental Material [24] for details). Hence, the probability density of finding the system in volume V in a particular N -particle configuration takes the form where the unconstrained partition function is given by (4) Notice that it would not be possible to implement completely open conditions with a single reservoir of ideal gas particles, since µ, P , and T are not independent in this case.
Given a system configuration C from which a new configuration C is generated in the simulation, we follow the Metropolis algorithm [2] using Eq. (3) and compute the acceptance probability of the new configuration as P acc (C → C ) = min [1, P(C )/P(C)]. MC moves in this case consist of displacements of particles, insertion and removal of particles, and changes of volume, yielding a potential energy variation ∆U = U(C )−U(C). A particle displacement is attempted by selecting a single particle at random with coordinates s and performing a random displacement from s to s . According to Eq. (3) and the Metropolis rule, this move is accepted with a probability P acc (s → s ) = min 1, e −β∆U . Similarly, the insertion of a particle at a random position and the removal of a random particle are accepted with respective probabilities Finally, trial moves that attempt to change the volume from V to V are accepted with probability Equations (5) and (6) correspond to particle insertion and removal acceptance probabilities in the grand canonical ensemble, while Eq. (7) is the acceptance probability for volume changes in the isothermal-isobaric ensemble [2]. Therefore, a consistent MC algorithm for simulations in the unconstrained ensemble can be obtained as a simple combination of the usual algorithms for these two ensembles. In the Supplemental Material [24] we describe the details of the procedure adopted to decide which kind of move is implemented at each MC step. We now illustrate two applications of the method. As a first example, consider a system with spatially constant, repulsive interactions φ(r i , r j ) = ε, for which the potential energy is , where ε > 0 is a coupling constant. This system is nonadditive because interactions are long-ranged, regardless of its size. Although interactions not depending on the interparticle distances may be difficult to justify physically, with this simple model we will show that repulsive interactions in nonadditive systems can withstand completely open conditions resulting in states of thermodynamic equilibrium. For the present case, the canonical parti- , from which the unconstrained partition function (4) becomes Since ε > 0, the probability of observing the system with N particles, given by the exponential in Eq. (8), has a maximum at with the control parameters taken such thatN > 0. Introducing T * = k B T /ε and x = N/T * , we rewrite the partition function (8) Hence, the distribution has a sharp peak aroundx in the limit T * → ∞, requiringN → ∞ withx fixed. Since E = −k B T ln Υ, we have E = − 1 2 εN 2 in this limit, which displays the dependence of the replica energy on µ, P , and T through Eq. (9). From this relation, it follows that the average volumeV = ∂E /∂P satisfies PV =N k B T , just like for an ideal gas. The basic feature introduced by the interactions, however, is that the size of the system can be controlled with independent µ, P , and T , which is impossible for a macroscopic ideal gas. Moreover, we have shown that a largeN can be realized for temperatures k B T ε, which can be seen as a weak coupling limit for small ε. In order to get a finite densityN 3 /V = P * /T * in the thermodynamic limit, the reduced pressure P * = P 3 /ε must be large as well, where is an arbitrary unit of length. Accordingly, the reduced chemical potential µ * = µ/ε − T * ln(Λ 3 / 3 ) has to be taken of the same order as T * and P * , since, from Eq. (9), we see thatN = µ * − T * ln(P * /T * ) + 1/2.
The average number of particles is known analytically in this simple model. To test the proposed MC scheme, we perform simulations in the unconstrained ensemble for this system taking µ * , P * , and T * as control parameters. In this case, the potential energy remains constant for random particle displacements and variations of volume, while ∆U = εN and ∆U = ε(1 − N ) for the insertion and removal of a particle, respectively. Particle displacements are only rejected when the particle leaves the simulation box (we do not use periodic boundary conditions here) and the remaining acceptance probabilities can be easily obtained from Eqs. (5)-(7), whose explicit expressions are given in the Supplemental Material [24]. The simulations are shown in Fig. 2 for fixed µ * while P * is varied, and varying µ * with fixed P * , in both cases for different values of T * . Solid curves in this figure correspond to the analytical expression forN .
Despite the fact that completely open conditions can be achieved by assuming the system coupled to two reservoirs, controlling µ, P , and T independently may be challenging in practice. In what follows we describe an example in which completely open conditions can be readily conceived and which shares key similarities with the previous example, even though this is not evident at first sight. We consider a system of N hard spheres of diameter σ at temperature T confined between parallel plates of area A and gap thickness H, so the volume of the system is V = AH. The dimensions of the plates are constant and much larger than the particle diameter σ, while H is comparable to σ. Neglecting hard-core interactions between particles but keeping A(H − σ) as the available volume due to particle-plate interactions, the system behaves as an ideal gas with free energy F given by [25] βF In this situation one can define a lateral pressure P lat = −H −1 ∂F/∂A (acting in directions parallel to the plates) and a transversal pressure P = −A −1 ∂F/∂H (acting in the direction perpendicular to the plates), which can be shown to be related to each other through [25] P = HP lat /(H − σ).
Clearly, P → P lat as H → ∞, so they are no longer independent in this limit. In addition, the chemical potential is given by µ = ∂F/∂N and can be written as Thus, the canonical replica energy When H and σ are comparable, the system is nonadditive because E = 0. In the limit H → ∞, P lat AH = P V → N k B T , so E → 0 and additivity is recovered. Assume now that, in this idealized approximation, the canonical and unconstrained ensembles are equivalent and let H and N be fluctuating quantities with aver-agesH andN , respectively. Ensemble equivalence does not hold, in general, when interactions are included. We suppose that the plates are surrounded by a fluid, acting as a reservoir in equilibrium with the system, with temperature T and chemical potential µ r which therefore fixes its pressure P r (as in a narrow pore with slit geometry [17,26]). Since the system is laterally open, equilibrium requires that µ = µ r and P lat = P r , while the transversal pressure P can be externally controlled regardless of the value of P r , provided the gap thickness is not too large. One possible way to control P is by applying weights to the plates, so these exert on the system a transversal pressure larger than the pressure P r of the surrounding fluid. Thus, µ, P , and T are independent and define the state of the system. A reason for this is that the pressure P is not an intensive quantity, since it depends on the size of the system as can be seen from Eq. (10). Although the above arguments presume an ideal system, they capture in first approximation the behavior of hard-sphere systems including the interactions, as we test below with simulations in the unconstrained ensemble.
We now concentrate on a more realistic treatment of the confined hard-sphere system by means of the MC method. In the simulations, lengths are measured in units of σ and energies in units of k B T . Accordingly, we introduce the reduced pressure P * = βP σ 3 and chemical potential µ * = βµ − ln(Λ 3 /σ 3 ), while the temperature is just a scaling factor which will be assumed constant. Periodic boundary conditions are implemented at the limits of the box in the transversal directions, while hard walls are assumed in the direction perpendicular to the plates. All MC moves in the unconstrained ensemble, consisting of particle displacements, insertion and removal of particles, and variations of the box length H at constant A, are rejected if they lead to an overlap between particles, between particles and the plates, and between the plates. When there is no overlap, particle displacements are always accepted and the remaining acceptance probabilities can be obtained by direct substitution from Eqs. (5)-(7) taking V = AH and a vanishing potential energy (see the Supplemental Material [24]). Scaling of particles coordinates for volume changes are performed in the direction perpendicular to the plates only. Here we further consider that the fluid surrounding the plates consists of hard spheres described by the Carnahan-Starling equation of state [27], for which the pressure is P * r (η r ) = (6η r /π)(1 + η r + η 2 r − η 3 r )/(1 − η r ) 3 in our reduced units, where η r is the packing fraction in this fluid. The corresponding chemical potential [28] reads µ * r (η r ) = ln(6η r /π) + (8η r − 9η 2 r + 3η 3 r )/(1 − η r ) 3 . We emphasize that this characterization of the reservoir serves only to evaluate the simulations, since the actual packing fraction of the system is determined by the control parameters µ * and P * .
To be able to compare the simulations with the analytical approximation discussed above, we restrict ourselves to small η r . We first choose η r = 0.04, 0.08, and 0.12 and fix the chemical potential of the system to µ * = µ * r (η r ), meanwhile an external pressure P * is applied on the plates (their area is fixed to A/σ 2 = 20 2 ). In Figs. 3(a) and 3(b), we show the average gap thick-nessH and the packing fractionη = π 6 σ 3 N/V as a function of P * . We observe thatH becomes large and approaches a macroscopic limit when the applied pressure approaches the pressure P * r (η r ) of the reservoir. In this limit, the freedom to control µ * and P * independently at fixed temperature is lost, as expected. Solid lines in Fig. 3(a) correspond to the approximation P * =HP * r (η r )/(H − σ) given by Eq. (10). Furthermore, a smooth kink aroundH ≈ 2σ is observed at P * = P * c ≈ 0.54 for the case η r = 0.12, which is enlarged in the inset of Fig. 3(a). For decreasing pressures P * > P * c , the gap slowly increases and the system is quasi-bidimensional because configurations with two particles aligned in the direction perpendicular to the plates are not realized. Fluctuations of the gap size at P * c allow such configurations to be realized and the gap grows faster for decreasing pressures P * < P * c . The average packing fraction consequently decreases for decreasing pressures P * > P * c . This behavior is also observed for η r = 0.04 and 0.08, but it is less pronounced in the plots. We next consider a situation in which P * is fixed to a value P * = P * r (η 0 ) for some packing fraction η 0 . Then, the chemical potential of the reservoir is varied (by changing η r ), so µ * = µ * r (η r ) also changes. In other words, we control µ * keeping P * constant. The results of the simulations are shown in Figs. 3(c) and 3(d) for η 0 = 0.04, 0.08, and 0.12. When µ * approaches µ * r (η 0 ), the average thicknessH approaches the macro-scopic limit. Solid lines in Fig. 3(c) correspond to the approximation µ * = µ * r (η 0 ) − ln H /(H − σ) , obtained by combining Eqs. (10) and (11). Therefore, we have shown that chemical potential and pressure can be independently controlled in this system at fixed temperature. We emphasize that the hard-core repulsion makes the system nonadditive for gap sizes comparable to σ, since E = 0 in this case. Similar to the previous example, a key point here is that repulsive interactions maintain equilibrium states under completely open conditions in a regime in which the system is nonadditive.
When restricting to the particular situation of a confined geometry as in the example above, we highlight that the MC scheme presented here is similar to the approach introduced in Ref. [17] describing the grand isostress ensemble. In that case, the replica energy E plays the role of the grand isostress potential considered there.
In conclusion, we have shown that a consistent MC scheme for simulations in the unconstrained ensemble can be obtained by combining the algorithms of the grand canonical and isothermal-isobaric ensembles. This scheme applies to nonadditive systems in which chemical potential, pressure, and temperature can be controlled independently. We have also shown with some examples that repulsive interactions can hold a nonadditive system in equilibrium under completely open conditions. While the implementation of the proposed scheme does not present any further difficulties other than those inherently associated with the isobaric-isothermal and grand canonical ensembles, we remark that nonadditivity is required to observe equilibrium states in the unconstrained ensemble. The proposed method paves the way for new developments in simulations of systems that exchange heat, work, and matter with their environment. Beyond long-range interacting systems and confined systems as considered here, such environmental conditions can be relevant, e.g., in small self-assembled aggregates [29][30][31] whose structures are stabilized by the competition of repulsive and attractive interactions.
S. R. thanks Peter Sollich for discussions and for suggesting Refs. [14,15]. A. C. acknowledges financial support from INFN (Istituto Nazionale di Fisica Nucleare) through the projects DYNSYSMATH and ENESMA.

REPLICA ENERGY AND THERMODYNAMIC RELATIONS
Here we explicitly integrate the energy balance equation for the ensemble, following the arguments in Ref. [22] (see also [16]). As stated in Eq. (1) of the main text, for a collection of N independent replicas with total energy, entropy, volume, and number of particles E t = N E, S t = N S, V t = N V , and N t = N N , respectively, one has the general expression Keeping all single-system properties constant, one ob- Now the integration is to be performed on the number of replicas only. Integrating between 0 and N yields and therefore we arrive at E = E −T S +P V −µN , which is always valid, in particular for any N .

PROBABILITY DENSITY IN THE UNCONSTRAINED ENSEMBLE
Here we present a derivation of the probability density in the unconstrained ensemble describing the system in a particular N -particle configuration, given in Eq. (3) of the main text. We follow Ref. [2], extending the standard arguments for other ensembles.
The system has N particles, volume V , and temperature T , and we recall that its canonical partition function takes the form where β = 1/k B T , Λ is the de Broglie thermal wavelength, and U(s N ; V ) is the potential energy, k B being the Boltzmann constant. In addition, reservoirs a and b are assumed to be ideal gases and their canonical partition functions are given by and respectively, where reservoir a has N a − N particles and volume V a and reservoir b has N b particles and volume We now consider the limit in which the reservoirs are infinite. Let us first focus on reservoir a. We take N a → ∞, V a → ∞ with N a /V a = ρ a , and use the limit N a /N → ∞ in the partition function such that In this limit, a change in the number of particles of the system does not change the chemical potential µ of reservoir a. Since the reservoir is an ideal gas, its chemical potential is given by µ = k B T ln(Λ 3 ρ a ) and hence In this limit, a change in the volume of the system does not change the pressure P of reservoir b. Since this reservoir is an ideal gas as well, the density ρ b = N b /V b can be written as ρ b = βP and hence, ( Using these formulas, the partition function Q b becomes Taking into account (5) and (7) to compute the partition function of the total system Q tot = Q a QQ b , in the limit N a , V b → ∞ the probability density P(N, where is the unconstrained partition function of the system. Here the factor βP is included to make Υ a dimensionless quantity, as it is usually done for the isothermal-isobaric partition function [2,13,16]. Using the canonical partition function (1), the probability density (8) can be rewritten as From this expression one directly gets the probability density P(N, V ; s N ) in a particular N -particle configuration, as given in Eq. (3) of the main text.

SIMULATION DETAILS
In this section we give more details about the simulations of the examples considered in the main text and the explicit expressions of the acceptance probabilities.
The number of Monte Carlo (MC) moves in a cycle is defined by m = N av + N ex + 1, where N av and N ex are fixed integers. To implement the algorithm, we generate a random integer R such that 1 ≤ R ≤ m and attempt a particle displacement if R ≤ N av , a volume change if R = N av + 1, and a particle exchange with the reservoir (insertion or removal with the same probability) otherwise. In this way, on average, per cycle the algorithm performs N av particle displacements, N ex particle exchanges, and one volume change. During all the simulations we set N ex = 1.
Before the production run, a calibration stage is performed in the simulations followed by a thermalization stage [32]. In the calibration stage, the maximum particle displacement and maximum volume variation are periodically updated to achieve an acceptance ratio of about 50%. Also, in this stage we periodically set N av = N to enforce that N av is close to the average number of particles which is a priori unknown. In the subsequent thermalization stage, simulations are carried out with all parameters fixed, and the average number of particlesN is computed. At the end of this stage, we set N av =N .
Finally, keeping all parameters fixed, we compute the averages in the production run, including the final value of the average number of particlesN . For each point in the plots of the simulations shown in the main text, the number of MC moves has been set to 10 6 during calibration, 10 9 during thermalization, and 3 × 10 9 or 8 × 10 9 in the production stage.
The acceptance probabilities for insertion and removal of particles and for volume changes in the unconstrained ensemble are obtained from expressions (5)-(7) of the main text. In the first example (system with spatially constant, repulsive interactions), using the reduced control parameters µ * , P * , and T * introduced in the main text, these equations become P acc (N → N + 1) = min 1, V * e (µ * −N )/T * (N + 1) , P acc (V → V ) = min 1, e −P * (V * −V * )/T * e −N ln(V * /V * ) , (13) where V * = V / 3 and V * = V / 3 , being an arbitrary length unit. For the second example (confined hard spheres), using the reduced control parameters µ * and P * defined for this case in the main text, when there is no overlap the acceptance probabilities take the form P acc (N → N + 1) = min 1, A * H * e µ * (N + 1) , P acc (N → N − 1) = min 1, N e −µ * A * H * , P acc (H → H ) = min 1, e −P * A * (H * −H * ) e −N ln(H * /H * ) , (16) where A * = A/σ 2 , H * = H/σ and H * = H /σ, σ being the particle diameter. We stress that in this case volume changes are obtained by performing variations of the gap thickness from H to H at constant area A.