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

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–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 ignor ensemble finite-size effects, i.e., reach the thermodynamic limit. In particular, for a system of total number of particles N₀ at temperature T in a volume V₀, it is possible to consider a subdomain of volume V with an average number of particles ⟨N⟩. The system can be considered in the grand canonical ensemble if V is of the order of 1% of the total volume V₀.⁴ This size constraint implies using huge simulation boxes which, in most cases, demand 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 (AT) sub-domain of volume V, defined within the simulation box, is in contact with a reservoir of coarse-grained (CG) particles.^5–9 A smooth interpolation between atomistic and coarse-grained 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–15 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 μVT ensemble requires beforehand knowledge of the chemical potential μ of the bath. This condition is analogous to the situation in the NVT 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.¹⁶ 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 (IG) reservoir. In the context of the adaptive resolution method, a semi-grand canonical method was proposed to illustrate this point.¹² The idea to interpret the adaptive resolution methodology in terms of a simulation with a particle reservoir was first formulated in Ref. 11. This identification of the coarse-grained system is completely general as illustrated in the case of atomistic/mesoscopic continuum adaptive resolution models^17,18 in which open boundary conditions can be readily enforced.^19–21 Nevertheless, the sampling resulting from particle insertion/deletion events, even if only applied in the coarse-grained region, becomes less representative as the density, concentration, and complexity of the system increase.

To consider these points, here we present a method that combines a particle insertion/deletion algorithm with the Hamiltonian adaptive resolution framework (H-AdResS).^22,23 Following the method suggested in Ref. 24, we replace the coarse-grained model by a reservoir of non-interacting thermalized 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.^25–27 

There is a permanent and ongoing interest in simulating non-equilibrium phenomena in inhomogeneous fluids.^28–42 With the open-boundary molecular dynamics method proposed here, it is possible to investigate such phenomena. 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 Fig. 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 conserves, on average, 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 is organized as follows: In Sec. 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 Sec. III, and finally, discuss, conclude, and outline research directions in Sec. IV.

In adaptive resolution simulations, it is possible to couple a target system with an ideal gas reservoir of particles at constant chemical potential.^24,43 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 open-system molecular dynamics simulations in equilibrium and, more importantly, nonequilibrium conditions. As an illustration, we study the prototypical example of liquid flow across a narrow channel. All simulations are performed with the LAMMPS simulation package,^44,45 where the method is implemented.

Before introducing the particle insertion algorithm, we validate H-AdResS to study a confined liquid between parallel walls.

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–7 atomistic (AT) and coarse grained (CG) representations of the system are concurrently present in the simulation box and coexist in thermodynamic equilibrium.^5–7 In particular, by using a position-dependent switching field λ(x) (see Fig. 1), it is possible to write a Hamiltonian (H-AdResS)^22,23 for the whole system. Since AT 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.^24,43 In this context, we write the Hamiltonian for a system composed of N = N_s + N_w solvent (s) and wall (w) molecules as⁴⁶ 

with (r, p) being the positions and momenta and $K=∑i=1Npi2/2mi$ being the total kinetic energy of the system. The term

is a harmonic, nearest-neighbour V^w(r > r_cut) = 0 potential used to restrain the position of the wall molecules. κ is the spring constant, r₀ 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 σ⁻³ with parameters κ = 1000 ϵσ⁻² and r₀ = 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. We impose this constraint to solve a technical problem, namely, due to the pressure imbalance present between the AT and IG regions, the wall might deform and its elastic response could affect the behavior of the solvent. Moreover, this constraint enforces no-slip conditions in the system. Indeed, thermal fluctuations of wall particles enhance the decoupling of wall–fluid interactions, thus increasing the slip length.³⁰ 

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 Fig. 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, is introduced in the Hamiltonian to compensate non-physical, drift forces due to gradients of the switching field. Moreover, it ensures momentum conservation on average and guarantees a uniform density throughout the simulation box. Note that the FEC is only applied to confined liquid particles instantaneously present in the hybrid region. 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. 45 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).²⁴ 

We model all the solvent–wall interactions with the truncated and shifted Lennard-Jones potential,

where the units of length, energy, and mass are defined by the parameters σ, ϵ, and m, respectively. In the following, we report the results in Lennard-Jones (LJ) units with time τ = σ(m/ϵ)^1/2, temperature ϵ/k_B, and pressure ϵ/σ³. 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 the 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.

