Open-Boundary Hamiltonian adaptive resolution. From grand canonical to non-equilibrium molecular dynamics simulations

We propose an open-boundary molecular dynamics method in which an atomistic system is in contact with an infinite particle reservoir at constant temperature, volume and chemical potential. In practice, following the Hamiltonian adaptive resolution strategy, the system is partitioned into a domain of interest and a reservoir of non-interacting, ideal gas, particles. An external potential, applied only in the interfacial region, balances the excess chemical potential of the system. To ensure that the size of the reservoir is infinite, we introduce a particle insertion/deletion algorithm to control the density in the ideal gas region. We show that it is possible to study non-equilibrium phenomena with this open-boundary molecular dynamics method. To this aim, we consider a prototypical confined liquid under the influence of an external constant density gradient. The resulting pressure-driven flow across the atomistic system exhibits a velocity profile consistent with the corresponding solution of the Navier-Stokes equation. In contrast to available computational methods in which external forces drive the system far from equilibrium, this approach conserves momentum and closely resembles experimental conditions. The presented method can be used to study various direct and indirect out-of-equilibrium conditions in complex molecular systems.


I. INTRODUCTION
Computational and experimental communities routinely cooperate by comparing the results obtained from their respective methods. However, such comparisons are intrinsically limited in scope because real systems approach the thermodynamic limit, whereas model systems usually have a finite number of particles. Indeed, standard computer simulations frequently consider closed systems whose fixed number of particles introduces finite-size effects due to the de facto simultaneous use of different statistical ensembles [1][2][3].
With the computational power nowadays available, it is tempting to ask whether it is possible to increase the size of the system and safely ignored ensemble finite-size effects, i.e. reach the thermodynamic limit. In particular, for a system of total number of particles N 0 at temperature T in a volume V 0 , it is possible to consider a subdomain of volume V with an average number of particles N . The system is in the grand canonical ensemble if V is of the order of 1% of the total volume V 0 [4]. This size constraint implies using huge simulation boxes which, in most cases, demand a tremendous computational effort.
A simplification of the physical representation of the particles in the reservoir alleviates this computational load. This idea is the essence of the adaptive resolution method: an atomistic sub-domain of volume V , defined within the simulation box, is in contact with a reservoir of coarse-grained particles [5][6][7][8][9]. A smooth interpolation between atomistic and coarsegrained forces, acting on molecules present at the interface between the two regions, ensures a consistent description of the whole system. Indeed, the adaptive resolution framework is a robust method to perform simulations in the grand canonical ensemble [10][11][12]. To build upon these results, we discuss two components that one should consider to make use of the adaptive resolution method to perform well-controlled open-boundary molecular dynamics simulations.
First, a simulation in the µV T ensemble requires to know beforehand the chemical potential µ of the bath. This condition is analogous to the situation in the N V T ensemble, in which one first requires to fix the number of particles N . In the original adaptive resolution framework, it is possible to obtain the difference in chemical potential between atomistic and coarse-grained representations of the system [13]. If one knows the chemical potential of the coarse-grained model, then it is straightforward to obtain and control µ for the atomistic one.
It is thus convenient to use the simplest possible coarse-grained representation such that the calculation of the chemical potential does not involve an additional external computation. Second, a system is in the grand canonical ensemble if it is in thermal and chemical equilibrium with an infinite particle reservoir. Computer memory limitations prevent the possibility to consider an infinite number of particles. Instead, a particle insertion/deletion algorithm coupled to the simulation setup effectively ensures this requirement by allowing the interchange of particles with an infinite ideal gas reservoir. In the context of the adaptive resolution method, a semi-grand canonical method was proposed to illustrate this point [14]. The idea to interpret the AdResS methodology in terms of a simulation with a particle reservoir was first formulated in Ref. [15] and subsequently generalised to atomistic/mesoscopic continuum adaptive resolution models [16,17] in which open boundary conditions can be readily enforced [18][19][20]. Nevertheless, the sampling resulting from particle insertion/deletion events, even if only applied in the coarse-grained region, become less representative as the density, concentration and complexity of the system increases.
To consider these points, here we present a method that combines a particle insertion/deletion algorithm with the Hamiltonian adaptive resolution framework (H-AdResS) [21,22]. Following the method suggested in Ref. [23], we replace the coarse-grained model by a reservoir of non-interacting thermalised particles (ideal gas). In this case, the applied external potential used to ensure a uniform density through the simulation box balances the excess chemical potential of the atomistic model. Therefore, the atomistic system is at constant chemical potential with a reservoir of ideal gas particles. We introduce the particle insertion/deletion algorithm, operating on the ideal gas reservoir, to overcome the limitations existing on available methods due to high density/concentration and system complexity conditions. The method is thus capable of performing constant µ molecular dynamics simulations without the necessity to include external forces and to compensate for depletion of particles in the reservoir [24,25].
It is possible to study non-equilibrium phenomena with this open-boundary molecular dynamics method. In particular, we consider a confined liquid such that its ideal gas reservoir is under the influence of a constant density gradient. Initially, a uniform density profile is enforced parallel to the surfaces (see Figure 1). Upon equilibration, a density gradient imposed in the reservoir induces a pressure-driven flow in the system with a velocity profile consistent with the corresponding solution of the Navier-Stokes equation. The method closely resembles experimental conditions and avoids artefacts present in current methods due to the non-conservation of linear momentum. Furthermore, the intrinsic arbitrariness of the ideal gas reservoir opens the possibility to study various direct and indirect non-equilibrium conditions.
The paper reads as follows: In section II, we validate the use of H-AdResS to study confined liquids and introduce the particle insertion/deletion algorithm. We present the results for pressure-driven flow in section III and finally discuss, conclude and outline research directions in section IV.

