Non-equilibrium molecular dynamics of steady-state fluid transport through a 2D membrane driven by a concentration gradient

We use a novel non-equilibrium algorithm to simulate steady-state fluid transport through a two-dimensional (2D) membrane due to a concentration gradient by molecular dynamics (MD) for the first time. We confirm that, as required by the Onsager reciprocal relations in the linear-response regime, the solution flux obtained using this algorithm agrees with the excess solute flux obtained from an established non-equilibrium MD algorithm for pressure-driven flow. In addition, we show that the concentration-gradient solution flux in this regime is quantified far more efficiently by explicitly applying a transmembrane concentration difference using our algorithm than by applying Onsager reciprocity to pressure-driven flow. The simulated fluid fluxes are captured with reasonable quantitative accuracy by our previously derived continuum theory of concentration-gradient-driven fluid transport through a 2D membrane [J. Chem. Phys. 151, 044705 (2019)] for a wide range of solution and membrane parameters even though the simulated pore sizes are only several times the size of the fluid particles. The simulations deviate from the theory especially for strong solute--membrane interactions relative to the thermal energy, for which the theoretical approximations break down. Our findings will be beneficial for molecular-level understanding of fluid transport driven by concentration gradients through membranes made from 2D materials, which have diverse applications in energy harvesting, molecular separations, and biosensing.


I. INTRODUCTION
Transport of liquid mixtures and solutions through porous membranes is pivotal to many applications, including water desalination and purification, 1 chemical separations, energy generation 2 and storage 3 , and biological and chemical sensing. 4But inadequate performance of current membranes limits widespread adoption of these technologies. 57][8][9] But gaps in knowledge of fundamental aspects of transport processes in 2D membranes, particularly those driven by the concentration gradients 10 that are central to many applications, hinders predictive design of 2D membranes.
The atomic-scale thickness of 2D membranes confers fundamentally different properties compared with conventional membranes that are highly beneficial.For example, 2D membranes can circumvent the permeabilityselectivity trade-off that often plagues conventional desalination or filtration membranes, 6 whereby selectivity against transport of an unwanted component of a mixture is generally associated with reduced permeability to a desired component, with computer simulations of 2D graphene membranes achieving almost complete salt rejection along with water fluxes orders of magnitude larger than those of current desalination membranes. 11The extreme thinness of 2D membranes also means that transmembrane gradients of e.g.pressure, concentration, or electric potential can be enormous, resulting in huge driving forces for fluid transport.The orders-of-magnitude higher power densities for osmotic power generation of single-layer MoS 2 membranes compared with conven-tional membranes was attributed partly to the massive salt gradient. 12In nanopore-based sensing or sequencing of macromolecules, which provides a cheap, fast, and portable method of chemical or biological sensing, the single-atom thickness of 2D membranes offers the possibility of discriminating single monomers in biopolymers due to the membrane thickness being comparable to the inter-monomer spacing. 6,9Achieving this goal, however, requires precise quantification of ion fluxes both in the absence and presence of the biomolecule.
Although numerical simulations of fluid transport across 2D membranes have been carried out in the presence of concentration gradients using both continuum 13,14 and molecular 11,12,[15][16][17][18][19] models, until our recent work, 20 no analytical theory existed to quantify the relationship between the solution and solute flow through a 2D membrane driven by a solute concentration gradient and basic parameters such as the membrane pore size and strength and range of solute-membrane interactions.This theory was derived from a continuum hydrodynamic model of the fluid and its quantitative accuracy was verified by comparison with computational fluid dynamics simulations. 20However, the continuum fluid approximation can break down when the dimensions of the confining pores becomes comparable to the size of the fluid molecules 21 , which could be an issue for 2D membranes even with large pores, due to their atomic-scale thickness.Furthermore, continuum models cannot easily account for the effects of mechanical flexibility of the membrane, which can be significant for 2D membranes 16,22,23 .They also do not readily predict from first principles some phenomena that can strongly impact nano-scale fluid transport, such as finite solid-fluid friction, inhomogeneous fluid properties, non-electrostatic ion-specific interactions, and non-ideal solutions.
But modelling fluid transport due to a concentration difference across a porous membrane using molecular simulations is not straightforward.This is because the periodic boundary conditions that are generally used in such simulations would create an artificial concentration jump across the periodic boundary when a concentration difference is applied across the membrane.Without constraints, mixing of the solution would occur across the periodic boundary instead of across the membrane.Until recently, a generally suitable molecular algorithm for simulating steady-state transport of liquid mixtures driven by concentration gradients has not existed.6][17][18][19] Such simulations in general do not enable the solution and solute fluxes for a given transmembrane concentration difference to be accurately quantified, since the concentrations change measurably over time in the relatively small systems that can be practically simulated.This can be especially problematic for electrolyte solutions, since electrostatic screening depends on concentration, which could have a significant impact for the large concentration gradients that occur across 2D membranes.
Methods for calculating concentration-gradient-driven transport from equilibrium simulations in the absence of a concentration gradient 24 also have deficiencies, particularly in the case of electrolytes, since they do not account for the effects of spatial variations of concentrationdependent properties such as electrostatic screening.On the other hand, stochastic algorithms that combine molecular dynamics (MD) with a grand canonical Monte Carlo (GCMC) method for maintaining constant reservoir concentrations by inserting, deleting, or exchanging particles in the reservoirs [25][26][27][28] suffer from low acceptance rates of stochastic particle moves at the high densities in liquids and unphysical particle insertions can perturb the flow. 2932]34 Significant advantages of these methods over stochastic algorithms are that they are efficient even at liquid densities and they do not involve unphysical particle insertion or deletion steps.However, most of these algorithms have significant deficiencies that limit their utility, such as not being applicable to complex geometries such as 2D membranes, 33 not explicitly modeling the concenentration difference (and thus not accounting for spatial variations in concentrationdependent properties), 33 or not enabling independent control of differences in solute concentration and pressure across the membrane , [30][31][32] which can lead to spurious effects, particularly at high solute concentrations, as discussed later on.One recent algorithm 34 has been developed that does not suffer from these issues, but it controls the solute chemical potential difference across the membrane, whereas the solute concentration difference is a more useful control parameter for direct comparison with experiments.
In this work, we have modified a previous deterministic NEMD algorithm 31 for constraining the solute concentration difference across a porous membrane in a periodic system to enable independent control of the transmembrane pressure difference.As a proof-of-principle of its ability to simulate and accurately quantify steady-state concentration-gradient-driven fluid transport through 2D membranes, we have applied this algorithm to a model system comprising a binary Lennard-Jones (LJ) liquid mixture and 2D LJ membrane.Using this model system, we have systematically investigated the effects of key parameters such as the membrane pore size, average solute mole fraction, and the strength and range of the solutemembrane interactions on the solution and solute fluxes.We have used these molecular simulation results to evaluate the accuracy of our previously developed continuum theory for predicting concentration-gradient-driven fluid transport through these atomically thin membranes.