In this example, the number of solvent and wall particles is fixed with N_s = 97 020 and N_$w$ = 48 000, 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 the Nose–Hoover thermostat and barostat for 5000 τ with time step size of δt = 10⁻³ τ. The temperature is fixed at k_BT = 2.0 ϵ with damping coefficient 10 τ and the pressure is fixed anisotropically at P = 2.65 ϵσ⁻³ 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 liquid particles to the total volume of the simulation box (⁠$ρeq*=NS/V$⁠), is reported in the inset of Fig. 2 for different fluid–wall interactions. This definition of density is useful, however, arbitrary. The volume accessible to the fluid cannot be uniquely determined because the definition of the contact surface between the fluid and the walls is somewhat unclear.³¹ 

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 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 Fig. 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 ϵ²⁴).

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 shown in Fig. 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.

It has been demonstrated that for a liquid in equilibrium at a given temperature T and composed of N molecules inside a container of volume V, the chemical potential equals the chemical potential of the uniform system.⁴⁷ However, as mentioned before, the definition of the volume available to the fluid is ambiguous.³¹ For such a reason, we avoid identifying the free energy compensation ΔH^s with the excess chemical potential of the uniform system. Hence, in Sec. II B, to precisely enforce a constant chemical potential in an open system, we introduce and validate the particle insertion algorithm for a bulk liquid system.

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 Sec. II A, as well as in Ref. 24. 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₀, and chemical potential μ, has exactly N particles follows the Poisson distribution⁴⁸ 

with N^* being the mean number of particles in the volume V₀. In the ideal gas case, N^* can be written in terms of the chemical potential of the system,

with β = 1/k_BT and λ being the mean thermal wavelength. In the limit N, N^* ≫ 1 with |(N − N^*)/N^*| ≪ 1, we obtain

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_BT with ρ^* = N^*/V₀. We are interested in fluctuations around a target density ρ^*; therefore, we rewrite P(N) in terms of ρ = N/V₀ as

where, in principle, k_μ = V₀/ρ^*, 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 that the present density ρ increases by ν is given by

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.

In conventional grand canonical simulations, the target system can interchange particles with an infinite ideal gas reservoir at constant chemical potential.⁴⁹ Alternatively, the atomistic/coarse-grained coupling present in adaptive resolution simulations provides a suitable playground to sample the grand canonical ensemble.^10–15 In our particular case, the size of the reservoir can be substantially increased since particles in the IG region are non-interacting.^24,43 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 ϵσ⁻³ and temperature k_BT = 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. (1) and the obtained FEC terms corresponding to the target density. The velocity and force of the ghost particles are set to zero, so they do not move during the H-AdResS parameterization runs. The on the fly calculation of compensation of the drift and thermodynamic forces are updated every 5 × 10⁴ time steps, during 10 × 10⁶ time steps (200 iterations). The resolution interval is divided into 200 bins of size Δλ = 0.005 and the length of the simulation box is uniformly discretized into slabs of size Δx = 0.5 σ. We employed values of c = 1 ϵ, σ_ρ = 2 σ, and ł = 4 σ for smoothing and scaling the thermodynamic force.

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 the other at higher, ρ > ρ^*, density.^50,51 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 Fig. 4, the evolution of density as a function of simulation time is presented in both cases, $ρL,RIG<ρ*$ (a) and $ρL,RIG>ρ*$ (b).

It is apparent from Fig. 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 behavior 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 behavior 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 Sec. II A. We remove the periodic boundary conditions along the x direction to decouple left and right reservoirs and introduce repulsive walls to confine the ideal gas particles. We use the Langevin thermostat to control the temperature in all regions at k_BT = 2 ϵ with damping coefficient 10 τ. The time step is δt = 10⁻³ τ, and for each case, the total simulation run is 14 × 10⁶ 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.

The combination of inducing a density gradient and enforcing a reference density generates an excess of frozen particles in the left reservoir. To keep the simulation running, we periodically relocate these frozen particles from left to right reservoirs. This technical aspect does not affect our results since frozen particles have zero momentum.