II. MODEL
In adaptive resolution simulations it is possible to couple a target system with an ideal gas reservoir of particles at constant chemical potential [23,26]. Particularly, it is possible to set the chemical potential of the target system by controlling the number density in the ideal gas reservoir. In this work, we implement a particle insertion algorithm that operates on the ideal gas region and permits fluctuations in the number density around a target value. Consequently, the standard H-AdResS setup now allows one to perform opensystem molecular dynamics simulations in equilibrium and, more important, nonequilibrium conditions. As an illustration, we study the prototypical example of liquid flow across a narrow channel.
Before introducing the particle insertion algorithm, we validate H-AdResS to study a confined liquid between parallel walls.

A. H-AdResS for a confined liquid
Let us consider a single component liquid confined between parallel walls with normal perpendicular to the x-axis and separated by a distance D. In the adaptive resolution method (AdResS) [5][6][7], atomistic (AT) and coarse grained (CG) representations of the system are concurrently present in the simulation box and coexist in thermodynamic equilibrium [5][6][7].
In particular, by using a position-dependent switching field λ(x) (see Figure 1) it is possible to write a Hamiltonian (H-AdResS) [21,22] for the whole system. Since AA and CG representations of the system follow different equations of state, an additional term is included in the Hamiltonian to guarantee a constant density (implying constant chemical potential)  across the simulation box. Finally, the choice of the coarse grained representation is rather arbitrary and can even be reduced to a reservoir of noninteracting -ideal gas (IG) -particles [23,26]. In this context we write the Hamiltonian for a system composed of N = N s + N w solvent (s) and wall (w) molecules as [27]: with (r, p) positions and momenta and K = N i=1 p 2 i /2m i the total kinetic energy of the system. The term is a harmonic, nearest-neighbour V w (r > r cut ) = 0 potential used to restraint the position of the wall molecules. κ is the spring constant, r 0 is the equilibrium length and r j,j is the distance between a given pair of wall particles. The wall particles are located on a fcc lattice of density 0.9 σ −3 with parameters κ = 1000 σ −2 and r 0 = 1.1626 σ. V w (r) is only applied during the initial equilibration. For production runs, the surface particles are frozen in their final equilibrium positions.
The last term in the Hamiltonian (1) is written as: The solvent molecules i interact via an intermolecular potential modulated by a switching field λ i ≡ λ(x i ) (see Figure 1). The switching field takes values 0 in the IG region where solvent particles only interact repulsively with wall particles and 1 in the AT region where solvent particles fully interact with both solvent and wall particles. A smooth interpolation between 0 and 1 is defined in the hybrid (HY) region. The free energy compensation (FEC) term, ∆H s , compensates non-physical forces due to gradients of the switching field and guarantees a uniform density throughout the simulation box. Note that the FEC is only calculated for the confined liquid particles. The starting point of the simulations corresponds to a fully atomistic case. Once these simulations are equilibrated, the adaptive resolution setup is turned on and ∆H s is computed on-the-fly following the procedure described in Ref. [28] included in the equilibration of the H-AdResS run. Finally, for homogeneous systems, ∆H s (λ = 1) corresponds to the system's excess chemical potential µ exc . This procedure has been identified with a thermodynamic integration in space, i. e. spatially-resolved thermodynamic integration (SPARTIAN) [23].
We model all the solvent-wall interactions with the truncated and shifted Lennard-Jones where the units of length, energy and mass are defined by the parameters σ, and m, respectively. In the following, we report the results in LJ units with time σ( /m) 1/2 , temperature /k B and pressure /σ 3 . For a given cutoff radius r cut the value U α,β (r cut ) is evaluated and subtracted from Eq. (3). The parameters used to describe all interactions between species (α, β) for different regions within the simulation box are presented in Table I with fixed solvent-solvent (s,s) and ideal gas-wall (ig,w) values. For solvent-wall (s,w) interactions in the AT region we consider the purely-repulsive Lennard-Jones potential (WCA) as well as truncated and shifted Lennard-Jones potential with r cut = 2.5σ with varying interaction strength modulated by the parameter η with η = 0.5, 1.0, 1.5, 2.0 and 2.5. Additionally to the purely repulsive case (not included in this table), the latter interactions are modulated in the AT region by the parameter η, whereas in the IG region is defined as purely repulsive.
In this example, the number of solvent and wall particles is fixed with N s = 97020 and N w = 48000, respectively. The size of the box is set by L x = 164.41 σ and L y = 49.32 σ while L z is fixed by the system's pressure with variations in the range of L z = 24.75 σ, · · · , 25.61 σ.
For the case of homogeneous liquid, i.e. no confining walls, L z = 18.74 σ. The initial fully atomistic equilibration is carried out in the NPT ensemble using Nose-Hoover thermostat and barostat for 5000 τ with time step size of δt = 10 −3 τ . The temperature is fixed at k B T = 2.0 with damping coefficient 10 τ and the pressure is fixed anisotropically at P = 2.65 σ −3 with damping coefficient 100 τ by applying a force normal to the walls. The final equilibrium density, which is defined as the ratio of number of solvent particle to the total volume of the simulation box (ρ * eq = N S /V ), is reported in the inset of Figure 2 for different fluid-wall interaction.
To validate the reliability of H-AdResS to study confined liquids, we verify that the FEC terms, ∆H s , converge. The size of AT and HY regions is 50 σ and 15 σ and the H-AdResS parameters are listed in Table II. As it was mentioned before, the application of the FEC on the system leads to a constant density profile across the simulation box. In this case, the system reaches an equilibrium state in which atomistic and ideal gas particles have equal chemical potential. The evolution of the FEC as a function of time is plotted in Figure 2.
After 200 iterations, which corresponds to 2000 τ simulation time, (see Table II for more detail), the algorithm converges to µ = 2.33 ± 0.008 for the bulk liquid (no confinement), in a good agreement with the previously reported value for LJ fluid (µ = 2.33 ± 0.01 [23]). Somewhat expected, the obtained values of ∆H s are inversely proportional to the To investigate the structure of the liquid inside the AT region, we calculate the equilibrium particle distribution function normal to the wall, g(z), and compare it with the results obtained from the fully atomistic simulation. The distribution function is obtained by binning the atomistic region along the z-axis and time averaging over all bins with equal z-coordinate. As it is shown in Figure 3, the surface-induced layering structure is well reproduced by the H-AdResS simulation. Moreover, the comparison of the bulk particle distribution (g(z); for |z| < 5σ) shows that the density of the atomistic region in the adaptive setup has converged to the atomistic reference density. These results confirm the suitability of our method to study confined liquids under the conditions here specified.
At this point, it is important to mention that the identification of ∆H s with the chemical potential of the solvent is not straightforward in the case of a confined liquid. To compute ∆H s we use as a reference density the nearly constant value obtained in the central region between the parallel plates, whereas a full identification of ∆H s with the chemical potential should explicitly include the dependence with the distance normal to the surface. Hence, in the next section, to unambiguously enforce a constant chemical potential in an open system we introduce and validate the particle insertion algorithm for a bulk liquid system.