A. System details
All MD simulations were performed using LAMMPS (3 Mar 2020 GPU-accelerated version), 35,36 with initial simulation configurations constructed using Moltemplate (version 2.16.1) 37 and visualization of simulation trajectories carried out using OVITO (version 3.0.0). 38The system comprised a binary liquid mixture (solute and solvent) and single-layer planar solid membrane parallel to the x, y plane containing an approximately circular hole centered at the origin (Figure 1(a)).The solid particles were placed in a single layer of the (111) surface of a close-packed face-centred cubic (fcc) lattice (lattice constant √ 2σ), and their positions were fixed throughout the simulations.The Lennard-Jones (LJ) potential, was used for the interactions between particles, where r ij is the distance between particles i and j, and ϵ ij , and σ ij are parameters quantifying the strength and range, respectively of their interactions.The solute-solute, solvent-solvent, and solute-solvent interaction parameters, ϵ ij , and σ ij , were set to the same values, ϵ and σ, respectively, in all simulations, while the parameters for the interactions between solute and solid (wall) particles, ϵ uw and σ uw , respectively, were varied for different simulations.Thus, the solute and solvent particles were identical in all respects except for their interactions with the membrane.The potential was cut-off at a distance of 4σ and all particles had mass m.LJ units are used throughout this work, with masses, distances, energies, temperatures, pressures, and times in units of m, σ, ϵ, ϵ/k B , ϵ/σ 3 , and τ = mσ 2 /ϵ, respectively, where k B is the Boltzmann constant.In all simulations, time integration was done with the velocity-Verlet integrator with a time step of 0.005τ and periodic boundary conditions were applied in all dimensions.Unless otherwise stated, simulations were carried out in the canonical (NVT) ensemble at a temperature of ϵ/k B using a Nosé-Hoover 39,40 style thermostat, 41 with only the velocity components perpendicular to the flow (z) direction thermostatted in NEMD simulations.

B. Constrained concentration-and pressure-difference algorithm
To carry out NEMD simulations of steady-state flow due to concentration and/or pressure gradients, we adapted an algorithm by Khalili-Araghi et al. 31 designed to maintain unequal solute concentrations across a membrane in a system with periodic boundary condition.Specifically, we modified it to enable independent control of concentration and pressure differences.Here, we use the convention that the membrane lies in the x, y-plane at z = 0 in the primary simulation box.
The algorithm in Ref. 31 applies a supplementary constant force f i in the direction perpendicular to the membrane to particles of type i in a small transition region of width d far from the membrane, as illustrated in Fig. 1(b).The algorithm was called the nonperiodic energy step method, since the addition of this force is equivalent to applying a nonperiodic energy step or ramp of size ∆ε i = −f i d across the transition region to the particles, which induces a concentration difference between the two sides of the membrane.The energy step ∆ε i needed to achieve a target ratio of the concentrations c + i0 and c − i0 of species i in the upper and lower fluid reservoirs, respectively, can be estimated from the relationship for a system of non-interacting (ideal) particles at infinite dilution, c + i0 /c − i0 = exp (∆ε i /(k B T )).The applied force required to maintain this concentration ratio in this case is which is used as the initial value of the applied force in the algorithm, i.e. f i (t = 0) = f ideal i .In general, an analytical relationship between the force and concentration difference or ratio does not exist for interacting particles.To account for the effect of particle-particle interactions, the forces are adjusted dynamically to converge to the target concentration ratio according to with where ∆t is the simulation time step, α and τ c are tunable parameters, c + i and c − i are the instantaneous concentrations in the upper and lower reservoirs, respectively, and ⟨• • • ⟩ denotes a time average over the interval (t − τ c , t], i.e. over a duration τ c immediately preceding the current time step.(We use the superscript "KA" to distinguish the equations in Ref. 31 from those in our modified algorithm, which are described below.)The instantaneous concentrations are measured in control regions each of width l b on either side of the membrane and far from it (shown in yellow in Fig. 1(b)) that is sufficiently wide to calculate the average in Eq. ( 4) accurately.From Eq. (3), the concentration ratio is expected to converge to the target ratio c + i0 /c − i0 over a duration on the order of ατ c .To avoid calculating the average in Eq. (3) every time step and to reduce the chance of ∆f i (t) becoming undefined due to c + i or c − i being zero, we have modified the force update algorithm from that in Ref. 31 such that the applied force is updated every τ c time steps instead of every time step and the instantaneous concentrations are averaged over time before taking the logarithm.Thus, we replace Eqs. ( 4) and (3) by and respectively, where the time average ⟨• • • ⟩ is taken over the interval [nτ c , (n + 1)τ c ) for n ≤ t/τ c < n + 1 and n ∈ Z.As in Ref.31, Eq.( 2) is used to initialize the applied force and convergence to the target concentration ratio is expected to occur over a duration on the order of ατ c .In principle, the force update scheme specified by Eqs. ( 3) and (4) from Ref. 31 or our modified version specified by Eqs. ( 5) and ( 6) can be used to constrain the concentration difference across of the membrane of any or all species in a multicomponent mixture.In Ref. 31, the external force was only applied to solute species (electrolyte ions in that case), whereas no external force was applied to the solvent (water).However, in general, applying an external force to solutes in the transition region without doing the same to the solvent will induce a hydrostatic pressure difference across the membrane that will affect the fluid fluxes, as we show below.Thus, an external force must also be applied to the solvent particles to achieve a desired pressure difference.Instead of constraining the solvent concentration using Eq. ( 6) to achieve this goal, we use a force update scheme for the solvent that directly controls the pressure difference, which more closely mimics how applied fields would be controlled experimentally.Thus, while we use Eqs.(5)  and (6) to update the applied force on the solute particles (which we label as particles of type i = u), for the solvent particles (type i = v), we replace Eq. ( 6) by where A is the cross-sectional area of the simulation box, ∆P 0 is the target pressure difference, N v (t + τ c ) is the instantaneous number of solvent particles in the transition region at time t + τ c , ∆P = P + − P − is the instantaneous pressure difference between the control regions on either side of the membrane, and the time average ⟨• • • ⟩ is computed identically to that in Eq (6).(Note that ⟨N v ⟩ may be a better choice than N v (t + τ c ) in Eq. ( 7) as its use would reduce the fluctuations in the applied force.)P ± is calculated from the diagonal components of the per-atom stress tensor summed over atoms in the control regions using the stress/atom compute in LAMMPS. 42We note that this method is known to be inaccurate for computing the local pressure in inhomogeneous systems, 43,44 for which more accurate but more computationally expensive and less widely implemented alternatives [43][44][45] exist.However, by keeping the control regions away from strong inhomogeneities in the fluid, this problem can be mitigated.For the initial applied force on the solvent particles, f v (t = 0) = 0 is used.

