Hybrid molecular-continuum simulation techniques afford a number of advantages for problems in the rapidly burgeoning area of nanoscale engineering and technology, though they are typically quite complex to implement and limited to single-component fluid systems. We describe an approach for modeling multicomponent hydrodynamic problems spanning multiple length scales when using particle-based descriptions for both the finely resolved (e.g., molecular dynamics) and coarse-grained (e.g., continuum) subregions within an overall simulation domain. This technique is based on the multiscale methodology previously developed for mesoscale binary fluids [N. D. Petsev, L. G. Leal, and M. S. Shell, J. Chem. Phys. **144**, 084115 (2016)], simulated using a particle-based continuum method known as smoothed dissipative particle dynamics. An important application of this approach is the ability to perform coupled molecular dynamics (MD) and continuum modeling of molecularly miscible binary mixtures. In order to validate this technique, we investigate multicomponent hybrid MD-continuum simulations at equilibrium, as well as non-equilibrium cases featuring concentration gradients.

## I. INTRODUCTION

Over the last two decades, numerous approaches for concurrently modeling phenomena across multiple length scales have been proposed, driven by a broad spectrum of emergent technologies at interfaces and the nanoscale. One tactic is to employ a molecular description for spatial regions that necessitate high detail (e.g., a surface where a nanoscale physical or chemical process occurs) and a simpler continuum approximation for other bulk regions where high resolution is not required.^{1–9} The first approach of this type was proposed by O’Connell and Thompson,^{1} who coupled molecular dynamics (MD) simulations to a continuum fluid, where the Navier-Stokes solution was obtained using finite-element methods, with exchange of fluxes between the two regions. Hadjiconstantinou and Patera^{3} provided another formulation through the Schwarz alternating method, making it possible to couple MD and continuum domains by iteratively matching their boundary conditions. This modeling strategy was subsequently applied to the moving contact line problem, resolving singularities that appear in the continuum solution by using a molecular description for the region near the boundary.^{2} More recently, we developed a different MD-continuum approach featuring a MD region coupled to a continuum region simulated via a fluctuating hydrodynamic solver called smoothed dissipative particle dynamics (SDPD).^{9} Unlike previous hybrid simulation methodologies, we used this stochastic *particle*-based hydrodynamic modeling technique for the continuum part of the simulation,^{10,11} allowing for a consistent and intuitive particle-based description throughout the entire system.^{12} Our MD-continuum multiscale strategy is based on the adaptive resolution approach “AdResS” originally developed by Praprotnik *et al.*,^{13–19} which interpolates between spatial regions of molecular and coarse-grained fluid descriptions, gradually switching additional degrees of freedom on and off as molecules traverse the intermediate zone between the two regions. Prior to our application of AdResS to the MD-continuum particle problem, it was applied to numerous systems, from liquids comprised of tetrahedral molecules^{13,14,18} to molecular-continuum simulations of water,^{17} and grand canonical simulations.^{20,21} Zavadlav and Praprotnik recently used AdResS to couple MD to dissipative particle dynamics (DPD)^{22} using the SWINGER algorithm for particle clustering.^{23}

However, a major shortcoming of all of these so-called “multiscale” hybrid simulation strategies is that they are limited to single-component scenarios, even though systems involving solvated species are ubiquitous in modern molecular and interfacial physics. A significant challenge is faithfully reproducing the dynamics of the dissolved species in solute-solvent systems across multiple spatiotemporal scales, including atomistic ones, which is essential to numerous fundamental problems and applications. For example, a number of time and length scales appear in drug delivery and particle transport in cardiovascular flow,^{24,25} where surface chemistry and functionalization is important over small scales while hydrodynamics dominates over larger ones. Other applications include electrolyte solutions at interfaces and corrosion,^{26,27} in which oxidation occurs at the molecular level and transport and the microstructure of the material are relevant over micrometer lengths and beyond. Hybrid multiscale techniques may also be useful in coarse-graining solvent in bulk regions when simulating biomolecules in solutions containing salt, including protein interactions in biomolecular assembly^{28} and ion channels in cellular membranes.^{29–31} Another such example is modeling surface nanobubbles and nanodroplets,^{32–36} where dissolved gas transport over potentially large length scales may occur in bulk solvents away from the bubble or drop. In light of the potential applications, the goals of this paper are (1) to outline an intuitive and generalizable framework for capturing the transport of solute in a miscible mixture due to diffusion and advection across multiple length scales and (2) to assess the applicability of the method in simple and appropriately devised test cases.

A major difficulty in multiscale multicomponent simulations is that at the sub-nanometer level, matter is discrete (e.g., a molecule is either a solute or a solvent), whereas over larger length scales and in the continuum approximation, the fluid is described in terms of fields that can assume any value between 0 and 1 (i.e., a fluid volume has a concentration associated with it). Therefore, in traditional MD-continuum hybrid methods, it is necessary to perform spatial or temporal averages that can translate properties in the “discrete” molecular region into continuum variables. Conversely, for free exchange of solute, mass, and momentum between the two regions, conveying the averaged information contained in the continuum fields to a collection of particles with discrete identities is also required. The major challenge is the communication between the molecular and continuum regions since they treat concentration in distinct ways, and this issue is illustrated in Fig. 1. Performing coupled MD-continuum simulations therefore requires being able to bridge these two fundamentally different concentration descriptions.

There has been prior work towards interfacing molecular regions containing a fluid mixture to continuum ones. A hybrid simulation strategy for bridging MD regions to finite-difference continuum solutions that allows for heat and solute transport was provided by Alexiadis *et al.*^{37,38} In that work, the authors couple MD simulations of an ideal fluid mixture to finite-element solutions of the transport equations by solving a Poisson-type equation for the continuum variables, where information from the atomistic region enters through the non-homogeneous term. Techniques for interfacing MD regions to numerical solutions of Navier-Stokes have also been extended to allow for modeling electroosmotic transport.^{39} For the most part, however, hybrid multiscale multicomponent simulation schemes are largely absent in the literature, in spite of an abundance of important systems.

In this paper, we provide a framework for coupling discrete particle concentration descriptions (MD or coarse-grained molecular) to a particle-based continuum description, opening the possibility for multiscale modeling of ideal mixtures, including simulations featuring atomically resolved regions. In particular, our ultimate goal is an algorithm for interfacing discrete particle models such as MD, DPD, and other high resolution descriptions applicable to finer scales, to a more coarse SDPD fluid representing the continuum region. As we show below, two key features are required to make a smooth transition from continuum to molecular. First, we use a particle-based continuum description in the form of SDPD, which naturally interfaces with the inherently particle-based molecular world. Second, we introduce a “discrete SDPD” region to bridge the transition from this continuum region in which the volume elements have a continuous range of solute concentrations to molecularly-sized SDPD fluid volume elements that are either pure solvent or solute. The design and implementation of the transition between these two regions is the key focus of the present work. The latter discrete SDPD region then offers a natural and convenient transition from SDPD to MD.

SDPD is a Lagrangian particle-based approach for solving the fluctuating hydrodynamic equations of Landau and Lifshitz^{40–42} and has been applied to a number of problems at the nano- and mesoscales.^{11,25,43–47} An attractive feature of SDPD is that the characteristic length scale (i.e., the degree of coarse-graining in the fluid) is a tunable parameter, controlled by the smoothing length *h*, which determines the size of the fluid volumes (or “particles”) into which the fluid is discretized. Thus, a fine-scale SDPD continuum region with large fluctuations can be coupled to increasingly coarse SDPD descriptions (i.e., regions with more massive SDPD particles having larger *h*, smaller fluctuations, and at a lower number density), giving significant computational savings due to the reduction in the number of particles.^{9,48,49} SDPD is a thermodynamically consistent fluctuating formulation of smoothed particle hydrodynamics (SPH), which is a numerical approach for solving the Navier-Stokes equations by decomposing the fluid into a collection of Lagrangian volumes, or interpolation nodes, whose positions evolve in time according to an equation of motion that approximates the hydrodynamic equations.^{50–52} SDPD is obtained by using GENERIC^{53–55} to incorporate thermal noise in the evolution equations for the SPH particles through the fluctuation-dissipation theorem. While the techniques described in this paper are in principle applicable to both SPH and SDPD as the continuum description, we focus on the latter.