B. Particle insertion algorithm
The proposed grand canonical molecular dynamics method consist of two parts: first, the AT/IG constant chemical potential coupling that has been already discussed in the previous section as well as in Ref. [23]. Second, we allow particle exchange between the IG region and an ideal gas reservoir used to control the chemical potential of the system. The details of the particle insertion algorithm applied on the IG region are the subject of this section.
We start by assuming that the IG region is in the grand canonical ensemble. The probability that the IG region, at temperature T , volume V 0 and chemical potential µ, has exactly N particles follow the Poisson distribution [29] with N * the mean number of particles in the volume V 0 . In the ideal gas case, N * can be written in terms of the chemical potential of the system with β = 1/k B T and λ the mean thermal wavelength. In the limit N, N * 1 with |(N − This is a normal distribution with mean value N * and standard deviation √ N * . In physical terms, this corresponds to the well known result for the isothermal compressibility κ of the ideal gas, i. e. κ = 1/ρk B T with ρ * = N * /V 0 . We are interested in fluctuations around a target density ρ * , therefore, we rewrite P (N ) in terms of ρ = N/V 0 as where in principle k µ = V 0 /ρ * but in the following we treat it as a free parameter related to the width of the distribution of possible values for the target density ρ * . With this probability distribution, we introduce the Metropolis algorithm used for particle insertion.
The probability to accept a move, namely, that the present density ρ increases by ν, is given by acc(ρ → ρ + ν) = min[1, exp(−k µ ν(ν + 2(ρ − ρ * )))] , and correspondingly Concerning the fluctuations around the target density, we select values of ν according to the normal distribution Finally, the particle insertion Monte Carlo moves are performed every ∆T exch time steps during which the number of particles is averaged.
FIG. 4. Evolution of the density in the left (ρ IG L ), right (ρ IG R ) ideal gas regions and atomistic (ρ AT ) regions when the system has initially lower (a) or higher (b) density than the equilibrium reference density (ρ * eq ).
In conventional grand canonical simulations, the target system can interchange particles with an infinite ideal gas reservoir at constant chemical potential [30]. Alternatively, the atomistic/coarse-grained coupling present in adaptive resolution simulations provides a suitable playground to sample the grand canonical ensemble [10][11][12]14]. In our particular case, the size of the reservoir can be substantially increased since particles in the IG region are non-interacting [23,26]. Alternatively, by introducing density fluctuations in the IG region, we also ensure interchange of particles with an infinite reservoir of ideal gas particles at constant chemical potential.
To verify grand canonical conditions, we run an equilibrium simulation for a bulk LJ liquid, i. e. with the confining walls removed, at a given pressure P = 2.65 σ −3 and temperature k B T = 2 . These conditions define the target density ρ * eq = 0.647 σ −3 and the corresponding chemical potential. Upon obtaining the equilibrated all-atom configuration, "ghost" particles are placed in each reservoir and then H-AdResS is performed using the Hamiltonian of Eq.  [28,31], where the method is implemented.
To verify that the particle insertion protocol drives the system to a target density, ρ * , we start with two versions of the system at the same temperature, but one at lower, ρ < ρ * , and one at higher, ρ > ρ * , density [32,33]. In both cases, we apply the FEC obtained from the target system to set the target chemical potential, and run the open boundary simulation using ρ * as an input for the particle insertion protocol. In Figure 4 the evolution of density as a function of simulation time is presented in both cases, ρ IG L,R < ρ * (a) and ρ IG L,R > ρ * (b). It is apparent from Figure 4 that the density in the three regions, left, right and AT converges to the reference density in all cases, independently of the choice of insertion frequency k µ . In general this result verifies that the open simulation setup described corresponds to a constant chemical potential molecular dynamics simulation. The behaviour at short times indicates that the FEC works to restore the density in the atomistic region (dashed lines) by depleting (a) / increasing (b) the number of particles in the reservoir. The effect of the particle insertion algorithm is thus to bring the density of the reservoir (solid and dotted lines) to the target value and equate the chemical potential across different simulation regions.
Finally, the behaviour of the system as a function of k µ is consistent with the interpretation given in terms of the width of the distribution of ρ * .
In the final section, we return to the confined liquid problem and use the particle insertion algorithm in a nonequilibrium molecular dynamics setup to induce a density gradient throughout the system.
In this section, we start with the equilibrated configurations for the confined liquid considered in the previous section. Here, we use the Langevin thermostat to control the temperature in all regions at k B T = 2 with damping coefficient 10 τ . The time step is δt = 10 −3 τ and for each case, the total simulation run is 14 × 10 6 time steps. To induce the density gradient between the atomistic region and ideal gas reservoirs, the particle exchange algorithm is applied independently in each reservoir. The reference number densities of particles in the left reservoir is set as the equilibrium number density ρ * L = ρ * eq . The reference number density in the right reservoir is increased with respect to the equilibrium number density.
We use k µ = 0.1 and ∆T exch = 100 δt for all simulation sets. The reservoir (frozen) particles are also exchanged between the two reservoirs every 100 time steps to balance the number of particles in the reservoirs. As both reservoirs are separated by repulsive walls, the number density of each reservoir is calculated over a smaller control volume than V IG which does not contain the depletion layer close to the right and left repulsive walls. Additionally, the width of the depletion layers changes when the solvent changes from LJ fluid (fully atomistic simulation) to ideal gas particles (H-AdResS) which changes the total available volume. This difference in volume can be compensated by increasing the hard core size of the repulsive LJ potential between the walls and ideal gas particle such that the number density far from the depletion zone equals the one in fully atomistic simulations. In our study, however, such difference is negligible.
Concerning statistics, all values are reported using the block averaging method. For each case, we divide the total simulation run into seven uniform blocks of 2 × 10 6 time steps each.
To remove the effect of the transient behavior, we do not consider the results of the first block (See inset of Figure 5). The average values of each block are calculated and then we report the average and standard deviation values of all blocks. Figure 5 shows the number density profile of particles when the system reaches the steady state (as it is shown in the inset). It is apparent the induced density gradient in the atomistic region and the ripples observed there are generated by the interaction with the surface. As expected, the difference in densities between the right and left reservoirs is equal to the actual nominal difference. The density profile at the IG L /HY interface exhibits bumps that become more distinct upon increasing the density gradient. These can be attributed to a mismatch in mobility due to particles changing identity from atomistic to non-interacting.
Thus, particles accumulate before entering the left reservoir, and once there the particle insertion algorithm flattens the profile to reach the reference density. For clarity, the profiles have been shifted by a constant quantity.
To effectively isolate the effect of the induced density gradient in the atomistic region, it is necessary to guarantee that the temperature is constant throughout the simulation box.
This is presented in Figure 6, clearly indicating a uniform temperature of the entire system.
As a matter of fact, the maximum deviation of temperature from the mean value occurs in the case of induction of largest density gradient and it is 0.015 (less than 1 percent of the mean temperature). This is found to be in the middle of the hybrid region where the resolution of the ideal gas changes effectively into LJ fluid.
At this point, we are in the position to compare our simulation results with hydrodynamics, namely, the Poiseuille flow equation [34]. A simple dimensional analysis justifies such a comparison. To this end, we emply the Knudsen number, defined as Kn = λ/D h with λ = 1/ √ 2ρπσ 2 the mean free path of the fluid particles and D h = 20σ the height of the channel.
Knudsen numbers for the system under consideration vary between Kn = 0.014, · · · , 0.016 indicating that a parallel can be drawn between continuum fluid dynamics and our simulation results [35].
A plane Poiseuille flow is created in a fluid with density ρ and dynamic viscosity µ confined between infinitely long parallel plates, separated by a distance D h , when a constant pressure gradient is applied along the axis of the channel x. The resulting flow is unidirectional. At low Reynolds number, the Navier-Stokes equation can be written as where the axial velocity u x (z) is only a function of the z-coordinate and the applied pressure gradient ∆p/∆x is constant. Using the boundary conditions u x (z = ±D h /2) = 0 (no-slip condition), we obtain Before comparing this result with the velocity profile obtained from simulations we emphasise two aspects. First, to obtain the result in Eq. (12) one assumes that the fluid is incompressible, and the LJ system under consideration is not. However, the Navier-Stokes equation for a compressible fluid contains a term proportional to ∇ · u that in our case is equal to zero because the flow is unidirectional along the x-direction and the magnitude of the flow velocity does not change along this axis. Second, the result in Eq. (12) was obtained by assuming that ∆p/∆x is constant. In our model we need to verify that the constant pressure gradient induced in between the right and left reservoirs creates a constant pressure gradient across the AT region. Plot of pressure profiles (Figure 7 (a)) for different induced densities, in the case of purely repulsive interaction with the walls, indicates that this is indeed the case.
As a matter of fact, we find that there is a linear relation between the pressure difference measured across the AT region, p AT R − p AT L and the nominal pressure difference p IG R − p IG L for The resulting velocity profile implies that the maximum velocity occurs at the center of the channel −u max = u x (z = 0). This maximum velocity is linearly proportional to the nominal density difference, i. e. u max ∝ (ρ * R − ρ * L ), as presented in Figure 9.