C. Simulation details
Two system sizes were simulated.Unless otherwise stated, a 50 × 50 unit-cell membrane and a total of 396 000 fluid particles were used.In addition, some nonequilibrium simulations were carried out with a larger 80 × 80 unit-cell membrane and a total of 1 622 016 fluid particles to verify that the simulations of the smaller system did not suffer from finite-size effects.
For each system size, fluid particles were placed on two cubic lattices of equal size on either side of the solid surface, which initially contained no pore, with the z dimension of the box sufficiently large that the particles did not overlap.All fluid particles were initially set to be solvent particles.The system was initially equilibrated in the isothermal-isobaric (NPT) ensemble for 10 6 time steps at a temperature of ϵ/k B and pressure of ϵ/σ 3 using a Nosé-Hoover 39,40 style thermostat and barostat 41 with only the z dimension barostatted for the fluid particles.The average box length in the z dimension measured over the last 10 5 time steps, by which time the instantaneous box length had plateaued, was used as the box length in all subsequent constant-volume simulations, which were at a temperature of ϵ/k B .The z dimension of the box was deformed at constant velocity over 10 3 time steps to reach this value.All solid atoms within a distance a of the origin were deleted to create an approximately circular pore of radius a in the membrane, as illustrated in Fig. S1 of the supplementary material.
Equilibrium MD and NEMD simulations were carried out for various combinations of the solute-membrane interaction parameters ϵ uw and σ uw , pore radius a, and average solute mole fraction χ, with ϵ uw /ϵ = 0.5, 0.8, 1.2 or 1.5, σ uw /σ = 0.8, 1.2 or 1.5, a/σ = 0 (equilibrium simulations only) 3, 4, 6, or 8, and χ = 0.05 or 0.2.In addition to NEMD simulations of concentrationgradient-driven flow, simulations of pressure-driven-flow without a concentration gradient were also carried out for selected systems using the algorithm in Ref. 46, which is similar in some respects to our constrained concentrationand pressure-difference algorithm, but in which a constant and equal force (f u = f v ) is applied to solute and solvent molecules in the transition region.Most simulations with a constrained concentration difference used a target transmembrane solute concentration ratio of c + u0 /c − u0 = 20, but ratios of 2, 3, and 5 were also used for selected systems to verify linear response of the fluid fluxes to the applied driving force.Unless otherwise stated, the transition region width d, control region width d b , and distance l b between the transition and control regions in the NEMD simulations were all 2σ (see Fig. 1) and the target pressure difference ∆P 0 was zero.Details of the simulated systems and their properties are given in Tables S1-S5 of the supplementary material.
Starting from the final simulation configuration from the previous NPT equilibration step, fluid particles were randomly converted into solute particles to give the desired average solute mole fraction.Additionally, for the NEMD simulations, the fluid particles were converted to achieve the maximum target solute concentration ratio of 20.Then, for each system geometry, a NEMD simulation was carried out using the constrained concentration-and pressure-difference algorithm with ϵ uw = ϵ and σ uw = σ, reducing the target solute concentration ratio at intervals of ∼ 8 × 10 6 time steps to obtain steady-state simulation configurations at each of the desired concentration ratios.Fig. 2 shows the variation of the solute concentration and pressure in the control regions with time from one of these simulations, which illustrates the ability of the algorithm to converge and maintain the concentration and pressure difference at the target values.The final configuration at each target concentration ratio was used as the starting configuration for simulations with other values of ϵ uw and σ uw at that target ratio.All these simulations used α = 50 and τ c = 50τ .
For the NEMD simulations, an automated equilibration detection method 47 was used to determine when each system had reached a non-equilibrium steady state and to estimate the effective number of uncorrelated samples in order to calculate steady-state averages and sta-tistical uncertainties (at the 95% confidence level) of fluctuating variables.Distribution functions such as solute concentration profiles, fluid density profiles, and pressure profiles were calculated only using data after the first 8 × 10 6 time steps of each simulation, which ensured the system was at equilibrium or in a non-equilibrium steady state.
It should be noted that changing solvent particles into solute particles with different solute-membrane interaction parameters results in deviations of the bulk pressure in the fluid reservoirs from the target pressure of ϵ/σ 3 in the NPT equilibration simulations, particularly at high concentrations, with average pressures in the equilibrium simulations varying from 0.93 to 1.00ϵ/σ 3 across the range of systems studied (Table S3 of the supplementary material).However, the bulk solution density remained approximately constant, varying from 0.782 and 0.787σ −3 across the range of solute-membrane interactions.