The work reported here was preceded by three previous studies from our group. In the first,^{48} we developed an algorithm to couple SDPD descriptions of a single-component fluid with different levels of resolution (i.e., different smoothing lengths corresponding to different size particles) in adjacent regions of space, thus allowing for multiscale modeling via SDPD. In the second,^{9} we constructed a strategy that couples MD to SDPD for simulating single-component fluid problems, allowing for multiscale modeling that includes a MD-based atomistic description of the fluid at the finest level of resolution. Finally, in the third study from our group, we derived a fluctuating continuum multiscale generalization of the SDPD method for two-component fluid mixtures.^{49} The resulting equations of motion are similar to the ones obtained by Ellero *et al.* for a fluid composed of solvated Hookean dumbbells.^{56} In our multicomponent SDPD simulation strategy, each particle (whose interpretation is a small volume of fluid) carries some amount of solute and solvent and therefore has a concentration associated with it. Hence, this gives a particle-based approach to solving the fluctuating hydrodynamic and solute transport equations for a miscible two-component solution. Here, we address the last step required for full multiscale simulation of fluid mixtures by introducing a MD-based description at the finest scale, where MD particles are either a solute or solvent, coupled to a SDPD fluid with each particle consisting of both solute and solvent and having a specified concentration.

We offer a brief review of these algorithms in Sec. II, but to maintain reasonable brevity, we do not discuss these methods in detail here, instead referring the interested reader to Ref. 9 for the MD-continuum method, and Ref. 49 for the multiscale multicomponent extension of SDPD. Here, we focus on a strategy for coupling molecular and continuum simulations of a homogeneous two-component fluid mixture, hence reconciling the single-component MD-continuum coupling techniques in Ref. 9 to the multicomponent generalization of SDPD in Ref. 49. To illustrate the resulting procedures, we perform simulations that probe correct reproduction of equilibrium properties, as well as non-equilibrium tests of quasi-1D diffusion. Although the ideas presented here are limited to ideal solutions, they are still valid for numerous problems where the ideal mixing approximation is appropriate due to dilute concentrations, i.e., any situation in which the solute dynamics is governed by the classic advection-diffusion equation. Moreover, the technique outlined in Secs. III–VII may be extensible to non-ideal solutions and non-isothermal conditions as well.

In the first part of this paper, we focus on the algorithm for coupling a discrete particle description such as MD (i.e., a fluid where each particle is identified as being either a solute or solvent) to a continuum fluid description having a concentration value associated with each particle (see Fig. 1). Although we ultimately want a continuous region modeled using SDPD and the discrete particles to be MD particles, we begin by considering a simpler problem: a multiscale simulation featuring a “fine” SDPD fluid where each particle assumes a discrete identity of either Φ_{i} = 1 (solute) or Φ_{i} = 0 (solvent), coupled to a more coarse-grained SDPD fluid where each fluid volume has a solute concentration that can have any value between 0 and 1. As discussed in Sec. VI, this kind of coupling between two SDPD fluids with different concentration descriptions is a crucial component of the transition from MD to SDPD continuum description. As discussed in Sec. VI, the latter will consist of a transition from MD to discrete SDPD, followed by the transition from discrete SDPD to the continuum SDPD description. The coupling of multiple scales in SDPD is achieved by allowing the larger, coarse-grained particles to split (or “refine”) into the finer discrete solute and solvent particles, and for these fine particles to combine (or “coarsen”) into the coarse ones. The obvious question that we address in the following is how to split continuous particles into smaller ones having discrete identities, and vice versa.

Rules have already been developed for splitting and combining single-component SDPD particles in a multiscale simulation that preserve mass, momentum, and mean particle position.^{48} Previously, we also described splitting and combining rules for a fluid mixture modeled using both fine and coarse SDPD regions,^{49} though in this case we used a continuous concentration description throughout, and this is not sufficient for coupling SDPD to a finely resolved molecular region where particles have discrete identities. Therefore, we only need rules for the splitting and combining steps that transition between the solute/solvent particle description and the continuum SDPD description with a concentration for each particle.

We begin by providing a brief overview of multiscale SDPD and the coupled MD-SDPD approach for single-component systems in Sec. II. The MD-SDPD techniques require modification for simulating fluid mixtures, and in Sec. III we consider the refining step, where a large, continuous SDPD particle enters the refining region and must split into smaller particles with a discrete solute or solvent identity. In Sec. IV, we then consider the reverse move and discuss the combining step, where two (or more) fine particles with discrete identities, which we identify using a concentration value of either Φ_{i} = 1 or Φ_{i} = 0, depending on whether they are a solute or solvent, respectively, are combined into a single coarse particle having a concentration Φ_{j}. In Secs. V and VI, we test these techniques by performing simple multiscale simulations (both continuum-continuum and MD-continuum), demonstrating that proper equilibrium solute distribution is recovered. Finally, in Sec. VII, we extend this multiscale multicomponent approach to non-equilibrium situations featuring concentration gradients and consider the simple problem of quasi-1D diffusion as a case study.

## II. MULTISCALE MODELING USING SMOOTHED DISSIPATIVE PARTICLE DYNAMICS

Smoothed dissipative particle dynamics (SDPD) is a particle-based approach for solving the hydrodynamic equations in the Lagrangian frame. For the case of an isothermal incompressible Newtonian fluid and ideal mixing, the appropriate transport equations governing the momentum and solute dynamics are^{57}

and

Here, *ρ* is the mass density of the fluid (combined mass density of the solute and solvent for the case of a binary mixture), Φ is the concentration defined as the ratio of solute mass to total mass, *η* is the fluid (or fluid mixture) viscosity, and *D*_{F} denotes the diffusivity of the solute with units ML^{−1} t^{−1}. In SDPD, the fluid domain is decomposed into a collection of fluid volumes or “particles” whose positions, concentrations, and velocities are evolved using forms of Eqs. (1) and (2) obtained from a discretization using an interpolant function. The resulting equations of motion for the SDPD particles are^{10,45}

for the particle momenta, while the concentration of each particle evolves according to^{58}

In the SDPD governing Eqs. (3) and (4), **v**_{i} denotes the velocity, Φ_{i} denotes the concentration, *m*_{i} denotes the mass, *ρ*_{i} denotes the density, and *p*_{i} denotes the pressure of the *i*th particle. Additionally, we have the relative position, velocity, and concentration of particles *i* and *j*, given by $eij\u2261rij/rij$ with **r**_{ij} ≡ **r**_{i} − **r**_{j}, **v**_{ij} ≡ **v**_{i} − **v**_{j}, and Φ_{ij} ≡ Φ_{i} − Φ_{j}. We denote the separation between particles *i* and *j* as *r*_{ij} ≡ |**r**_{i} − **r**_{j}|. *W*_{ij}(**r**_{i} − **r**_{j}, *h*) ≡ *W*_{ij} is the smoothing kernel. The level of coarse-graining in SDPD is an input parameter controlled by the smoothing length *h*, which determines the extent of the function *W*_{ij} and controls the size of the SDPD particles.

Small values of *h* give mesoscopic fluid volumes that stochastically exchange momenta and solute with each other due to thermal fluctuations. Therefore, in addition to the reversible and dissipative interactions between SDPD particles in Eqs. (3) and (4), we have the random contributions to the momentum and solute equations of motion,^{10,45,49}