IV. DISCUSSION AND CONCLUSIONS
In computer simulations, the investigation of a large closed system allows one to sample the grand canonical ensemble. A subdomain of volume V and an average number of particles N inside a system with a fixed number of particles N 0 and volume V 0 is considered to be in the grand canonical ensemble provided the condition V /V 0 ≈ 0.01 holds [4]. For the Lennard-Jones systems investigated here, a subdomain of volume V ∼ (10 σ) 3 , with the estimate for the correlation length of the system as 10 σ, would require V 0 ∼ 100V . With a density ρ = 0.647 σ −3 , the simulation would need N 0 ∼ 10 6 . This number is one order of magnitude larger than the number of particles considered here, which already highlights the advantage of using the proposed method. Additionally, we emphasise that a detailed assessment of the relative atomistic, hybrid and ideal gas region sizes might suggest that we could further reduce the size of the system without affecting our main results.
Instead of the direct approach discussed above, it is perhaps more convenient to use one of the available grand canonical molecular dynamics methods [30,32,33,[36][37][38][39][40][41][42]. A common ingredient in most approaches consists of inserting particles with a given system-dependent probability. In general, when the system under consideration is simple in terms of the force field used for its description, or if the system is at low density/concentration conditions, this is the method of choice. Far away from such conditions, the particle insertion protocol becomes highly inefficient. In this respect, the adaptive resolution framework constitutes an alternative for existing methods. In particular, for the Hamiltonian adaptive resolution discussed in Ref. [23], the target system is in contact with a reservoir of ideal gas particles at constant temperature, volume and, by ensuring a uniform density across the simulation box, also at constant chemical potential. The combination with a particle insertion algorithm operating in the ideal gas region guarantees a reservoir of infinite size, thus completing the definition of grand canonical ensemble. Therefore, the method proposed here is a robust strategy to perform open-boundary molecular dynamics simulations, mainly when the system under consideration is dense, or highly concentrated in the case of mixtures, and regardless of the complexity introduced by the force field description.
A straightforward change in the geometry and periodic boundary conditions in the simulation box allows one to decouple the ideal gas reservoir. Hence, it is possible to simultaneously impose different temperature, density and concentration conditions on the system. In the particular case of induced density gradients [35,[43][44][45][46], current non-equilibrium molecular dynamics methods introduce external forces that might introduce artefacts in the simulation resulting from non-conservation of linear momentum. Instead, the method proposed here conserves momentum on average. Moreover, the simplicity of the reservoir gives the possibility to study different out-of-equilibrium conditions for complex molecular systems, which constitutes a significant improvement over state-of-the-art simulation methods.
To conclude, in this work, we presented a method to perform open-boundary molecular dynamics simulations. We used the H-AdResS framework where the atomistic target system is in physical contact with a reservoir of non-interacting particles at constant temperature,