III. RESULTS AND DISCUSSION
A. Application of constrained concentration-and pressure-difference algorithm Fig. 3 depicts the solute concentration, centerline solute concentration (calculated for solute particles within a distance σ of the axis passing through the pore center), and total (solute + solvent) fluid density profiles perpendicular to the membrane for a system to which our constrained concentration-and pressure-difference algorithm was applied, either with or without enforcing the pressure constraint ∆P 0 = 0 using Eq. ( 7).(Results for the highest concentration ratio without the pressure constraint are not shown because the concentration ratio never converged to a steady state.)The method without the pressure constraint is equivalent to the original constrained concentration-difference algorithm of Khalili-Araghi et al. 31 The system in Fig. 3 had a pore radius a = 6σ, solute mole fraction χ = 0.2, and overall repulsive solute-membrane interactions, as indicated by the solute depletion near the membrane in Fig. 3(a).Qualitatively similar results were obtained for other systems, as illustrated in the supplementary material for a system with same solute-membrane interactions but lower ( χ = 0.05) solute mole fraction in Fig. S2 and for a system with the same solute mole fraction but with attractive effective solute-membrane interactions in Fig. S3.
The solute concentration profiles with and without the pressure constraint in Fig. 3   brane pressure difference, as shown in the pressure profiles for the same system in Fig. 4 when the pressure difference is not constrained.On the other hand, constraining the pressure difference to ∆P 0 = 0 results in equal pressures and total fluid densities on either side of the membrane.Fluid flow driven by the pressure difference polarizes the solute concentration near the membrane, Pressure profile in non-equilibrium simulations with (fv ̸ = 0, solid lines) and without (fv = 0, dashed lines) the transmembrane pressure difference constrained to be zero for the same conditions in Fig. 3.
resulting in the differing solute concentration profiles in Fig. 3(a) and (b) with and without the pressure constraint.As shown by the fitted curves in Fig. 3(a), the centerline solute concentration when both the concentration and pressure difference are constrained is consistent with theoretical predictions under conditions in which there is a concentration difference but no pressure difference and the solute-membrane interaction range is small compared with the pore radius, in which the concentration profile is expected to be an inverse tangent function of the axial coordinate. 20luid fluxes were determined by counting the number of solute and solvent particles crossing the membrane as a function of time, as illustrated in Fig. 5 for two systems to which the constrained concentration-and pressuredifference algorithm was applied with zero pressure difference.(In practice, we measured fluxes at the boundary of the simulation box, which was in the center of the transition region (see Fig. 1), which gives the same result as any plane parallel to the membrane at steady steady due to particle conservation.)The effective solute-membrane interactions in Fig. 5(a) and (b) are repulsive and attractive, respectively, i.e. solute is depleted and enhanced near the membrane relative to the bulk, respectively.The linearity of the curves vs time shows both systems are in the steady state with constant fluid fluxes for most of the simulation (Similar behavior was observed for all systems studied.)Consistent with expectations for systems with a transmembrane concentration difference but no pressure difference and for which solute diffusion dominates advection (Péclet number Pe < 1), the solute flux is in the direction of decreasing concentration, while the total solution flux due to concentration-gradient-driven diffusiosmosis is opposite in direction for membranes that repel and attract the solute, with flow towards increasing and decreasing concentration, respectively. 10he flux of particles of type i was calculated as a numerical derivative, Ṅi = ∆N i /∆t, of the cumulative particle number N i crossing the membrane vs time t in the steady state.The simulation trajectory was divided into intervals of 2 × 10 5 timesteps = 10 3 τ for which the flux was calculated, from which the average flux with statistical uncertainties was obtained using the method described in Sec.II C. (We verified that the calculated uncertainties were insensitive to halving or doubling this time interval.)Fig. 6 shows the total ( Ṅ = Ṅu + Ṅv ) and solute ( Ṅu ) fluxes vs the target concentration ratio for a system with repulsive effective solute-membrane interactions for low ( χ = 0.05) and high ( χ = 0.2) solute mole fractions with or without the ∆P 0 = 0 pressure constraint.Consistent with expectations for this system when no transmembrane pressure difference is applied, the total flux is towards increasing concentration whereas the solute flux is towards decreasing concentration when the pressure constraint is enforced.On the other hand, both the total and solute flux are towards decreasing concentration without the pressure constraint, due pressure-driven flow as a result of the non-zero transmembrane pressure difference.The relative discrepancy between the total solution fluxes with and without the pressure constraint is similar for both solute mole fractions and different concentration ratios, highlighting the general importance of applying the pressure constraint to obtain accurate fluxes in constrained concentrationdifference simulations.The solute flux is significantly different with and without the pressure constraint, but the relative discrepancy is much smaller at the lower solute mole fraction, due to the greater importance of solute diffusion over advection (lower Péclet number) at the lower mole fraction, suggesting that there may be circumstances in which the pressure constraint may not greatly affect the solute flux.
The pressure difference was not constrained in the original constrained concentration-difference algorithm in Ref. 31.The study for which the algorithm was developed focused on measuring the ionic current through biological membrane channels due to a concentration difference for a relatively dilute aqueous KCl electrolyte.A concentration ratio of 0.1:1 mol L −1 , pore radius a ≈ 0.55 nm, and pore length L ≈ 6 nm was considered for the OmpF porin that was studied.To maintain the concentration ratio across the pore an external force f = 0.557 kcal mol −1 was applied to the ions within a d = 2.5 Å transition region.As explained in the next section, this force would create a pressure difference of magnitude |∆P | = N u f /A across the membrane, where N u is the number of ions in the transition region and A is its cross-sectional area.Using A = 123.5 Å × 123.5 Å and the average ion concentration of cu = 1.1 mol L −1 in these simulations gives N u ≈ 50, and thus |∆P | ≈ 2.5 kPa.Ignoring the hydraulic resistance of the pore ends for simplicity, which would reduce the flux further, the total solution flux due to this pressure difference can be estimated from the Hagen-Poiseuille equation, 48 Q = πa 4 ∆P/(8ηL), where η is the solution shear viscosity, which we have taken to be that of pure water, η = 8.94 × 10 −4 Pa s. 49 Using the parameters above gives an estimate of the convective ion flux of cu Q ≈ 10 4 s −1 .The lowest total ionic current that was measured in Ref. 31 was ≈ 10 pA, which gives a lower bound (corresponding to a perfectly ion-selective channel) on the total ion flux of ≈ 6 × 10 7 s −1 .Thus, the convective ion flux due to the induced pressure difference would have been negligible compared with that due to the applied concentration difference in this study, and so the application of a pressure constraint would have made little difference to the results.