where *d***W**_{ij} and *dV*_{ij} are tensorial and scalar increments of the stochastic Wiener process, respectively, and $dW^ij$ is the traceless and symmetric part of *d***W**_{ij}. *A*_{ij} and *G*_{ij} are noise amplitudes chosen such that the fluctuations introduced into the particle momenta and concentrations by Eq. (5) balance the dissipative terms in Eqs. (3) and (4) in accordance with the fluctuation-dissipation theorem. The noise amplitudes are given by^{10,45,49}

where *k*_{B} is Boltzmann’s constant, *T* is the system temperature, and *m*_{0} is the mass of a single solute or solvent molecule or atom, which for simplicity we assume all have the same mass. Hence, *m*/*m*_{0} gives the number of atoms or molecules contained within a single SDPD fluid particle. We have also introduced Θ_{i} = Φ_{i}(1 − Φ_{i}).

The smoothing length *h* (i.e., SDPD resolution) can vary spatially in order to selectively coarse-grain parts of the simulation domain by allowing large SDPD particles to split into small ones and small particles to combine into more massive ones (see Refs. 48 and 49). The rules for splitting particles are constructed such that particle momenta, mass, and center of mass are conserved,^{49} e.g., for the case of a single particle *i* splitting into two particles *j* and *k* each having half the mass of the parent, we have

where the random displacement vector Δ**r**_{ij} = **r**_{i} − **r**_{j} is generated according to Δ**r**_{ij} = *random*(*h*_{i}). Note that particles can split into integer values greater than two, though we describe this case for simplicity. Similarly, for the reverse move where particles *j* and *k* combine into a single particle *i*, the rules are

These rules are used to interface two SDPD fluids featuring a different resolution using a transition zone that is subdivided into three parts: (1) splitting, (2) overlap, and (3) combining subregions^{49} (see Fig. 3 in Sec. IV). In this approach, the coarse region is discretized into large fluid volumes (particles), whereas the finely resolved one is composed of a collection of less massive particles at a higher number density. If a large particle is transported to the refining region, it splits into fine SDPD particles. Similarly, a fine particle that crosses into the combining region is merged with another nearby fine particle into a single coarse one. The fine and coarse particles coexist in the overlap region, which must be sufficiently large to ensure a smooth transition of fluid properties across the interface between the two different scales;^{59} if this region is too small, particles near the interface separating the two resolutions will interact with different numbers of neighbors from particles in the bulk, resulting in error. As the finely resolved continuum domain is coupled to increasingly coarse SDPD regions with more massive particles and ever larger values of *h*, the fluctuations introduced by Eq. (5) become less significant.

While a purely continuum method, the equations of motion governing the SDPD particles [Eqs. (3) and (5)] are similar to those in molecular dynamics, and hence we can perform coupled MD-continuum simulations by allowing particles to change the type from MD to SDPD and vice versa depending on their location inside the simulation box using an adaptive resolution approach.^{13–19} In previous work,^{9} we coupled MD to SDPD for the single-component case using a switching function *s*(**r**) that is unity in the MD region, zero in the SDPD region, and smoothly and monotonically transitions between these values inside a buffer zone between the atomistic and continuum domains. Then, the reversible pair interaction between particles *i* and *j* is written as a sum between MD and SDPD conservative forces with weighting given by *λ* = *s*(**r**_{i})*s*(**r**_{j}),^{9}

Here, $FijSDPDrev$ denotes the pressure force in the SDPD equation of motion [i.e., the first term on the right-hand side of Eq. (3)]. Equation (9) guarantees that particles interact through MD forces in the atomistic region, SDPD interactions in the continuum one, and through a linear combination of both force types inside the intermediate (i.e., “buffer”) zone between the MD and SDPD fluids. The SDPD dissipative forces [second term on the right-hand side of Eq. (3)] and random forces [Eq. (5)] are applied to all particles inside the buffer region without weighting in order to thermostat the fluid.

Note that there is a work associated with transforming particles from MD to SDPD type due to the different fluid descriptions, which is supplied by a so-called “thermodynamic force” that acts on particles in the buffer to counter chemical potential gradients that emerge due to multiple resolutions.^{9,19} By incorporating a thermodynamic force and using the MD-SDPD force interpolation approach outlined above, we previously showed that a MD region can be coupled to an atomically scaled SDPD fluid (i.e., one where continuum particles have the same mass as the MD ones, ensuring mass conservation). This finely resolved SDPD domain can in turn be coarse-grained by coupling to SDPD regions having more massive particles and larger values for the smoothing length *h* by allowing the fine particles to combine into more coarse ones and coarse particles to split into fine ones. This approach has been applied to single-component systems as detailed in Ref. 9 and is extended to binary fluid mixtures in Secs. III–VII.

## III. REFINING CONTINUOUS PARTICLES INTO PARTICLES WITH DISCRETE IDENTITIES

We now turn to the issue of refining a continuous particle with solute concentration Φ_{j} into “discrete” particles that are either a solute or solvent (i.e., that have a zero or one value for Φ_{j}). Note that in prior work,^{49} we used continuum SDPD particles for the entire problem domain, with particles of different sizes (different *h*) in separate parts of the fluid domain. In this case, rules for combining or splitting particles were straightforwardly determined from conservation laws such that the total mass, momentum, and solute remain constant. In the present work, however, we consider the transition from a region in which the particles have a discrete identity as either a solute or solvent, to a region where the particles are fluid volumes containing both solute and solvent with a solute concentration between 0 and 1. The rules in Eqs. (7) and (8) for the particles’ masses and momenta upon splitting or combining are unchanged and are still used here. In this work, we focus on the concentration variable alone and how concentration values are assigned upon splitting a continuum particle or combining discrete particles that may either be a solute or solvent. This is important for situations such as ones with a continuum fluid coupled to a MD domain since in the molecular region particles have discrete identities. It is clear that when a continuum particle splits into two or more discrete particles, some of the dissolved species may be lost or gained and will not be precisely conserved globally. Therefore, in the following, we do not preserve the total solute number exactly but instead enforce a constant chemical potential within the transition region where particles divide and combine. In other words, we assume that the fluid within the transition region is locally in equilibrium with a bath and hence the concentrations of particles crossing this region are determined by the effective chemical potential of this bath. Our approach is related to grand canonical sampling that has a long and well-studied history^{60} and has already been applied in the context of AdResS.^{20,21}

For simplicity, we first consider the case where coarse particles are twice as massive as the fine ones, i.e., we have fine particles (both solute and solvent) with mass and smoothing length *m*_{1} and *h*_{1}, and large particles with *m*_{2} = 2*m*_{1} and *h*_{2} > *h*_{1}. Thus, when a coarse particle enters the refining region, it splits into two particles, each having half the mass of the parent. Later, we generalize this approach for situations involving large fluid volumes that are *n* times more massive than the fine ones, where *n* is an integer greater than one. For the present case, when a large particle *i* enters the refining region, it divides into two discrete daughter particles *j* and *k*. When assigning concentrations to the daughter particles, we have the following possibilities (Fig. 2): (1) The large particle *i* splits into two solute particles (Φ_{j} = Φ_{k} = 1), (2) the large particle splits into one solute and one solvent particle (e.g., Φ_{j} = 1 and Φ_{k} = 0), or (3) the large particle splits into two solvents (Φ_{j} = Φ_{k} = 0). Note that for case (2), there are actually two possibilities: either Φ_{j} = 1 and Φ_{k} = 0 or Φ_{j} = 0 and Φ_{k} = 1. Therefore, this splitting move has a degeneracy of two (Ω_{2} = 2), whereas there is no degeneracy in the other two cases (Ω_{1} = Ω_{3} = 1).

Next, we construct rules for selecting from these three possibilities such that a constant chemical potential is enforced. Consider a mixture of two types of particles within a box (i.e., one where the *j*th particle has a concentration of either Φ_{j} = 1 or Φ_{j} = 0) with average bulk solute concentration $\Phi $ where $0\u2264\Phi \u22641$, and suppose we remove a particle from the box. The probability of drawing a solute particle (Φ_{j} = 1) is simply proportional to the average solute concentration, $P\Phi j=1\u221d\Phi $. Therefore, the probability for drawing two solute particles (i.e., particles with concentration Φ_{j} = Φ_{k} = 1) is given by $P(\Phi j=1)P(\Phi k=1)\u221d\Phi 2$. Note that this is only exactly true if the first particle is replaced before drawing the second. However, for large systems comprised of many SDPD particles, this is a good approximation in the case without replacement as well (a simple calculation shows that the percent change in the average concentration due to removing a single particle from a box containing *N*_{p} SDPD particles scales as ∼1/*N*_{p}). Moreover, when the chemical potential is fixed, this becomes exact since the average concentration is maintained constant by contact with a bath. Hence, the probability for case (1) is

where we have included the degeneracy factor Ω_{1}. Similarly, the probability for drawing a solvent molecule is $P(\Phi j=0)\u221d1\u2212\Phi $. Therefore, the probability for drawing one solvent and one solute is *P*(Φ_{j} = 1)*P*(Φ_{k} = 0), and thus for case (2), we have

Finally, the probability for a large particle splitting into two solvent particles [case (3)] is

Equations (10)–(12) with degeneracies Ω_{1} = Ω_{3} = 1 and Ω_{2} = 2 give the probabilities for each of the three cases in Fig. 2. When a particle crosses into the refining region and divides into two daughter particles, we choose from one of the three possibilities using the above probabilities. Note that they are already normalized due to the use of mass fractions, $\u2211lPl=1$. In practice, selecting from one of these three outcomes is done by generating a uniformly distributed random number *X* [*X* ∼ *U*(0,1)]. Then, if *X* < *P*_{1}, the particle splits into two solutes, *P*_{1} ≤ *X* < *P*_{1} + *P*_{2} results in a solute and solvent, and *X* ≥ *P*_{1} + *P*_{2} gives two solvent particles.

It is straightforward to extend these rules for large fluid volumes dividing into more than two small particles. For the general case of a coarse particle splitting into *n* fine particles, the probabilities become

Here, the *l*th degeneracy factor is the $nl\u22121$ binomial coefficient. By coupling to a fluid with coarse particles significantly more massive than the fine ones, there is greater reduction in the total number of particles used, giving further computational savings. However, note that if the discrepancy in the particle masses is too large, this leads to artifacts in properties such as the density distribution, as well instability due to heating upon particle insertion.^{49,59} In addition, compressible systems or ones comprised of more complex molecules may experience issues when combining or refining particles, which can be remedied by bundling the finer particles together when they are near the transition zone using strategies such as SWINGER.^{22,23} The SWINGER algorithm could also be used to remove the need for a splitting/combining region that is spatially distinct from the MD-to-continuum transition zone, in which case the rules described in this paper can be used to assign concentrations to particles that enter the transition zone. In this work, we present results from the *n* = 2 (Sec. V) and *n* = 6 (Secs. VI and VII) cases.

## IV. COARSENING DISCRETE PARTICLES INTO PARTICLES WITH CONTINUOUS CONCENTRATION LABELS

We now consider the converse situation when two small, discrete particles are combined into a single continuous particle. This combining step occurs inside the coarsening region (Fig. 3) and is the reverse move to the particle splitting illustrated in Fig. 2. Note that for the particle splitting discussed in Sec. III, the amount of the solute is not conserved exactly, and instead we enforce the desired chemical potential (i.e., we assume that the transition region is locally in equilibrium with a constant chemical potential bath). Similarly, we construct a coarsening procedure that imposes the same target chemical potential (i.e., average concentration) as the splitting rules. Note that while it is possible to combine particles in a way that conserves solute, this is not possible when splitting particles, and it is important to instead impose a constant chemical potential condition to give symmetry between the splitting and combining moves. This guarantees that at equilibrium on average there is no loss or gain of solute due to the splitting and combining steps, even if the total amount of the solute is not conserved precisely.

The concentration inside a small, finite volume within a binary mixture can assume a number of different values due to thermal fluctuations since this region is free to exchange the solute with its surroundings. At equilibrium, this results in a distribution of concentration values. For a collection of SDPD particles, each having dimensionless mass *m*/*m*_{0}, these concentration values are described by a Gaussian distribution with variance