As both reservoirs are confined 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⁶ time steps each. To remove the effect of the transient behavior, we do not consider the results of the first block (see the inset of Fig. 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 shown in the inset). It is apparent that 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.

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 Fig. 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 the largest density gradient and it is 0.015 ϵ (less than 1% 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.⁵² A simple dimensional analysis justifies such a comparison. To this end, we employ the Knudsen number, defined as Kn = λ/D_h with $λ=1/2ρπσ2$ being the mean free path of the fluid particles and D_h = 20 σ being 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.⁵³ 

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 emphasize 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.

To compute the pressure, we discretize the simulation box in cubic cells of linear size 1 σ and build the pressure profile by considering a sequence of cells aligned along the x-axis and located at a distance 5 σ from the rigid surface. Hence, the structured layers of fluid close to the wall are not taken into account, which results in effectively counting the bulk contribution to the pressure. We compute the stress tensor per-atom S for particles instantaneously present in the cell. Finally, based on the Irvin–Kirkwood derivation,⁵⁴ we obtain the pressure tensor for every cell of volume V as P = −(S_xx + S_yy + S_zz)/3V⁵⁵ and average contributions along the y and z directions to obtain the pressure as a function of x. We emphasize here that in case one requires a more precise calculation of the pressure tensor, especially close to inhomogeneous regions, it is essential to use the method of planes introduced in Ref. 31.

Plot of pressure profiles [Fig. 7(a)] for different induced densities, in the case of purely repulsive interaction with the walls, indicates that the pressure gradient is indeed constant. As a matter of fact, we find that there is a linear relation between the pressure difference measured across the AT region $pRAT−pLAT$ and the nominal pressure difference $pRIG−pLIG$ for different values of the induced density, as indicated in Fig. 7(b). The slope of the line is 0.47. This value coincides with the constant pressure gradient across the resolution, namely, $LATLAT+2LHY+Lcorr=510$⁠. The quantity L_Corr = 20 σ measures the distance from the interface of the HY region to the left reservoir up to the point at which the pressure reaches the equilibrium pressure. The constant pressure gradient obtained in the AT region can be understood in terms of the virial expansion of the pressure in terms of the density for a LJ system, namely, relatively small gradients in density (of maximum 0.08 σ⁻⁴ here) induce ideal-gas-like response in the atomistic system.

To compare with the result expected from Eq. (12), we calculate the corresponding average of u_x(z) in the atomistic region and obtain the velocity profiles plotted in Fig. 8 for all liquid–wall interactions. The solid lines represent the parabolic functions fitted to the data points. In all cases and apart from relatively small statistic fluctuations, the simulated velocities follow closely the expected parabolic profile.

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., $umax∝(ρR*−ρL*)$⁠, as presented in Fig. 9.

Finally, we can use u_max to compute the viscosity η as

For the WCA case, we use all the values of ΔP/Δx (Fig. 7) and u_max (Fig. 9), corresponding to the density imbalances considered in our study. We estimate D_h = 19 σ from the density profiles in Fig. 3. Finally, we obtain η = 6.3 ± 0.6 ϵσ⁻³τ. This result is comparable with the value obtained using Müller–Plathe’s method⁵⁶ (η = 5.2 ± 0.3 ϵσ⁻³τ). The discrepancy might be due to lack of statistics and the uncertainty related to the ambiguous definition of the volume accessible to the fluid, counted in the value of D_h.

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₀ and volume V₀ is considered to be in the grand canonical ensemble provided that the condition V/V₀ ≈ 0.01 holds.⁴ For the Lennard-Jones systems investigated here, a subdomain of volume V ∼ (10 σ)³, with the estimate for the correlation length of the system as 10 σ, would require V₀ ∼ 100 V. With a density ρ = 0.647 σ⁻³, the simulation would need N₀ ∼ 10⁶. 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 emphasize 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.^{49–51,57–63} 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. 24, 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,^{29,53,64–66} current non-equilibrium molecular dynamics methods include external forces that might introduce artifacts in the simulation resulting from non-conservation of linear momentum.⁶⁷ Instead, the method proposed here conserves, on average, momentum. 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.

In summary, 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, volume, and chemical potential. In addition to the straightforward calculation of the chemical potential, the use of H-AdResS allows one to study liquid mixtures directly. In this context, we introduced a particle insertion/deletion algorithm that operates, at minimal computational expenses, on the ideal gas reservoir. Approaches exploiting similar ideas are available in the literature;^58,68 however, they lack the flexibility provided by the coupling to the ideal gas system. The proposed method allows one to perform constant chemical potential simulations under various conditions. More importantly, by studying the pressure-driven flow through a channel, we showed that it is also possible to perform well-controlled non-equilibrium molecular dynamics simulations.

R.P. and R.C.-H. contributed equally to this work.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors thank Pietro Ballone and Dominic Spiller for critical reading and insightful comments. M.H., K.K., R.P., and R.C.-H. gratefully acknowledge funding from SFB-TRR146 of the German Research Foundation (DFG). This work was supported by the European Research Council under the European Union’s Seventh Framework Program (No. FP7/2007-2013)/ERC (Grant Agreement No. 340906-MOLPROCOMP). This project received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant No. 758588).