B. Transport coefficients: verification of linear response and Onsager reciprocity
From now on, we focus on NEMD simulations using our constrained concentration-and pressure-difference algorithm in which the pressure difference has been constrained to be zero.To quantify the concentrationgradient-driven fluid fluxes for all of the systems studied, we define two transport coefficients -the diffusioosmotic mobility, and the solute permeance, which characterize the total volumetric solution flux Q and solute flux J u = Ṅu , respectively for a given transmembrane osmotic pressure difference ∆Π at temperature T .These definitions follow the notation in our previous work, 20 in which we derived a theory of concentration-gradient-driven flow through 2D membranes for dilute solutions, but the equations above generalize them to arbitrary solute concentrations. 33,50qs. ( 8) and ( 9) reduce to the corresponding equations (Eqs.( 33) and ( 34)) in Ref. 20 in the dilute solution limit, where ∆Π = k B T ∆c u .50 The solution flux can be determined from the simulations as Q = Ṅ ρ , where ρ is the bulk total fluid density, which we calculated as the average of the total fluid density in the upper and lower control regions, i.e. ρ = (ρ + + ρ − )/2.
As described recently for a similar NEMD simulation algorithm, 34 the transmembrane osmotic pressure difference ∆Π and hydrostatic pressure difference ∆P can be calculated by considering the balance of applied forces of the fluid particles in the transition region.Decomposing the applied force on each solute particle and each solvent particle as f u = f + δf u and f v = f + δf v , respectively, such that N u δf u + N v δf v = 0, the force due to the osmotic pressure difference is −A∆Π = N u δf u = −N v δf v , while the force due to the hydrostatic pressure difference is −A∆P = (N u + N v )f = N f , where N u and N v are the number of solute and solvent particles, respectively, in the transition region, N is the total number of fluid particles in the transition region, and A is its cross-sectional area. 34(Note that, following the convention in Sec.II B we have defined ∆Π and ∆P as differences across the membrane rather than across the transition region as was done in Ref. 34, so the sign in the previous equations is opposite that in the equivalent equations in Ref. 34.) From these equations, ∆Π and ∆P can be calculated in terms of the number of particles and applied force on particles of each type in the transition region as and Alternatively, ∆Π and ∆P can be calculated from the solute concentration and pressure, respectively, in the control regions.In this case, ∆P = P + − P − , i.e. the pressure difference evaluated in the constrained concentration-and pressure-difference algorithm.On the other hand, ∆Π can be estimated using the osmotic pressure of an incompressible ideal binary mixture, which can be derived from the entropy of mixing at any solute concentration (assuming the same solute and solvent molecular volume v) to be 33,51,52 where χ is the solute mole fraction and ρ = 1/v is the total fluid density (which is assumed to be independent of χ).This equation reduces to the standard van't Hoff equation, Π = c u k B T , for χ ≪ 1.We have evaluated ∆Π = Π + − Π − from Eq. ( 12) using the solute mole fraction, χ + and χ − , and bulk fluid density, ρ + and ρ − , in the upper and lower control regions, respectively.We have compared ∆Π and ∆P calculated from the applied force balance in the transition region and from the concentration/density or pressure in the control regions in Figs.S6 and S7, respectively, in the supplementary material for all our non-equilibrium simulations.This includes simulations in which the osmotic pressure difference was non-zero and the pressure difference was constrained to be zero, simulations in which the osmotic pressure difference was non-zero and the pressure difference was unconstrained and thus non-zero, and simulations in which the osmotic pressure difference was zero and pressure difference was non-zero.This comparison shows perfect agreement between ∆P calculated using either method, but ∆Π calculated from the control regions using Eq. ( 12) overestimates (by up to ≈15%) the value obtained from the transition region force balance using Eq. ( 10) for larger ∆Π.The origin of this discrepancy is likely the absence of a well-defined "bulk" solute concentration in the upper or lower fluid reservoirs to unambiguously define the osmotic pressure in Eq. ( 12) in simulations with a non-zero concentration difference, since the concentration varies throughout the system (see e.g.Fig. 3).By contrast, even with a transmembrane pressure difference, the pressure profile is essentially flat in either reservoir except in the immediate vicinity of the transition region or membrane (see e.g.Fig. 4), so the ∆P between the control regions is representative of the transmembrane pressure difference.
Using ∆Π from the force balance in the transition region, we have verified that the transport coefficients in Eqs. ( 8) and ( 9) were independent of the system size and transition region width for selected systems in Fig. S8.We have have also verified that the transport coefficients were independent of ∆Π, i.e. the systems were in the regime of linear response of the fluid fluxes to the applied osmotic driving force, for selected systems as shown in Fig. 7 for κ DO and in Fig. S9 of the supplementary material for P s .The three selected systems encompassed those most likely to deviate from linear response, namely those with the three highest solution flux magnitudes (which includes those with the two highest solute flux magnitudes) at the highest target concentration ratio c + u0 /c − u0 = 20 for the pore radius a = 6σ that was used in most of simulations.As indicated by solid symbols in these figures, κ DO and P s are independent of ∆Π up to the highest ∆Π, except for the system with the largest fluxes (ϵ uw = 1.5ϵ, σ uw = 1.5σ, χ = 0.2), for which linear response appears to hold up to ∆Π ≈ 0.15ϵ/σ 3 (c + u0 /c − u0 = 3).Although a few simulations were carried out for systems with a larger pore radius of a = 8σ, for which the fluid fluxes were greater than those with a = 6σ for the same solution and membrane properties, the fluxes were well within the range in which the system with the largest fluxes was still in the linear-response regime.Thus, we can be confident that all systems besides that one were in the linear-response regime for the conditions simulated.
In the linear response regime, the fluid fluxes due to transmembrane differences in the hydrostatic pressure ∆P and osmotic pressure ∆Π (or, equivalently, chemical potential, ∆µ), are given by 53 or using ∆µ = cu ∆Π, 34 with L 11 = Λ 11 , L 12 = Λ 12 /c u , L 21 = Λ 21 /c u , and is the volumetric solvent flux, and cu and cv are the average bulk solute and solvent concentrations, respectively, which we have calculated as the average of the concentrations in the upper and lower control regions, i.e. ci = (c + i +c − i )/2.Note that a number of previous studies on concentration-gradientdriven transport 33,34,50,[55][56][57][58] have not distinguished between the total volumetric solution flux Q and volumetric solvent flux Q v in Eqs.(13) or (14), although this distinction is clear in the derivation by de Groot and Mazur 53 and in equations used in other studies. 59,60For the dilute solutions investigated in most of these studies, this distinction is not important, since Q ≈ Q v in this regime, but Q and Q v can differ significantly at high solute concentrations such as studied here, especially when the solute and solvent fluxes are in opposite directions.
From the definitions of the transport coefficients in Eqs. ( 8) and ( 14), L 12 = κ DO /(k B T ), and given that T = ϵ/k B in all our simulations, this means that L 12 and κ DO have the same numerical value in reduced LJ units in this study.For two of the systems for which L 12 was measured in Fig. 7, we also measured L 21 using the definition in Eq. ( 14) in simulations of pressure-driven flow in the absence of a transmembrane concentration difference for several pressure differences using the NEMD algorithm in Ref. 46.These results are presented in Fig. 7, and show that the L 12 = L 21 reciprocial relation is ver-ified, at least for the lowest applied pressures, with the consistency between our algorithm and the previously established and widely applied method in Ref. 46 demonstrating the validity of our method for quantifying nonequilibrium concentration-gradient-driven transport.Fig. 7 also shows that the pressure-driven flow simulations deviate from linear response at much lower values of ∆P than the ∆Π values at which deviations occur in the concentration-gradient-driven flow simulations.This means that much longer simulations were required to obtain roughly comparable statistical uncertainties in the linear-response regime for the pressure-driven flow simulations compared with the concentration-gradient-driven flow simulations, highlighting the computational benefits of directly measuring the diffusioosmotic mobility in NEMD simulations with an applied concentration gradient, at least for 2D membrane systems.