where Δ*μ* denotes the difference between the solute and solvent chemical potentials, a consequence of the constant mass constraint for the SDPD particles. Equation (14) is obtained following a derivation similar to the one for concentration fluctuations in Ref. 61. Assuming an ideal solution, i.e., $\Delta \mu =kBTln\Phi \u2212kBTln1\u2212\Phi $, the variance in Eq. (14) becomes simply $\Delta \Phi 2=m0\Phi 1\u2212\Phi /m$.^{49} In this case, concentration distribution is

Here, *m* is the mass of the coarse SDPD particle and *m*_{0} is the mass of an individual solute or solvent molecule; hence, the dimensionless mass *m*/*m*_{0} is also equivalent to the total number of molecules contained in a single fluid particle. We can therefore assign a concentration value to the newly formed daughter particle upon coarsening by sampling this distribution. Similar to the splitting move, the amount of the solute contained in the parents is irrelevant to the concentration of the newly formed fluid volume, which only depends on the chemical potential of the bath. Note that by drawing values at random from Eq. (15), the new continuous particle is automatically equilibrated with the surrounding fluid. Furthermore, this is readily extensible to coarsening steps where the large particle is arbitrarily more massive than the fine particles since the parameter *n* is proportional to the coarse particle’s mass, *m*, which affects the distribution in Eq. (15). Equation (14) gives a basis for determining the distribution one needs to sample in a coarsening step when dealing with non-ideal solutions where the chemical potential features a different functional form.

Alternative approaches for combining particles include giving the resulting daughter particle a concentration that is equal to the equilibrium concentration ($\Phi i=\Phi $) or assigning a value based on the concentrations of the parent particles such that the solute is conserved (e.g., a pair of discrete particles with concentrations Φ_{j} = 1.0 and Φ_{k} = 0.0 combine into a single coarse particle with concentration Φ_{i} = 0.5). However, we find that these approaches bias the solute distribution inside the coarsening region such that it no longer represents the appropriate equilibrium one, given by Eq. (15). When coarse particles that do not comply with the equilibrium distribution are then destroyed to form smaller fluid volumes, this may result in a net loss or gain of solute over time. In other words, when considering the distributions of solute left over after a collection of coarsening steps, and after a collection of refining steps, the distributions are perfectly symmetric when using Eq. (15) and the solute is conserved on average. On the other hand, if a different approach is used for the coarsening step, the solute left over from the combining steps may not balance the leftover solute from the refining steps and there will be a net gain or loss of solute over time (depending on the average concentration). Therefore, it is necessary to couple both the refining and combining regions to the same chemical potential bath, and in this work, we assign coarse particle concentrations by sampling Eq. (15).

## V. TEST OF DISCRETE-CONTINUOUS SDPD CONCENTRATION COUPLING AT EQUILIBRIUM

We use simple two-scale SDPD simulations to test this approach at equilibrium (see Fig. 3, Sec. IV). Specifically, we perform a multiscale simulation having a finely resolved SDPD fluid coupled to a coarse one where the particles are twice as massive, and the number density is halved. Concentrations are defined as dimensionless mass fractions, and all values are reported with the reduced unit convention in Refs. 48 and 49. The particles in the fine region are “discrete” and can only take on values of Φ_{j} = 1 (solute) or Φ_{j} = 0 (solvent), whereas the particles in the coarse region have a concentration assigned to them that can have any value between 0 and 1. The simulation of a discrete SDPD region coupled to a more coarse continuous SDPD region is a toy model for testing the discrete-coupling splitting and combining rules in Secs. III and IV, and an essential step towards coupled MD-continuum simulations of a fluid mixture, since coupling to MD requires an atomically scaled discrete SDPD fluid. Therefore, in order to test this discrete-continuum coupling approach, we construct a simple test system with a normal SDPD region coupled to one where we do not solve Eq. (4), hence restricting particles there to only have concentrations of either 1 or 0. Coupling to this artificial discrete SDPD construction is necessary for hybrid MD-continuum simulations, as shown in Sec. VI. The splitting probabilities for the transition zone between the continuum SDPD region and the discrete one are given by Eqs. (10)–(12), and the concentration of a coarse particle resulting from a combining step is drawn from Eq. (15).

Initially, 1000 fine particles are placed on a 10 × 10 × 10 cubic lattice inside a volume with dimensions 25 × 25 × 25. The coarse region is prepared by creating 512 particles on an 8 × 8 × 8 lattice inside a 25 × 25 × 25 region. This coarse fluid is placed adjacent to the fine one such that the normal vector to the interface between fine and coarse regions points in the *z*-direction. The system features periodic boundary conditions in all directions, and the two transition regions are located between *z* = 0 and 6, and *z* = 25 and 31, with each transition zone being divided into refining, overlap, and coarsening subregions having a width of 2. Hence, the total simulation box has dimensions 25 × 25 × 50, with the fine-grained region located at 6 < *z* < 25, and the coarse-grained one at 31 < *z* < 50. Before starting the equilibration, the particle lattice positions are translated in the positive *z*-direction by half the transition region width so that initially the interfaces between fine and coarse particle domains are precisely centered inside the overlap subregions.

For the SDPD model, we choose parameters that correspond to a Lennard-Jones fluid at reduced temperature *T* = 1.0 and mass density *ρ* = 0.8. With this mass density and lattice number density, the resulting particles have mass *m*_{1} = 12.5 in the fine region and *m*_{2} = 25.0 in the coarse (henceforth the 1 subscript denotes the fine region, and 2 is the coarse-grained one). With our reduced unit convention, these masses imply fluid particles where each one is comprised of a cluster of ∼13 Lennard-Jones atoms in the fine region and 25 atoms in the coarse one. The corresponding smoothing lengths are *h*_{1} = 3.00 and *h*_{2} = 3.75, respectively. At the selected state point, the fluid bulk and shear viscosities are *ζ* = 0.9^{62} and *η* = 1.9.^{63} Finally, in SDPD, it is necessary to specify the equation of state for the fluid, and thus its thermodynamic properties, in order to calculate the pressure force between particle pairs. Rather than using the full Lennard-Jones equation of state, we use a simple linearized version with the form $pi=cs2\rho i$,^{48,51,52} where *p*_{i} is the pressure at particle *i*, *ρ*_{i} is the density at the *i*th particle, and *c*_{s} is the speed of sound. Incompressible fluids in SDPD are approximated as quasi-incompressible by choosing a large value for the speed of sound such that particles are highly repulsive and density variations are low,^{52,58,64} typically within 3% of the target density.^{64} Hence, we select *c*_{s} = 5.0,^{48} which gives small fluctuations in density without requiring a prohibitively small time step. The particle concentrations are updated using the Euler-Maruyama integrator^{65} with a time step of Δ*t* = 0.01.

An important detail is that the self-diffusivity of the fine, discrete particles (MD or SDPD) must equal the solute diffusion constant in the continuum part of the simulation, i.e., the diffusive transport of solute should be uniform throughout the system independent of the degree of coarse-graining or solute description. While velocity and concentration fluctuations are statistically independent at equilibrium,^{61} velocity fluctuations can enhance solute transport in a non-equilibrium setting.^{66,67} We presently only consider an equilibrium case but still carefully match the diffusivities in anticipation of Sec. VII, where we impose a concentration gradient across the system. In the fine particle region, solute diffusion occurs purely due to the Brownian motion of the discrete particles, and hence it is described by the discrete particles’ self-diffusion constant *D*_{1}. However, the solute diffusivity in the coarse, continuous SDPD region has two distinct contributions since particles not only move but also exchange the solute with each other in a pairwise fashion. Thus, solute transfer occurs in the coarse SDPD region even if the particles are stationary due to solute fluxes between particle pairs [described by Eqs. (4) and (5)], and the actual thermal motion of the SDPD particles may enhance solute transport further. Matching diffusivities requires some care since for small fluid volumes in non-equilibrium conditions, velocity fluctuations lead to an enhancement of diffusive transport.^{66,67} Specifically, the diffusion constant characterizing the total or “effective” rate of solute transfer in the continuous SDPD region is written as a combination of these two contributions: (1) the “bare” Fickian diffusion constant *D*_{F} for the solute transferred from one particle to a neighbor (either due to a difference in the chemical potential between the two particles or due to a stochastic flux of solute) and (2) an additional quantity *D*_{2} that accounts for the random motion of the fluid particles due to velocity fluctuations. *D*_{2} is simply the self-diffusivity of the continuous SDPD particles due to thermal motion. The effective diffusion constant for the continuum SDPD fluid is therefore the sum of these two distinct contributions^{68}

From mean-square-displacement calculations, we find that the self-diffusion constant of the SDPD fluid particles with mass *m*_{2} = 25.0 is *D*_{2} = 0.0162. *D*_{1} = 0.0325 for SDPD particles with *m*_{1} = 12.5, which is the self-diffusivity of a discrete SDPD particle in the fine region [note that in the discrete SDPD domain, we do not solve the diffusion equation Eq. (4), and hence *D*_{F} = 0 and solute transport is purely due to the thermal motion of the SDPD particles].

Because the fine and coarse SDPD particles in a multiscale simulation feature a distinct self-diffusivity, the solute dynamics must be modified to equalize the rate at which the solute diffuses throughout the simulation box. To maintain a uniform diffusivity across the simulation box, the finely resolved fluid self-diffusion constant in the discrete particle region *D*_{1} therefore must equal the effective diffusion constant in the continuum SDPD region, *D*_{eff} [Eq. (16)], i.e., *D*_{1} = *D*_{eff}, which is achieved by tuning *D*_{F} for the continuous SDPD particles. Equation (16) suggests a choice of *D*_{F} = *D*_{1} − *D*_{2} = 0.0162 for the constant that enters the discretized diffusion equation for the coarse-grained particles in order to achieve this matching. Note that the diffusion constant in the discrete part of the simulation box *D*_{1} must be greater than the self-diffusivity of the SDPD particles in the continuum region (*D*_{1} > *D*_{2}); if this condition is not met, having uniform diffusivity would require a negative value for the bare diffusion constant *D*_{F}. If a simulation features a hierarchy of continuum SDPD domains, each having a different level of coarse-graining (i.e., smoothing length),^{49} then the Fickian diffusion constant *D*_{F} may need to be adjusted for each region such that the effective diffusivity is uniform [although for large particles, *D*_{2} is vanishingly small compared to *D*_{F}, and thus Eq. (16) becomes simply *D*_{eff} ≈ *D*_{F}. Hence, modifying the value of *D*_{F} for more coarse fluids is not necessary].

Equilibration of the fluid mixture features two stages. This is necessary since initially there is a slight transfer of mass from the fine region into the coarse one, a consequence of the higher self-diffusivity of the smaller SDPD particles relative to the coarse ones, *D*_{1} > *D*_{2}.^{48} After the distribution of particles across the system is equilibrated, and the number of particles in the fine and coarse regions stabilizes and fluctuates around an equilibrium value, the concentrations are reset to the average and thermally equilibrated. In other words, we equilibrate the total mass distribution first, and the solute distribution second. The combined preparation process is over 1 × 10^{6} time steps, and the production run is over 1.6 × 10^{7} time steps. Three separate cases, each having a different average concentration, are studied: (i) $\Phi =0.4$, (ii) $\Phi =0.5$, and (iii) $\Phi =0.6$.

The equilibrium concentration profiles from the production runs are shown in Fig. 4. The profiles are flat, without significant deviations due to the interface between the discrete and continuous regions, with errors in the three profiles of 1.37%, 0.93%, and 0.48% in the discrete region for cases (i), (ii), and (iii), respectively. In the continuous region, the concentrations for each of the three cases are within 0.23%, 0.19%, and 0.15% of the target value. By construction, the continuous SDPD region reproduces the correct fluctuations,^{49} i.e., the distribution of concentrations is given by Eq. (15). Throughout the simulation for case (i), there are a total of 123 023 time steps where refining events occur and 123 027 instances of coarsening ones, with particle splitting or combining occurring on average every 130 time steps. The rate for particle splitting is identical to the rate for combining since the simulations are at equilibrium. The average solute lost to the bath in case (i) due to a refining step over the course of the production run is 0.0269, while 0.0268 is the solute gained for particles combining; hence, there is approximately no net gain or loss of the solute within the interface region. Similarly, the average solute lost and gained per splitting and combining move is 0.0484 and 0.0499, respectively, for case (ii), and 0.0446 and 0.0470 for case (iii). In all tests, the net solute gained is balanced by the solute lost due to the reverse step.

## VI. HYBRID MD-CONTINUUM SIMULATIONS OF AN IDEAL FLUID MIXTURE AT EQUILIBRIUM

We now demonstrate how this technique enables MD-continuum coupling for fluid mixtures. In particular, we consider a binary fluid modeled across multiple resolutions, with a molecular region for the small scales, and a continuum description for the large ones. The hybrid molecular-continuum approach for single-component systems described in Ref. 9 is an essential element in this test. The method involves a MD region coupled to a fluid composed of molecularly sized SDPD particles, using an adaptive resolution approach^{13–19} that gradually interpolates between MD and SDPD interactions across a buffer region using Eq. (9). Due to the particle description employed across all length scales, momentum and mass conservation are intuitively prescribed.

Note that to mediate the MD-to-SDPD transition, this technique requires SDPD fluid “volumes” that have the same size as the individual molecules/atoms, which is unphysical.^{9} Despite being atomically sized, these particles are Markovian due to their coupling to a momentum bath, similar to MD fluids with a dissipative particle dynamics thermostat.^{69} This kind of highly resolved SDPD fluid still satisfies conservation laws and fluctuation-dissipation and is used to bridge MD to coarser continuum domains. Importantly, our initial tests demonstrate that it is possible to shrink the atomically resolved SDPD region described in Ref. 9 to zero size, minimizing the number of such ultrafine SDPD particles required for a multiscale simulation (results not shown). Thus, the only parts of the simulation box where atomically sized SDPD particles appear are in the refining and overlap regions in the interface between the MD and SDPD fluids (see Fig. 5).

Since the single-species MD-continuum approach briefly summarized in Sec. II has been discussed at length,^{9} we focus only on how the concentration is treated across the simulation box. In generalizing this MD-continuum approach to multicomponent systems, there is no ambiguity in how the molecular region is handled; we simply have a MD region featuring more than one species. Similarly, the continuum region is unambiguously modeled using the multicomponent SDPD formulation described in Ref. 49, where each fluid particle carries solute, has a concentration value assigned to it, and these concentrations are evolved according to the discretized diffusion equation. The MD-continuum coupling is handled as follows: using Eq. (9), a solute MD atom/molecule becomes a discrete SDPD particle with concentration Φ_{i} = 1 upon crossing the MD-to-continuum transition zone (orange area, Fig. 5) into the SDPD region, and vice versa. A solvent MD atom/molecule becomes a SDPD particle with Φ_{i} = 0. Then, these fine, discrete SDPD particles are coupled to more coarse, continuous ones across the fine-to-coarse SDPD transition zone (blue area, Fig. 5) using the discrete-continuous coupling described in Secs. III–V of this paper. Note that there is no pure atomically scaled SDPD fluid region as in Ref. 9; by removing this domain, the MD-to-continuum and fine-to-coarse SDPD transition regions are placed immediately next to each other. Hence, the overall interface is divided into these two parts: one where MD particles transition to continuum ones and a second part where these fine SDPD particles are coarsened (see Fig. 5).

A test for this coupling methodology is illustrated in Fig. 6(a) and is set up as follows: a MD wall is created by initializing 845 atoms on a 13 × 13 × 5 lattice inside a volume with dimensions 13 × 13 × 5. Immediately next to this region, we place the MD fluid domain by creating a 12 × 12 × 12 lattice of atoms inside a 13 × 13 × 13 volume. This region is located adjacent to the MD wall such that the *z*-direction is normal to the liquid-solid interface. Next, the transition zone is created by placing 1056 particles on an 11 × 12 × 8 lattice inside a region with side lengths 13 × 13 × 8. Once again, this part of the simulation box is positioned next to the MD one such that the interface between them has a normal in the *z*-direction. Finally, we initialize the coarse SDPD fluid by generating a region with dimensions 13 × 13 × 26 and placing a 6 × 7 × 14 lattice of particles inside. This fluid domain is next to the transition region, where the normal vector to the interface between them is in the *z*-direction. The system features an explicit MD boundary between *z* = 0 and *z* = 5 and a MD fluid between *z* = 5 and *z* = 18. The region between *z* = 18 and 23 is the MD-to-continuum transition zone, where MD and SDPD interactions are turned on and off using the adaptive resolution approach from Ref. 9, Eq. (9). Between *z* = 23 and *z* = 29 is the fine-to-coarse particle transition described in Secs. III and IV, as well as Ref. 49, where fine SDPD particles are transformed into coarse ones, and vice versa. This part is itself divided into a refining region (located 23 < *z* < 25), overlap region (25 < *z* < 27), and coarsening region (27 < *z* < 29). Particles situated between *z* = 29 and *z* = 46 are coarse SDPD fluid particles. The SDPD particles between *z* = 46 and *z* = 52 constitute the wall, and their positions and concentrations are not evolved in time. For more details about boundaries in SDPD, see Ref. 64.

We investigate the same thermodynamic state point (*T* = 1.0 and *ρ* = 0.8) and hence use identical values for the shear and bulk viscosities as in Sec. V. In order to approximately match both the compressibility and absolute pressure in the MD region, we use a linear equation of state with the form $pi=cs2\rho i\u2212\rho 0+p0$ (here *ρ*_{0} and *p*_{0} are the target density and pressure, respectively, where for the present case, *p*_{0} = 1.4).^{9} This is a less computationally intensive linear approximation to the full Lennard-Jones equation of state.^{70} Matching the absolute pressure of the SDPD fluid to the one in the MD region is required since a pressure mismatch will result in a transfer of mass across the interface separating the domains with different resolutions. For this test, we couple the MD fluid to a SDPD region with particles six times more massive; hence, splitting probabilities are given by Eq. (13) with *n* = 6, and the concentration of a newly formed coarse particle is drawn from a numerically pre-computed probability distribution of a single-scale SDPD fluid with *m*/*m*_{0} = 6 (discussed in more detail below). Given the initial lattice and choice of mass density, the coarse SDPD particles have a mass of 6 and smoothing length *h* = 2.35, while the fine particles (MD and discrete SDPD) have a mass of unity and *h* = 1.30. Since we have an ideal mixture, the SDPD model corresponds to a reference MD fluid where atom interaction energies and sizes are uniform, and the distinction between solute and solvent atoms is purely in terms of labels.

Recall that for very small fluid particles, matching diffusivities requires additional care since velocity fluctuations result in enhanced diffusion in non-equilibrium settings (see Sec. V). In other words, the self-diffusivity of the solute MD molecules should match the effective diffusion constant in the continuum part of the simulation, which is composed of a Fickian and self-diffusion contribution, as described by Eq. (16). Using mean-square-displacement calculations, the self-diffusion constant of the SDPD fluid volumes with *m* = 6.0 is calculated to give *D*_{SDPD} = 0.0486. Our MSD calculations are in agreement with the results of Kulkarni *et al.*,^{48} who found that SDPD particles with mass between *m* = 2.0 and 6.0 approximately match the self-diffusivity of their Lennard-Jones atom counterparts at the current thermodynamic state point. *D*_{LJ} = 0.0737 in the MD region, which is the self-diffusivity of a Lennard-Jones particle at the present conditions in a fluid with interaction energy and length *ε* = 1 and *σ* = 1, respectively, and potential cutoff *r*_{c} = 2.6. Uniform diffusivity across the simulation box requires that the MD diffusion constant *D*_{LJ} equals the effective diffusion constant in the SDPD region, *D*_{eff} = *D*_{F} + *D*_{SDPD}, i.e., *D*_{LJ} = *D*_{eff}, where *D*_{F} is the diffusion constant that enters the discretized SDPD form of the diffusion equation. Equation (16) therefore suggests a choice of *D*_{F} = 0.0252 for the Fickian diffusion constant that ensures this matching.