C. Comparison of simulation vs theory
We have previously derived a theory of fluid transport through a circular pore in an infinitesimally thin planar membrane due to a transmembrane concentration difference 20 by solving the continuum hydrodynamic (Stokes, advection-diffusion, and continuity) equations for low-Reynolds-number steady-state flow of a dilute solution of an incompressible Newtonian fluid, under the assumptions that solute diffusion dominates solute advection (Péclet number Pe ≪ 1) and that the effective solute-membrane interaction potential U is small compared with the thermal energy k B T .This theory is straightforwardly generalized to arbitrary solute concentrations, by analogy with a related theory derived for concentration-gradient-driven fluid transport parallel to a planar surface at high solute concentration: 33,50 the equations derived in Ref. 20 for the transport coefficients quantifying the fluid fluxes at low concentration apply at high concentration, while the distinction between low and high solute concentrations is manifested in the concentration dependence of the osmotic pressure driving force in Eqs. ( 8) and (9).Thus, the diffusioosmotic mobility κ DO quantifying the total solution flux is 20 and the solute permeance P s quantifying the solute flux (evaluated at the pore mouth at z = 0) is 20 where a is the pore radius, η is the solution shear viscosity, D is the solute diffusivity, and the oblate-spheroidal coordinates ζ and ν are defined in terms of the radial and axial coordinates by r = a (1 + ν 2 )(1 − ζ 2 ) and z = aνζ, respectively.The effective solute-membrane interaction potential U , which includes contributions both from direct solute-membrane interactions and indirect solvent-mediated interactions, is defined by where c(ζ, ν) is the solute concentration distribution and c u∞ (ζ, ν) is a hypothetical solute concentration distribution for the same boundary conditions but with U = 0.The theory yields simple scaling relationships for the transport coefficients as a function of the pore radius a and strength ϵ and range λ of U in the limits of weak interactions (ϵ ≪ k B T ) and a small (λ ≪ a) or large (λ ≫ a) interaction range relative to the pore size. 20lthough none of our simulations correspond strictly to the λ ≪ a or λ ≫ a limits, λ is significantly smaller than a for all but the smallest pore studied.Our simulation results are consistent with the predicted scaling in the λ ≪ a limit, in which κ DO and P s are both expected to be proportional to the pore radius a, 20 as shown in Fig. 8. Fig. 8 also shows that κ DO is approximately independent of the average solute mole fraction for the simulated conditions, whereas P s depends significantly on the average solute mole fraction.The behavior of P s appears to be due to a subtle interplay of the opposing diffusive and advective solute fluxes for the systems in Fig. 8, which is not captured by the theory as it assumes that solute advection is negligible.
The dependence of κ DO and P s on the solutemembrane interaction strength and range parameters, ϵ uw and σ uw , are given in in the supplementary material in Figs.S10 and S11, respectively.It should be noted that the strength and range of U is not simply proportional to these parameters, due to the complex manybody contributions to the effective interactions.Interestingly, Fig. S10 shows a linear dependence of κ DO and P s on ϵ uw for fixed σ uw , which is consistent with the predicted scaling with the effective solute-membrane interaction strength ϵ for ϵ ≪ k B T , even though the direct solute-membrane interactions are certainly not weak in all cases.The dependence of κ DO and P s on σ uw for fixed ϵ uw shown in Fig. S11 is more complex, in part because σ uw controls not only the range but also the strength of the direct solute-membrane interactions, since a given solute particle interacts with more membrane particles as σ uw increases; thus, P s varies non-monotonically while κ DO even changes sign with increasing σ uw .
We have used Eqs.( 15)-( 17) to predict κ DO and P s for all the simulated systems without using any information from the NEMD simulations.We used analytical equations fitted to equilibrium MD simulation data for the diffusivity and shear viscosity of the LJ fluid over a wide range of density and temperature for the same interaction cutoff distance in our simulations 61 to obtain D = 0.0697σ 2 /τ and η = 1.84ϵτ /σ 3 for the temperature T = ϵ/k B and total bulk density ρ ≈ 0.787σ −3 in all the simulations.To obtain U from Eq. ( 17), we used the solute concentration profile from an equilibrium MD simulation with the same solute-membrane interaction parameters, pore radius, and average solute mole frac- tion as the NEMD simulation for which the transport coefficients were being predicted.In this case, c u∞ (ζ, ν) in Eq. ( 17) is a constant and equal to the bulk solute concentration in the equilibrium simulation.Distributions of the solute concentration and total fluid density in all the equilibrium simulations, in one dimension (1D) as a function of the axial (z) or radial (r) coordinate and in two dimensions (2D) as a function of both r and z, are given in Figs.S29-S38 of the supplementary material.The integrals in Eqs. ( 15) and ( 16) were computed using the quad and nquad functions, respectively, in the SciPy Python package and 2D solute concentration distributions were interpolated using a bivariate spline with the RectBivariateSpline SciPy function.
The pore radius a in Eqs. ( 15) and ( 16) corresponds to where the hydrodynamic boundary conditions are applied in the continuum theory and does not necessarily correspond to the definition of the pore radius used up to this point, which was the distance from the pore center within which the centers of solid atoms were absent.The Gibbs dividing surface from equilibrium MD simulations has previously been shown previously to describe the hydrodynamic boundary position in NEMD simulations of fluid flow accurately, 62 and so we have used this prescrip-tion to define an effective pore radius a h to replace the actual pore radius a in Eqs. ( 15) and ( 16), given by (18)   where the total fluid density ρ(r, z = 0) in the plane of the membrane pore at z = 0 was obtained from the 2D equilibrium distribution ρ(r, z) by the same bivariate spline interpolation described above for the solute concentration distribution.(In practice the second integral was calculated up to finite value of r beyond which the total fluid density ρ(r, z = 0) was zero.)As shown in Fig. S39 of the supplementary material, a h is within a few percent of a, so either a or a h could be used in Eqs. ( 15) and ( 16) with little difference.(The dependence of a h on the solute-membrane interaction strength and range parameters, ϵ uw and σ uw , for fixed a are also given in the supplementary material in Figs.S40 and S41.) Fig. 9 compares the diffusioosmotic mobility κ DO from the simulations with that calculated from the theory for all the simulated systems (which includes variations in pore radii, average solute mole fraction, and the strength and range of solute-membrane interactions) except for those with the strongest attractive solute-membrane interactions (ϵ uw = 1.5ϵ, σ uw = 1.5σ).The latter data are not included because they occur on a very different scale from the rest of the data, making visualization on the same figure difficult, and because the assumption in the theory of weak solute-membrane interactions certainly breaks down for these systems.A comparison of all the simulated systems in the linear-response regime is given in Fig. S12 of the supplementary material.The simulations and theory are also compared in Fig. S13 of the supplementary material for calculations using the actual pore radius a instead of the effective pore radius a h in the theory, showing that this choice makes little quantitative difference.
There is good quantitative agreement between the theory and simulations for most of the data in Fig. 9 , although deviations are evident for larger magnitudes of κ DO , especially for positive κ DO .Discrepancies between the theory and simulations are not entirely surprising, given the number of approximations made in the theory, namely weak effective solute-membrane interactions, a small Péclet number, and an infinitesimally thin planar membrane.The discrepancies appear to be largely due to the break down of the assumption of weak effective solute-membrane interactions, as indicated by the correlation between the deviation of the theory from the simulation results and the surface solute excess Γ used to color the data points in Fig. 9, which quantities the degree of adsorption (Γ > 0) or depletion (Γ < 0) of the solute at the membrane surface relative to the bulk.Γ in Fig. 9 was calculated from the solute concentration profile c u (z) perpendicular to a membrane containing no pore in an equilibrium MD simulation of a system with otherwise identical fluid and membrane properties to the

NEMD simulation using
where, in practice, the upper integration limit in Eq. ( 19) was taken to be the maximum value of z in the simulation box (this choice was not crucial as c u (z) → c u∞ within the simulation box).
On the other hand, there does not appear to be a clear correlation between the discrepencies between the theory and simulation for κ DO and other potentially relevant parameters such as the Péclet number Pe, pore radius a, or solute-membrane interaction strength and range parameters, ϵ uw and σ uw , as shown in Figs.S14-S17 of the supplementary material, in which the data points have been colored by the value of these parameters.If anything, the discrepancies show the opposite trend vs Pe to that expected, with the deviation between theory and simulation increasing with decreasing Pe.We estimated Pe, which measures the relative magnitude of solute advection to solute diffusion, by where we have used the solvent volumetric flux Q v instead of the total volumetric flux Q to quantify the advective solute flux because the total flux Q includes the diffusive component.The assumption of an infinitesimally thin membrane is expected to become less accurate as the aspect ratio of the pore decreases, which corresponds to decreasing pore radius.As noted above, although the effective solute-membrane interaction strength depends on ϵ uw and σ uw , the dependence on either parameter is not straightforward, and thus the lack of correlation of the discrepancies in κ DO with either parameter is not unexpected.
Fig. 10 shows a similar comparison between theory and simulation to Fig. 9, but for the solute permeance P s .Since solute advection is assumed to be neglible (Pe ≪ 1) in the theory, but is clearly is significant in many of the simulations, for which 0.02 ≲ Pe ≲ 0.5, in Fig. 11 we have also compared P s from the theory with the solute permeance from the simulations calculated from the diffusive flux only, which we define as As with Fig. 9, we have excluded data for the strongest attractive solute-membrane interactions (ϵ uw = 1.5ϵ, σ uw = 1.5σ) simulated from Figs. 10 and 11, but the corresponding plots containing all the simulation data can be found in Figs.S18 and S19, respectively, in the supplementary material.As for κ DO the use of the actual pore radius a or effective pore radius a h in the theory makes little difference, as shown in Fig. S20 in the supplementary material.
The agreement between the theory and simulations is not as good for the solute permeance as for κ DO , but the agreement appears to improve by excluding the advective flux from the simulation definition of the solute permeance.As for κ DO , the discrepancies between the the- ory and simulation appear to be most strongly correlated with the strength of the effective solute-membrane interactions as quantified by the surface solute excess used to color the symbols in Figs. 10 and 11 (similar plots for all the simulated systems colored by the Péclet number Pe, pore radius a, or solute-membrane interaction strength and range parameters, ϵ uw and σ uw , are given in Figs.S21-S28 in the supplementary material).

IV. CONCLUSION
We have developed a constrained concentration-and pressure-difference algorithm for non-equilibrium molecular dynamics simulations of steady-state fluid transport driven by concentration and/or pressure differences across a porous membrane in a system with periodic boundary conditions.Our algorithm adapts a previous algorithm by Khalili-Araghi et al., 31 which controls the transmembrane concentration difference by applying an external force to solute particles in a transition region far from the membrane, by also applying an external force to solvent particles in the transition region to control the transmembrane pressure difference.Applying this algorithm to a model system comprising a binary Lennard-Jones liquid mixture and a 2D Lennard-Jones membrane containing a circular pore, we have simulated steady-state concentration-gradient-driven fluid tranport across a 2D membrane with molecular resolution for the first time, enabling accurate quantification of the solution and solute fluxes due to a given applied concentra-tion difference.We have shown that the application of the pressure-difference constraint has a significant effect on the solute concentration distribution across the membrane and the steady-state solution flux due to a transmembrane concentration difference for both low and high average solute concentrations, although the solute flux is less affected by the pressure-difference constraint at low solute concentrations.We have also shown that the solution flux due to an applied concentration difference generated by our algorithm is consistent with Onsager reciprocity in the linear-response regime by comparison with fluid fluxes due to an applied pressure difference.Furthermore, we have shown that directly simulating a transmembrane concentration difference is far more efficient for quantifying the concentration-gradient-driven solution flux in the 2D membrane systems studied than the indirect approach of applying the Onsager reciprocal relations to pressure-driven flow simulations.This is because only very small pressure differences can be applied before the pressure-driven fluid fluxes deviate from linear behavior.Finally, we have shown that our recently developed theory of concentration-gradient-driven flow across a 2D membrane, 20 although derived for a continuum fluid model, gives reasonably good quantitative agreement with the molecular simulations, especially for the total fluid flux, demonstrating its utility for quantifying the fluid transport even in molecular systems.Nevertheless, deviations from the simulation results are evident particularly for strong solute-membrane interactions, for which the assumptions of the theory break down.

SUPPLEMENTARY MATERIAL
The supplementary material contains details of the parameters and properties of all the simulated systems, additional non-equilibrium simulation results (concentration, density, and pressure distributions; comparison of transmembrane osmotic pressure and hydrostatic pressure differences calculated by different methods; verification of linear response for the solute permeance; verification of the independence of measured transport coefficients of system size and transition region width; additional plots of transport coefficients vs system parameters; additional comparisons of transport coefficients obtained from simulation and theory), and solute concentration and total density distributions from equilibrium simulations.

This work was supported by the Australian Research
Council under the Discovery Projects funding scheme (Grant No. DP210102155).This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government, and from the University of Adelaide's Phoenix High-Performance Computing service.

Electronic Supplementary Information
Non-equilibrium molecular dynamics of steady-state fluid transport through a 2D membrane driven by a concentration gradient Daniel J. Rankin and David M. Huang Department of Chemistry, School of Physics, Chemistry and Earth Sciences, The University of Adelaide, SA 5005, Australia

S1 Systems
Table S1.System sizes in equilibrium simulations: L x/y = system x and y dimensions, L z = system z dimension, a = pore radius, χ = average solute mole fraction, N u = number of solute particles, N v = number of solvent particles, and N w = number of solid particles.Table S3.Systems properties in equilibrium simulations: a = pore radius, ε uw = solute-membrane interaction strength parameter, σ uw = solute-membrane interaction range parameter, χ = average solute mole fraction, t total = total simulation duration, and a h = effective pore radius (calculated from Gibbs dividing surface).The bulk solute concentration c u∞ , bulk total fluid density ρ ∞ , and bulk pressure P ∞ were calculated by averaging the axial profiles of the solute concentration (c u (z)), total fluid density (ρ(z)) and pressure (P(z)) for distances > 30σ from the membrane, where the profiles were constant.S4.Systems properties in non-equilibrium concentration-gradient-driven flow simulations: d = width of transition region, L x/y = system x and y dimensions, a = pore radius, ε uw = solute-membrane interaction strength, σ uw = solutemembrane interaction range, χ = average solute mole fraction, c + u0 /c − u0 = target concentration ratio, ∆P 0 = target pressure difference, t total = total simulation duration, ρ = average total fluid density in control regions, cu = average solute concentration in control regions, ∆c = average concentration difference between control regions, and ∆Π and ∆P = average osmotic pressure difference and average pressure difference between reservoirs calculated from the applied force balance in the transition region.S5.Systems properties in non-equilibrium pressure-driven flow simulations: d = width of transition region, L x/y = system x and y dimensions, a = pore radius, ε uw = solute-membrane interaction strength, σ uw = solute-membrane interaction range, χ = solute mole fraction, f u = force on solute particle, f v = force on solvent particle, t total = total simulation duration, ρ = average total fluid density in control regions, cu = average solute concentration in control regions, ∆c = average concentration difference between control regions, and ∆Π and ∆P = average osmotic pressure difference and average pressure difference between reservoirs calculated from the applied force balance in the transition region.S6.Average osmotic pressure difference calculated from the applied force balance in the transition region vs that calculated from the concentration difference between control regions using Eq. ( 12) of the main paper for the high ( χ = 0.2, circles) and low ( χ = 0.05, squares) average solute mole fractions for concentration-gradient-driven flow simulations with (blue symbols) and without (orange symbols) the pressure different across the membrane constrained to be zero and for pressure-driven flow simulations without an imposed concentration difference across the membrane (green symbols).Error bars are smaller than the symbols.  .Solute permeance vs osmotic pressure difference ∆Π across the membrane from concentration-gradient-driven flow simulations with target pressure difference ∆P 0 = 0 for the high average solute mole fraction ( χ = 0.2) and a membrane with pore radius a = 6σ and various solute-membrane interaction parameters: ε uw = 1.5ε, σ uw = 1.5σ (circles); ε uw = 0.5ε, σ uw = 0.8σ (squares); ε uw = 0.5ε, σ uw = 1.5σ (triangles).The translucent bands are the 95% confidence intervals for the lowest value of ∆Π for each set of solute-membrane interaction parameters.The fact that all the data points in each set lie within the corresponding band, except for the circle for the highest ∆Π, indicates that all the systems except that one are in the linear-response regime.          .Effective pore radius vs solute-membrane interaction range parameter σ uw for fixed interaction strength parameter (ε uw = 0.5ε) and pore radius (a = 6σ ) and different average solute mole fractions ( χ = 0.2 (circles) and χ = 0.05 (squares))..