As in Sec. V, three separate equilibrium simulations are performed, each having a different average concentration ($\Phi =0.4$, $\Phi =0.5$, and $\Phi =0.6$). For the continuum part of the system, a Dirichlet boundary condition is applied at the solid boundary such that the concentration at this location equals the average, equilibrium value. In addition, a no-slip condition is applied to the velocity at the continuum solid-liquid interface. The system is equilibrated for 1 × 10^{6} time steps, and data is collected for 1.3 × 10^{7} time steps. We choose a time step value of Δ*t* = 0.001. Care must be taken to use a sufficiently small time step due to precision issues associated with very small particles that otherwise result in SDPD particles having a negative concentration or concentration greater than one.

An issue in the MD-SDPD coupling is that the continuum SDPD particles in the region adjacent to the MD region are very small, whereas the distribution for the coarsening step [Eq. (15)] assumes large particles ($m/m0>>1$) and is only Gaussian in the large-system thermodynamic limit. When particles are small such that concentration fluctuations are significant, the probability distribution is no longer Gaussian and it is not possible to directly use Eq. (15). Therefore we approximate the true ideal mixing probability distribution for small systems using a Gaussian. Hence, the distribution for the solute $P\Phi $ is numerically precomputed beforehand from a single-scale SDPD simulation, and a truncated Gaussian fit of this distribution is used in place of the analytical result given by Eq. (15). This approximation is a construction necessary for mediating the MD-continuum coupling.

Concentration profiles from the three equilibrium simulations are shown in Fig. 6(b). The errors in the profiles for the molecular, continuum, and transitional regions are summarized in Table I. Error in the MD part of the simulation is below 1.4%, whereas in the continuum region, it is below 1.0%, and in the transitional region, the largest error is 0.5%. The table also includes the total number of splitting and combining events over the length of the production run. These two numbers are always approximately the same, implying no net mass transfer between regions at equilibrium. The frequency for splitting and combining particles averaged over the three tests is every 370 time steps. Finally, we include the average solute created or destroyed per splitting and combining event. As before, the net solute exchanged with the bath per refining/coarsening step is approximately zero, i.e., the solute lost on average during particle splitting or combining is balanced by the solute gained from the bath by the reverse move.

. | MD region percent error . | SDPD region percent error . | Transition region percent error . | Number of refining events . | Number of coarsening events . | Average solute change per event (refining) . | Average solute change per event (coarsening) . |
---|---|---|---|---|---|---|---|

〈Φ〉 = 0.4 | 0.65 | 0.54 | 0.43 | 35 197 | 35 197 | −0.0075 | 0.0063 |

〈Φ〉 = 0.5 | 1.43 | 0.51 | 0.45 | 35 062 | 35 064 | −0.0207 | 0.0211 |

〈Φ〉 = 0.6 | 0.40 | 0.98 | 0.31 | 35 146 | 35 140 | 0.0033 | −0.0045 |

. | MD region percent error . | SDPD region percent error . | Transition region percent error . | Number of refining events . | Number of coarsening events . | Average solute change per event (refining) . | Average solute change per event (coarsening) . |
---|---|---|---|---|---|---|---|

〈Φ〉 = 0.4 | 0.65 | 0.54 | 0.43 | 35 197 | 35 197 | −0.0075 | 0.0063 |

〈Φ〉 = 0.5 | 1.43 | 0.51 | 0.45 | 35 062 | 35 064 | −0.0207 | 0.0211 |

〈Φ〉 = 0.6 | 0.40 | 0.98 | 0.31 | 35 146 | 35 140 | 0.0033 | −0.0045 |

Note that while we focus on MD-continuum simulations in this section, this approach is general to any situation involving a particle description for the fluid, even ones where molecular regions are not present altogether. For example, the widely used, bottom-up, coarse-grained particle technique known as dissipative particle dynamics (DPD)^{71,72} has been applied to qualitative multicomponent simulations where DPD particles are assigned a unique identity (e.g., type A or type B).^{73,74} Using the discrete-continuum technique described here, it is possible to couple a binary DPD fluid region to a more coarse continuum one modeled using SDPD. Another possibility is to interface a fine binary fluid composed from particles with potentials obtained from Boltzmann inversion or other bottom-up coarse-graining techniques to a more coarse SDPD continuum fluid.

## VII. NON-EQUILIBRIUM MULTISCALE SIMULATION OF A BINARY MIXTURE

We now extend the tests to non-equilibrium situations, where the chemical potential is no longer uniform across the system, giving rise to fluxes of solute. We specifically focus on the problem of quasi-1D diffusion as a case study, where a concentration gradient is imposed across the system perpendicular to the interface between the MD and continuum parts of the simulation box. The system is set up identically to the one in Sec. VI, with a fluid bounded by two walls [see Fig. 6(a)]. As before, a MD wall is placed next to a MD fluid, which is coupled to a SDPD fluid where particles are six times as massive. Previously we assumed a fixed chemical potential, and thus some equilibrium concentration, when performing the splitting and coarsening steps. Recall that the interface region between the discrete and (larger) continuum SDPD particles is the portion of the simulation box that is coupled to a bath. When concentration gradients are present, we need to adjust the chemical potential of the bath to match the chemical potential that would exist at this location in the absence of multiple degrees of coarse-graining, i.e., if a uniform fluid description was used throughout. In order to generalize this approach to non-equilibrium simulations, we therefore calculate the average concentration in a small bin to the left of the transition region, as well as the average concentration in a bin to the right of this interface, and then perform a linear interpolation to determine the average inside the transition zone (where this average is used for the splitting and refining steps). In other words, the average concentration that appears in Eqs. (13) and (15) is no longer a constant value across the simulation box but is calculated from a linear interpolation across the interface region, and the bath enforces a chemical potential based on this value. This interpolation is performed whenever it is necessary to split or combine particles. Therefore, we make a local equilibrium assumption and linearize the concentration profile near the interface. Note that if the length scale characterizing the concentration gradient is similar to the width of the interface region, this would introduce error.

As in Sec. VI, the probability distribution for the coarsening step is not given by Eq. (15) owing to the small particle size and requires additional care. Moreover, since the average concentration is no longer fixed, the distribution will change as the average concentration inside the interface varies, while the system evolves towards its equilibrium state. Hence, we first set up a series of equilibrium SDPD simulations at a number of different average concentrations at regular intervals in order to numerically determine the concentration probability distribution for each one. Next, we use a truncated Gaussian function to fit these results. Once the fitting parameters for a number of different average concentrations are known, we can interpolate between these values to approximate the solute probability distribution for any average concentration that may appear in the transition region. In other words, we determine $P\Phi $ for any $\Phi $ that may potentially appear within the interface region between the discrete and continuous domains, and therefore we are able to sample the probability distribution appropriate to the instantaneous average concentration inside the transition zone. In our test below, we know that the concentration inside the transition region should only vary between 0.4 and 0.6, and therefore we only pre-calculate distributions for this range, though in principle we could parametrize distributions for the full spectrum of concentrations between 0 and 1. This procedure assumes that the time scale for the evolution of $\Phi $ inside the transition region is slow relative to other time scales (e.g., the motion of the individual particles), and that the bins used to calculate this average are sufficiently large to give an accurate result and remove the effects of statistical nose.