FIG. 1 .
FIG. 1.(a) Snapshot of a typical MD simulation system, consisting of a binary LJ liquid mixture (solvent particles are translucent) and a single-layer planar fcc (111) lattice membrane with a circular pore.(b) Schematic of constrained concentration-and pressure-difference algorithm: a force fi is applied perpendicular to the membrane to each particle of type i (solute or solvent) within a transition region (red) of width d; the concentration ratio and pressure difference is measured between control regions (yellow) of width d b on either side of the membrane a distance l b from the transition region; the applied forces are dynamically adjusted to converge the concentration ratio and pressure difference to target values.

FIG. 2 .
FIG. 2. (a) Solute concentration cu and (b) pressure P in the upper (+) and lower (−) control regions vs time when applying the constrained concentration-and pressure-difference algorithm for a target concentration ratio c + u0 /c − u0 = 20 and target pressure difference ∆P0 = 0, starting from a configuration in which solutes were uniformly distributed on either side of the membrane at the target concentration ratio (a = 6σ, ϵuw = ϵ, σuw = σ, χ = 0.2).
Fig.3depicts the solute concentration, centerline solute concentration (calculated for solute particles within a distance σ of the axis passing through the pore center), and total (solute + solvent) fluid density profiles perpendicular to the membrane for a system to which our constrained concentration-and pressure-difference algorithm was applied, either with or without enforcing the pressure constraint ∆P 0 = 0 using Eq.(7).(Results for the highest concentration ratio without the pressure constraint are not shown because the concentration ratio never converged to a steady state.)The method without the pressure constraint is equivalent to the original constrained concentration-difference algorithm of Khalili-Araghi et al.31The system in Fig.3had a pore radius a = 6σ, solute mole fraction χ = 0.2, and overall repulsive solute-membrane interactions, as indicated by the solute depletion near the membrane in Fig.3(a).Qualitatively similar results were obtained for other systems, as illustrated in the supplementary material for a system with same solute-membrane interactions but lower ( χ = 0.05) solute mole fraction in Fig.S2and for a system with the same solute mole fraction but with attractive effective solute-membrane interactions in Fig.S3.The solute concentration profiles with and without the pressure constraint in Fig.3(a) and (b) are dramatically different.At first glance, the total fluid density profiles in Fig.3(c) are similar.But closer inspection, as shown in the inset, reveals a transmembrane density difference without the pressure constraint, which is induced by the net force exerted on the solution by the applied force on the solute particles used to constrain the concentration difference.The net force is manifested in a transmem-

FIG. 3 .
FIG. 3. (a) Solute concentration, (b) centerline solute concentration, and (c) total fluid density vs z coordinate in non-equilibrium simulations with (fv ̸ = 0, solid lines) and without (fv = 0, dashed lines) the transmembrane pressure difference constrained to be zero for a = 6σ , ϵuw = 0.5ϵ, σuw = 0.8σ, χ = 0.2, and various target solute concentration ratios c + u0 /c − u0 .The inset in (a) zooms in on z values near the membrane, whereas the inset in (c) zooms in on density values around the bulk density.The dotted black lines in (b) are least-squares fit of the solid line to a function of the form b0 + b1 tan −1 (z/b2) with fit parameters b0, b1, and b2.
Fig. S1.Snapshots of membrane pores of the various radii a that were simulated.

Fig
Fig. S2.(a) Solute concentration, (b) centerline solute concentration, and (c) total fluid density vs z coordinate in concentration-gradient-driven flow simulations with ( f v ̸ = 0, solid lines) and without ( f v = 0, dashed lines) the transmembrane pressure difference constrained to be zero for a = 6σ , ε uw = 0.5ε, σ uw = 0.8σ , χ = 0.05, and c + u0 /c − u0 = 20.The inset in (a) zooms in on z values near the membrane, whereas the inset in (c) zooms in on density values around the bulk density.The dotted black line in (b) is a fit of the solid line to a function of the form b 0 + b 1 tan −1 (z/b 2 ) with fit parameters b 0 , b 1 , and b 2 .

Fig. S3 .
Fig. S3.(a) Solute concentration, (b) centerline solute concentration, and (c) total fluid density vs z coordinate in concentration-gradient-driven flow simulations with the transmembrane pressure difference constrained to be zero for a = 6σ , ε uw = 1.5ε, σ uw = 1.5σ , χ = 0.2, and various target solute concentration ratios c + u0 /c − u0 .The inset in (a) zooms in on z values near the membrane, whereas the inset in (c) zooms in on density values around the bulk density.The dotted black lines in (b) are least-squares fit of the solid line to a function of the form b 0 + b 1 tan −1 (z/b 2 ) with fit parameters b 0 , b 1 , and b 2 .

Fig
Fig.S4.Pressure profile in pressure-driven flow simulations without an imposed concentration difference across the membrane for various values of the applied force f u = f v on solute and solvent molecules for the high average solute mole fraction ( χ = 0.2), repulsive effective solute-membrane interactions (ε uw = 0.5ε, σ uw = 0.8σ ), and pore radius a = 6σ .

Fig. S5 .
Fig.S5.Pressure profile in pressure-driven flow simulations without an imposed concentration difference across the membrane for various values of the applied force f u = f v on solute and solvent molecules for the high average solute mole fraction ( χ = 0.2), attractive effective solute-membrane interactions (ε uw = 1.5ε, σ uw = 1.5σ ), and pore radius a = 6σ .
Fig.S7.Average pressure difference calculated from the applied force balance in the transition region vs that calculated between control regions for the high ( χ = 0.2, circles) and low ( χ = 0.05, squares) average solute mole fractions for concentration-gradient-driven flow simulations with (blue symbols) and without (orange symbols) the pressure difference across the membrane constrained to be zero and for pressure-driven flow simulations without an imposed concentration difference across the membrane (green symbols).Error bars are smaller than the symbols.

Fig
Fig. S8.(a) Diffusioosmotic mobility and (b) solute permeance vs pore radius a for fixed average solute mole fraction ( χ = 0.2) and solute-membrane interaction parameters (ε uw = 0.5ε, σ uw = 0.8σ ), and different reservoir dimensions L x/y and transition region widths d, showing that these transport coefficients are independent of L x/y and d.

Fig
Fig. S18.Solute permeance from simulation vs theory for all systems for c + u0 /c − u0 = 20 except for those with the strongest attractive solute-membrane interactions and high solute mole fraction (ε uw = 1.5ε, σ uw = 1.5σ , χ = 0.2), for which c + u0 /c − u0 = 2 was used instead.Symbols are colored by the surface solute excess Γ (units: σ ) and different symbol shapes distinguish high ( χ = 0.2, circles) and low ( χ = 0.05, squares) average solute mole fractions.The simulation and theory values are equal along the solid line.
Fig.S39.Effective pore radius (calculated from the Gibbs dividing surface in equilibrium simulations) vs pore radius (radius within which there were no membrane atom centers) for different average solute mole fractions ( χ = 0.2 (circles) and χ = 0.05 (squares)).

Table S2 .
System sizes in non-equilibrium simulations: a = pore radius, χ = average solute mole fraction, L x/y = system x and y dimensions, L z = system z dimension, N u = number of solute particles, N v = number of solvent particles, and N w = number of solid particles.