For the present test, the concentration is initially 0.4 across the entire simulation box, and after equilibration, the concentration at the SDPD wall is increased to 0.6, and a gradient begins to propagate across the box. The MD wall naturally provides a no-flux boundary condition. Since the concentration can vary between 0.4 and 0.6 inside the interface region, prior to starting we perform a separate series of equilibrium pure SDPD simulations at the following concentrations: 〈Φ〉 = 0.40, 0.45, 0.50, 0.55, 0.60. For each one, the distribution $P\Phi $ is calculated, fit to a truncated Gaussian function, and the values for the fitting parameters are tabulated. Then, if we are running the full hybrid simulation and the value of the concentration inside the buffer is determined for example to be 0.47, we perform a linear interpolation using the fitting parameters for the 0.45 and 0.50 distributions, giving the distribution for that instant $P\Phi =0.47$. This distribution is used in place of Eq. (15) if it is necessary to perform a coarsening step. The splitting probabilities are the same as before, except the average concentration in Eq. (13) is replaced with the interpolated value. To perform the interpolation and determine the average transition region concentration, the average concentration is calculated using bins of width 5 immediately to the left and to the right of the transition zone. We find that this is a sufficient bin size to give a converged average for the present test.

Prior to switching on the concentration gradient, the system is equilibrated at an average concentration of $\Phi =0.40$ for 1 × 10^{6} time steps. Due to significant statistical and thermal noise, we average the time-dependent concentration profile over 10 separate trajectories, each initialized using a different random seed. The concentration profile at several instants is shown in Fig. 7, compared to the analytical solution from the diffusion equation. The average error per bin over the length of the simulation relative to the analytical solution is 3.02% in the MD part of the simulation box, 2.98% within the transition region, and 2.10% in the continuum part. Note that while we capture the correct overall diffusion dynamics, the numerical solution may advance slightly ahead or behind of the analytical one since the linear interpolation inside the buffer region will overestimate/underestimate the average concentration depending on the curvature of the concentration profile. Choosing a more sophisticated interpolation approach that takes into consideration the instantaneous concentration profile curvature would yield more accurate results, though for simplicity we use a linear approximation in this work. An additional source of error in the solution is that the discrete atomically scaled SDPD particles appearing within the refining and overlap subregions have a slightly higher diffusivity than their MD counterparts,^{9} which is typical for particles having softer coarse-grained interactions. Moreover, particles within the MD-to-SDPD transition zone that feature mixed atomic/continuum interactions exhibit a spatially dependent self-diffusivity. These factors would affect the dynamics of the system being investigated, and it may be necessary to modify the solute dynamics within the transition region (e.g., by adjusting the strength of the coupling to the thermostat^{16}) to maintain a uniform diffusivity, though this is beyond the scope of this paper, and we find that without this correction, the solution is sufficiently accurate for the present case.

## VIII. CONCLUSIONS

In this work, we described an approach for multicomponent hydrodynamic simulations with bridging between high-resolution molecular regions where particles represent individual solute or solvent molecules/clusters, to lower-resolution continuum ones where each particle is a fluid volume that may contain some amount of solute and solvent. This approach is general to any discrete concentration model for the high-resolution region (MD or coarse-grained molecular), and any continuum particle approach for the more coarse-grained part, though in our work we focus on using SDPD. This is achieved by allowing large continuum particles to split into finer discrete ones, and vice versa. The splitting and combining rules are constructed such that an average concentration (i.e., chemical potential) is enforced, allowing for straightforward multiscale equilibrium simulations of fluid mixtures. Moreover, by performing an interpolation of the concentration profile across the interface region separating the fine and coarse-grained domains, it is possible to adjust the bath chemical potential on-the-fly and perform non-equilibrium simulations as well. A simple binning average was used in Sec. VII to determine the concentration profile values to the left and right of the interface region, which gives an estimated value for the concentration inside the transition zone. This approach can be generalized to more complicated concentration profiles with variation in all directions by calculating bin averages in three dimensions or performing local SPH averages instead, though the length scale characterizing gradients in concentration should be larger than the transition zone width to minimize errors due to this local linearization. Importantly, we demonstrated the applicability of this method by performing both equilibrium and non-equilibrium MD-continuum simulations of a binary ideal mixture.

While this approach now allows multicomponent hybrid simulations that may have many useful applications, there are certainly directions along which it can be improved. First, the methodology here is limited to ideal mixing, where the distinction between the solute and solvent is only based on particle identities, and interactions between all particles in the fine (e.g., molecular) region are identical. Of course, ideal mixing is a good approximation for many situations, and this approach is still applicable to any problem where the solute flux is proportional to the gradient of concentration. The simple splitting and refining rules can be adjusted for more complicated non-ideal fluid mixtures if the functional form of the entropy and chemical potential or activity coefficients are known (e.g., a regular solution having an enthalpy of mixing); this would require deriving a new model for the concentration dynamics through GENERIC and concurrently solving an equation for the entropy production. Another potential limitation is that the self-diffusivity of the atomically scaled, discrete SDPD particles in a MD-continuum simulation may not match the diffusivity of the MD atoms precisely since coarse-grained particles typically have softer interactions and therefore larger diffusion constants. This puts a limit on the allowed level of coarse-graining for the SDPD region because an unphysical, negative Fickian diffusion constant is required if the self-diffusivity of the SDPD particles is larger than the MD ones. Finally, as the size of the SDPD particles shrinks, fluctuations in concentration become increasingly large. In such situations, it is no longer possible to use Eq. (15) for the coarsening steps, and instead the probability distribution must be pre-computed numerically prior to starting the multiscale simulation.

Importantly, the framework provided in this paper opens the possibility for a number of applications, including the transport and margination of porous drug particles in intravascular drug delivery.^{24,25} Here particles experience complex interactions due to surface chemistry at the nanoscale, van der Waals and lubrication forces at micron separations, and long-range hydrodynamics over larger scales. Using multiscale multicomponent techniques, it is possible to perform simulations that also capture the dynamics of the active ingredient lost to the surrounding fluid. Another important system is surface nanobubbles and nanodroplets, nanoscale bubbles and drops that form along hydrophobic boundaries and stably exist over time scale orders of magnitude beyond classical predictions.^{32–36} All current models that explain this unusual stability feature multiple length scales and multiscale multicomponent techniques, such as the one in this work, are essential for capturing the key physics at all levels.

While the systems discussed here are uncharged, the adaptive resolution approach for bridging MD and coarse-grained fluid representation is also applicable to problems featuring electrostatic interactions,^{15–17} and coarse particle fluid descriptions similar to the ones used here have been modified to include the evolution of the electrostatic potential for dilute electrolyte solutions.^{75} Therefore, the approach presented in this work is potentially extensible and useful for a host of important electrokinetic phenomena, including corrosion,^{26,27} modeling of biomolecules solvated in ionic solution,^{28} and bioprocesses featuring ion transport (e.g., ion channels).^{29–31} These techniques allow for molecular detail in the region surrounding the surface, protein, or membrane, while having a more coarse-grained description for the bulk fluid further away. Ultimately, multiscale simulation methods have been historically difficult to implement, and we believe this general extension to multicomponent problems provides a basis not only for large-scale modeling of problems involving binary mixtures but also for future advancements in the field towards methods that are more efficient and refined.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge the support of the National Science Foundation (Award No. CBET-1256838) and a Dow Discovery Fellowship. Computational resources were provided by the Center for Scientific Computing from the CNSI, MRL: an NSF MRSEC (No. DMR-1121053) and NSF CNS-0960316. Molecular visualization was performed with the UCSF Chimera package. Chimera is developed by the Resource for Biocomputing, Visualization, and Informatics at the University of California, San Francisco (supported by NIGMS P41-GM103311).