THE JOURNAL OF CHEMICAL PHYSICS 142, 194504 (2015)

Modeling nanoscale hydrodynamics by smoothed dissipative particle dynamics Huan Lei,1,a) Christopher J. Mundy,1,b) Gregory K. Schenter,1,c) and Nikolaos K. Voulgarakis2,d)

1 2

Pacific Northwest National Laboratory, Richland, Washington 99352, USA Department of Mathematics, Washington State University, Richland, Washington 99372, USA

(Received 7 February 2015; accepted 6 May 2015; published online 19 May 2015) Thermal fluctuation and hydrophobicity are two hallmarks of fluid hydrodynamics on the nano-scale. It is a challenge to consistently couple the small length and time scale phenomena associated with molecular interaction with larger scale phenomena. The development of this consistency is the essence of mesoscale science. In this study, we use a nanoscale fluid model based on smoothed dissipative particle dynamics that accounts for the phenomena associated with density fluctuations and hydrophobicity. We show consistency in the fluctuation spectrum across scales. In doing so, it is necessary to account for finite fluid particle size. Furthermore, we demonstrate that the present model can capture the void probability and solvation free energy of nonpolar hard particles of different sizes. The present fluid model is well suited for an understanding of emergent phenomena in nano-scale fluid systems. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4921222]

I. INTRODUCTION

The connection between theoretical frameworks operating at different spatio- and temporal scales is an active area of research that is currently being identified as mesoscale science. Herein, we refer mesoscale science as identifying the correct principles at the lower scale (say atomic) that need to be exported to frameworks at a higher scale or a reduced description (say continuum) to explain new, emergent phenomena that cannot be adequately described by either of the individual frameworks. Identifying the mesoscale principles that enable the connection of well established theoretical frameworks to elucidate hydrodynamic phenomena at the nanoscale will require fundamentally different type of mathematical modeling and numerical approaches than those that are currently present in microscopic and macroscopic frameworks. To this end, the mesoscale principles of surface tension and compressibility have already been identified as the molecular scale phenomena to inform our understanding of the hydrophobic effect.1–4 Specifically, at large length scales, the hydrophobic effect is determined by the depletion of hydrogen bonds in the vicinity of the soluteliquid interface. The natural tendency of a system to reduce large surface tension established at such interfaces can also drive drying phenomena between the solutes interfaces, which result in particles aggregation and phase separation. On the other hand, for small length scales, the hydrophobic effect is determined by a finite compressibility governed by the local interaction between individual fluid molecules that allow for rearrangement of the solvent molecules to accommodate small hydrophobic solutes. The crossover from small to large a)Electronic mail: [email protected] b)Electronic mail: [email protected] c)Electronic mail: [email protected] d)Electronic mail: [email protected]

0021-9606/2015/142(19)/194504/9/$30.00

scale hydrophobicity occurs on nanometer length scales. The hydrophobic effect has been implicated in many important processes such as creating and maintaining phase separation critical in many chemical (e.g., segregation of oil and water) and biological processes (e.g., self-organization of biomolecules). We refer to Refs. 5 and 6 for a review on this subject. At small enough scales, where thermal fluctuations dominate over surface tension, capillary instabilities lead to theoretically interesting and technologically important phenomena, such as spreading of nano-droplets on solid surfaces7 and break up of nano-jets.8 Nanoscale fluid dynamics will be dominated by both hydrodynamic fluctuations and hydrophobicity. Thus, any theoretical description of aforementioned processes must contain the interplay and coupling between compressibility falls between a molecular framework and a continuum description parameterized by transport coefficients and Navier Stokes or Fickian diffusion equation. The multi-facet nature of the interplay between hydrophobicity, thermodynamic fluctuations, and hydrodynamic coupling on different length scales imposes difficulties that make the formulation and solution of the governing equation of hydrodynamics a very challenging task. The situation becomes particularly complex as the size and time scale of fluctuations become comparable to scale size associated with boundary conditions. Several theoretical works have been proposed to bridge the gap between the multiple length scales. Remarkably, the first attempt to account for thermal fluctuations into the deterministic hydrodynamic equations was made by Landau and Lifshitz.9 By imposing a stochastic stress term with the magnitude determined by thermal temperature via fluctuation-dissipation theorem, they established the fluctuating hydrodynamics (FHD) framework and proposed the Landau-Lifshitz-Navier-Stokes (LLNS) equation that determines the spatiotemporal evolution of fluid density fields under the influence of thermal fluctuations.10 Although the

142, 194504-1

© 2015 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 22 May 2015 05:26:48

194504-2

Lei et al.

J. Chem. Phys. 142, 194504 (2015)

initial FHD framework was based on this phenomenological argument, it was later justified and derived directly from kinetic theory11,12 and has been applied to model multiphase and complex fluid systems.13–17 Several numerical methods, such as finite volume method,18–21 finite different method,13,22,23 and Lattice-Boltzmann method24 have been applied to solve the stochastic LLNS equation at small scale. Remarkably, one fundamental question in FHD is whether the equations governing the dynamics could be parameterized to reproduce all-atom MD simulations of molecular fluids. In Ref. 13, it was shown that if the finite size of solvent molecules, d mol, is taken into account, the FHD framework can adequately describe the compressibility of fluids at as small scales as 5 Å. Regarding the hydrophobic interaction, Lum-Chandler-Weeks and subsequently Chandler and co-workers proposed field theories1,2 that incorporate the coupling density fluctuations on both large and small length scales to describe the aqueous solvation of nonpolar hard particle. The field-theory Hamiltonian was numerically resolved by mapping to lattice models.2,25,26 In a similar spirit, Borgis et al.27,28 proposed a molecular density functional theory to model the hydrophobic solvation in water. A limitation of the aforementioned studies is that although they contain the correct description of the hydrophobic effects, the coupling to the hydrodynamic modes is not present. An attempt to couple solvation and hydrodynamics is contained in Ref. 29. As an Eulerian grid discretization, there are numerical artifacts such as negative density15 and spurious currents of fluid.30 In this work, we aim to use a reduced computational model to accurately capture the collective motion of fluids while maintaining molecular fidelity by resolving both hydrophobicity and hydrodynamics. We establish the connection between the framework of molecular detail and the framework of continuum description through a Lagrangian particle fluid model based on smoothed dissipative particle hydrodynamics31,32 (SDPD). We show that by extending the standard SDPD model to impose prescribed surface tension while establishing proper scaling principles, our model can resolve the relevant nanoscale capillary wave spectrum and local density fluctuations. Thus, we will show that the dynamic model presented herein can resolve both the nanoscale interfacial fluctuations and solvation free energy across different length scales. II. MESOSCOPIC FLUID MODEL AND METHOD

In this work, we employ the Lagrangian based framework of the smoothed dissipative particle dynamics31,32 to model the isothermal hydrodynamics. The fluid system is represented as discretized particles following the flow field similar to the smoothed particle hydrodynamics (SPH).33,34 In this model, the number density of each fluid particle is defined by  ni = W (r i j , h), (1) j

  where r i j = ri − r j is the distance between particle i and j, and W (r, h) is normalized bell-shaped smooth function with finite support h. In this work, we adopt the quintic spline function proposed in Ref. 35. The deterministic part of the

dynamics of the fluid particle, derived from the discretized Navier-Stokes equation, is given by r˙ i = vi , (2)  Pi P j ∂W * + + e mi v˙ (d) i j i = − n2 n2j - ∂r i j j , i   1 1 ∂W  5µ + vi j + (vi j · ei j )ei j , 3 j ni n j r i j ∂r i j where ri , vi , mi , and Pi represent the position, velocity, mass, and pressure of particle i, respectively. µ represents the viscosity of the fluid system. Equation (2) is closed by relating pressure Pi with density ρi = mi ni through equation of state,   ( )γ ρ −1 , (3) P = P0 ρ0 where parameter P0 and γ can be chosen to satisfy prescribed sound of speed or density fluctuation. The thermal fluctuation is modeled by imposing random force term in a thermodynamically consistent way via the GENERIC framework,36,37  ˆ i j ei j , mi dv(s) Ai j d W i = j

20 1 1 ∂W Ai j = − µk BT 3 ni n j r i j ∂r i j (

(4)

) 1/2 ,

ˆ i j = (dWi j where k BT is the thermal temperature, d W + dWTij )/2 and dWi j represent a matrix of independent increments of the Wiener process. The combination of Eqs. (2) and (4) recovers the standard SDPD equation in bulk. Moreover, to study the interface fluctuations, we need to impose the proper surface tension near the interface. Rather than directly imposing the surface stress through numerical computation of the surface curvature,38 we impose an additional pairwise force between the fluid particles as is developed in Refs. 39–41,  mi v˙ (p) f (p)(r i j )ei j , (5) i = j,i

where f (p)(r i j ) is short-range repulsive and long-range attractive by taking the form f (p)(r i j ) = f α0 βr i j (Aψr a − ψr b ), ψr0(r) = e



r2 ij 2r 2 0

(6) ,

where r a and r b represent the width of the two Gaussian kernels ψr a and ψr b , respectively. f 0 determines the magnitude of the force and A determines the ratio between the repulsive and attractive force terms. In general, the force term (for A > 1, r a < r b ) is repulsive at short range and attractive at long range. The long-range attractive force term contributes to the interfacial energy and the repulsive force term further helps to prevent the particle clumping together. In this work, we choose A = 32, r a = h/7.0, and r b = 2r a such that the contribution to virial pressure from f (p) is close to zero.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 22 May 2015 05:26:48

194504-3

Lei et al.

J. Chem. Phys. 142, 194504 (2015)

In general, the pairwise force f (p) can be imposed to generate surface tension between multiphase fluids, determined by Tαα + Tβ β − 2Tα β , where Tα β represents the interfacial energy between fluid phases α and β. In this study, we aim to model the nanoscale hydrophobicity (e.g., the void probability, see Sec. III B). We model single fluid phase coupled with vacuum and the interfacial energy is determined by the fluid phase. Following the theory of surface tension by Lord Rayleigh,42 the surface tension can be related to the interfacial energy T (t) of the fluid, e.g., the work required to separate the bulk fluid into two independent parts. Different from the multiphase fluid systems (e.g., Ref. 40), two terms contribute to the interfacial energy T (t): the pressure gradient term and the pairwise force between the SDPD particles, namely, T (t) = T (d) + T ( f ),  h 1 (f ) r 4⟨n⟩2 f (p)(r)dr = ζ⟨n⟩2 f α0 β , T =− π 8 0  1 h P(n(z))dz = λP0, T (d) = − 2 0    h ( z) 2 r W (r)dr , n(z) = ⟨n⟩ 1/2 + 2π 1 − r z

(7a) (7b)

A. Interfacial fluctuations on nanoscale

Studies presented in Ref. 47 show that the SDPD shows consistent scaling of thermal fluctuations in reduced scale in bulk fluid. In this study, we further illustrate the ability of the present model based SDPD to describe nanoscale interfacial fluid fluctuations. Here, we model the water system as a slab of fluid within a volume of 20h × 20h × 16h with periodic boundary conditions in both x and y directions, where the discretization length h varies between 0.994 and 3 nm. For each length scale h, the number density of the SDPD particle is 32.768h−3. SDPD particle corresponds to fluid volume from about 0.31 × 0.31 × 0.31 nm3 to 0.9375 × 0.9375 × 0.9375 nm3. For each snapshot, we construct the instantaneous density field n(r) by

(7c) n(r) = (7d)

where ⟨n⟩ is the average number density of the fluid in bulk, n(z) represents the fluid surface density when the two slabs of fluid are separated by z, and λ ≈ 0.065 for quintic spline W (r) and ζ = π(−Ar a6 + r b6 ). As a brief summary, we model the fluid system via the Lagrangian-based particle model by dv = dv(d) + dv(s) + dv( f ),

Ref. 46, and the shear viscosity µ = 0.88 cP. In general, for surfaces of partial hydrophobic and hydrophilic, we need to further tune the model parameters by imposing proper interfacial energy. The temperature for all simulations is T = 298 K.

(8)

where dv(d), dv(s), and dv( f ) are defined by Eqs. (2), (4), and (5). The surface tension is determined by Eq. (7). Also we remark that for the systems we considered in the present work, f (p) will not substantially modify (within 3%) the compressibility of the fluid system. We refer to Refs. 42 and 43 for derivation of Eq. (7b). In Appendix A, we include a detailed derivation of Eq. (7), appropriate for our single component system. Equation (8) is integrated by the explicitly velocity Verlet algorithm.44 To maintain numerical stability, time step ∆t is chosen to satisfy the constraints,35,45    ∆t ≤ min 0.25h/3cs, 0.125h2/9ν, 0.25 min h/3fs , (9) where ν is the kinetic viscosity; cs is the fluid speed of sound; and f s is the acceleration on any fluid particle.

III. APPLICATION TO WATER

In this section, we apply the fluid model based SDPD described above to study the interplay between hydrodynamic fluctuations and hydrophobicity in water. Numerical results are compared directly with atomistic MD simulation of SPC/E model of water whenever possible. We choose the parameters such that the surface tension of the fluid agrees with the watervapor interface and the compressibility of the fluid agrees with the bulk water. In what follows, we set speed of sound c = 1500 m/s, the surface tension σ = 0.0633 J/m2 following

N 

W (|r − rl | ,r c )nl

(10)

l

on lattice grid r(i,j,k) B (x i , y j , z k ), where r c = h, (x i , y j ) h L and z k = k × dz, k = − dz , = (i, j) × dp, i = 1, 2, . . . , dp h h − dz + 1, . . . , dz . In this study, we set dp = h/3 and dz = 0.01h; we have conducted intensive sensitivity studies and verified that the numerical results are independent of dz given dz < 0.012h. Following Willard and Chandler,48 the instantaneous fluid ˜ y) is defined as the isosurface of the fluid density height h(x, by ˜ y)) = nbulk/2. n(x, y, h(x,

(11)

˜ y), we compute the Fourier modes of h(q) ˆ With  h(x, ˜ y)ei(q x x+q y y)dxd y. For comparison, we also con= L12 h(x, ducted MD a simulation of SPC/E water model in a volume of 20 × 20 × 16 nm3 with 213 120 water molecules at thermal temperature T = 298 K using Nose-Hoover thermostat. The instantaneous fluid height can be obtained from the similar procedure by constructing the density field using Eq. (10) with r c = 0.994 nm. To decrease statistical error, 12 identical independent MD simulations are performed. Regarding computational cost, the running time to create the 1.2 ns trajectory for the MD simulation requires 6144 central processing unit (CPU) hours compared to 1934, 534, and 76 CPU hours for the present model with h = 0.994, 1.25, and 1.988 nm, respectively, on AMD CPU Opteron 6272 [email protected] GHz. The spectrum in Fig. 1(a)  the ensemble average of  denotes 2 ˆ the interfacial fluctuations h(q) obtained from the present simulations for different resolution numbers in direct comparison with the prediction of MD simulations. For low  wave 2 2π ˆ number (i.e., |q| . 5h ), the magnitude scales h(q) ∼ |q|1 2 with excellent quantitative agreement with the prediction from capillary wave theory (CWT).49,50 These results indicate that the modeled aqueous fluid system approaches the continuum regime on the length scale of O(5h). As wave number increases, the spectra (including MD result) gradually deviate from the CWT prediction beyond the wave number qc . This

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 22 May 2015 05:26:48

194504-4

Lei et al.

FIG. 1. (a) Capillary wave spectrum of the present fluid model with varying resolution length h compared to the MD simulation of SPC/E water model. Inset plot: capillary wave spectrum of fluid model with h = 0.994 nm. The ˆ 2 dashed vertical line represents the critical wave number q c where h(q) deviates from CWT prediction. (b) The critical wave number q c obtained from the fluid model with varying spatial resolution h, the solid line represents a fitting curve q c ∼ r1c .

is because at high wave number, the fluid interface is governed by the local organization of the fluid particles on small scale. The capillary wave theory based on the assumption the local interfacial energy solely determined by surface tension and local curvature is not valid at small length scale. The deviation regime further depends on the spatial resolution of the fluid model. For low resolution (large h), it is insufficient to resolve the interface on the length scale comparable to MD simulation. Therefore, the spectra deviate from MD results at relatively small qc (i.e., the SDPD spectra deviate from CWT above qc while MD spectrum is still consistent with CWT). As we increase the resolution, qc increases and spectrum agrees with the CWT prediction at intermediate wave number, as shown in Fig. 1(a). In particular, as the resolution approaches the molecular scales (h = 0.994 nm), qc increases, and the spectrum approaches to the MD simulation results throughout the entire wave number regime. To quantify the fluid density distribution across the interface, we adopt two different definitions of the fluid interface.

J. Chem. Phys. 142, 194504 (2015)

First, we investigate the fluid density n(z) defined by Eq. (10) where z refers to the distance to the average position of the interface noted as Gibbs Dividing surface (GDS). To compare the density field interpolated in smooth kernel function, we also compute the density profile of the MD simulation results by Eq. (10) using h = 0.994 nm. As shown in Fig. 2(a), the fluid density decreases smoothly from bulk density to 0 across the interface, where the width of the transition region further depends on the resolution length of the present model. Furthermore, to unmask the fluid structure which is smeared by the fluctuating interface, we compute a coarse grained fluid particle density n( ˜ z˜) with respect to the distance between the SDPD particle and the instantaneous fluid height h˜ following the prescription of Willard and Chandler.48 In contrast with n(z), the density distribution n( ˜ z˜) exhibits structure within the picture of the instantaneous interface. Moreover, the peak structure is sharper and shifted into the bulk region as we increase the resolution length h. This result can be understood as follows: for larger h, each coarse grained fluid SDPD particle represents a larger number of water molecules. Therefore, the first layer of the fluid particles from the instantaneous interface is shifted into the bulk region. Furthermore, the position of each coarse grained SDPD particle, corresponding to the center of the multiple water molecules, is more localized and therefore exhibits larger density peak. On the other hand, as we increase the resolution, the n( ˜ z˜) approaches the results of the MD simulation, as shown in Fig. 2(b). As expected, we emphasize that both the instantaneous ˜ y) defined by Eq. (10) and the capillary fluid height h(x, wave spectrum depend on the smoothed length h, we adopt to construct the density field in Eq. (11).51,52 Therefore, one question arises naturally: Does the present fluid model conserve the fluctuation spectrum as we refine the resolution of the fluid model? To this end, we construct the density field of fluid systems with different resolution length h using fixed smooth length r c = 1.988 nm by Eq. (10). Instantaneous fluid height is again computed using Eq. (11). Fig. 3 presents the fluctuation spectrum obtained from the present fluid model with resolution length scales h = 0.994, 1.25, 1.58, and 1.988 nm. By taking larger smooth length r c , the interfacial fluctuation is further smeared out and dampened, resulting in lower magnitude of the spectrum at high frequencies. By using the same smooth length to construct the density field, the fluctuation spectra collapse onto a single curve. This result validates the scalability of the present model: as we decrease the resolution length h to molecular level, the present model conserves the fluctuation spectrum on longer length scales. Importantly, the present model conserves the fluctuation spectrum while continuing to exploit fluctuations on a small length scale that cannot be captured by the low resolution model either due to the insufficiency of spatial resolution or other numerical issues, namely, interface overhangs and bubbles as noted in Ref. 48. Remarkably, given the same smooth length, the obtained spectra agree well with the MD simulation results throughout the entire wave number regime. Since the high wave number spectrum is governed by the structure factor of the fluid,51 this result suggests that the present model can capture the compressibility of the target fluid system under proper scaling rules, as will be discussed in detail below.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 22 May 2015 05:26:48

194504-5

Lei et al.

J. Chem. Phys. 142, 194504 (2015)

FIG. 2. (a) Normalized density profile of water near the interface defined as GDS. (b) Normalized density of the fluid particle near the instantaneous interface defined by Eq. (11).

B. Hydrophobicity at small and intermediate length scales

Having verified that the present model captures the fluid fluctuations from the continuum to molecular limit, we further explore whether our coarse grained description of water can be extended to characterize hydrophobicity at the nanoscale. At the molecular scale, the hydrophobic effect is governed by local fluid structure and density fluctuation and can be well characterized by Gaussian field model53 and information theory.3,4 For a small nonpolar hard particle, the solvation free energy of volume v, with a one-basis set approximation, is  2 β∆µ = nbulk v 2/2σ v + ln 2πσ v ,   (12) σ v = d 3r d 3r′ξ(r, r′), v

v

where ξ(r, r ) represents the fluid density variance and can be approximated by ′

2 ξ(r, r′) = nbulkδ(r − r′) + nbulk (g(|r − r′|) − 1),

(13)

where g(r) represents the pair correlation function.

FIG. 3. Capillary wave spectrum of the fluid model obtained from h˜ by computing the density field using fixed smoothed length r c = 1.988 nm.

Following the Gaussian field approximation, we compute the solvation free energy of an nonpolar hard particle (radius R < 0.35 nm) simulated by the present fluid model, as shown in Fig. 4(a). In particular, we consider two different resolution scales h = 0.994 and 1.25 nm, where each SDPD particle corresponds to Nm = 1 and 2 water molecules, respectively. However, unlike the capillary wave spectrum, the solvation free energy directly computed by the present model exhibits further dependence on the resolution length scale. For high resolution scale Nm = 1, the solvation free energy predicted by the present model agrees well with the mean field theory, indicating that the present model can capture the local fluid density fluctuations on molecular scale. In contrast, for lower resolution (e.g., Nm = 2), the solvation free energy shows significant deviations from the MD predictions. To reconcile the above inconsistency, we further examine the definition of solvation free energy. For water molecule, the solvation free energy of a volume v is defined as the energy required to form a cavity of v, where the positions of the water molecules are treated as simple points3,4,53 (e.g., the positions of the oxygen atoms). Following we compute  this definition, p (n p − ⟨n p ⟩)2 within a target the density fluctuation δn = spherical region VT with radius Rv , where n p represents the particle number density within the target volume. Fig. 4(b) shows the numerical results of the present model with scale Nm = 1 and 2. For both cases, the density fluctuation δn p approaches the thermodynamic limit for a large target volume (e.g., large Rv ). In contrast, the density fluctuation δn p exhibits a deviation from the thermodynamic limit prediction for smaller volume. In particular, δn p simulated by Nm = 2 yields the largest deviation from the theoretical prediction, which is also consistent with the underestimated solvation free energy for Nm = 2 as shown in Fig. 4(a). The above result is not unexpected. As shown in Ref. 13, the overestimation of the density fluctuation δn p is due to the ambiguous size of fluid molecules in nanoscale. By approximating the fluid molecules as simple points, we ignore the density contribution from additional particles in the vicinity of the specified volume resulting in the overestimation of the density fluctuation. This suggests that we need to consider

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 22 May 2015 05:26:48

194504-6

Lei et al.

J. Chem. Phys. 142, 194504 (2015)

FIG. 4. (a) Solvation energy of nonpolar hard sphere on molecular scale predicted by Gaussian field model with different modeling scale N m = 1 and N m = 2. The theoretical prediction from Lum-Chandler-Weeks (LCW) theory is also shown for comparison. (b) Local density fluctuation within a spherical volume simulated by the present model with modeling scale N m = 1 and N m = 2. The filled symbols represent the numerical results by approximating the fluid particle as simple points. The open symbols represent the numerical results by defining the fluid density using Eq. (14).

the effective size of the coarse grained SDPD particle when computing solvation free energies. Fig. 5 shows a sketch of three fluid particles near the spherical cavity volume (solid line circle). Particle I represents a fluid particle corresponding to single fluid molecule while particles II and III represent fluid particles corresponding to multiple fluid molecules. For a fluid modeled by particle I, the cavity size can be determined by particle position consistent with a explicit water molecule model (oxygen atoms) and therefore yields consistent solvation free energy. However, for particle II, if we determine the cavity size by the SDPD particle position, the actual “void space” is determined by the “sub-molecule” near the cavity center (the inner dashed-dotted circle), which yields an actual solute cavity size that is smaller than the target cavity size. As a result, the computed solvation free energy for this smaller effective cavity differs from the result of the target cavity. To maintain a consistent cavity size for different coarse grained lengths of the SDPD particle, the effective fluid particle size must be known and taken into consideration as depicted by particle III. This is done in practice shifting the center of the SDPD particle position outward by δR (the outer dashed circle).

FIG. 5. Sketch of the fluid particle near a spherical cavity. Particle I represents a fluid particle corresponding to single fluid molecule. Particles II and III represent fluid particles corresponding to multiple fluid molecules.

To determine the particle size, namely, δR at modeling scale Nm , it is tempting to model this quantity as the “Gaussian core” similar to Eq. (1). However, we note that this definition corresponds to the fluid density field at the hydrodynamic scale.1 Instead, similar to Ref. 13, we determine the effective SDPD fluid particle size by revisiting the microscopic definition of the fluid compressibility. As shown in Fig. 4, δn p deviates from thermodynamic limit prediction using the simple point approximation. On the other hand, the compressibility (density fluctuation) recovers the thermodynamic limit prediction provided we take the molecular size of the fluid into consideration. Therefore, we denote each SDPD fluid particle as sphere with radius r mol and define the fluid density by 1  mol nmol = X (ri ), VT i i (14)   X i = VT ∩ Vimol(ri ) /Vimol(ri ), where X imol represents the fraction of the particle volume overlapping with the target volume VT . The fluid particle size is determined by tuning r mol such that δnmol recovers the local fluid compressibility prediction, as shown in Fig. 4(b). With the obtained fluid particle radius r mol, we approximate the cavity size correction 1 mol δR = α(1 − )r . (15) Nm In principle, the coefficient α further depends on the target cavity size Rv . However, in the present work, we find that α ≈ 1 provides a good approximation for Rv > r mol. With the scaling correction by Eq. (15), we can further validate the present model by investigating the void probability distribution within a cubic volume 1.2 × 1.2 × 1.2 nm3 and T = 298 K. Here, we consider three different scales with h = 0.994, 1.25, and 15.8 nm, corresponding to Nm = 1, 2, and 4, respectively. Fig. 6 shows the normalized fluid particle number distribution P(∆N) where ∆N = (N − ⟨N⟩) / ⟨N⟩, and ⟨N⟩ represents the average particle number within the cubic volume. For comparison, we also present the MD simulation results of extended simple point charge (SPC/E) model from

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 22 May 2015 05:26:48

194504-7

Lei et al.

J. Chem. Phys. 142, 194504 (2015)

is likely due to approximation in Eq. (15) and insufficient accuracy to resolve the probe volume. IV. CONCLUSION

FIG. 6. Normalized fluid particle number distribution within a cubic space 1.2 × 1.2 × 1.2 nm3 simulated by the present model with N m = 1, 2, and 4 compared with the MD simulation result from Refs. 25 and 54.

Refs. 25 and 54. As expected, for Nm = 1, P(∆N) agrees well with the MD simulation result throughout the whole regime. For Nm = 2 and 4, P(∆N) deviates from MD simulations at intermediate regimes (∆N ∼ −0.5) since ⟨N⟩ refers to the number of “point” fluid particle within the cubic space. Nevertheless, for all three cases, the void probability distributions agree well with the MD simulation results, yielding consistent solvation free energies −k BT ln P(−1). In contrast, we also compute the P(∆N) for Nm = 4 without the scaling correction of Eq. (15) (cyan dashed-dotted line), yielding overestimated void probability P(−1). Finally, we compute the solvation free energy of a spherical particle with various radius within the bulk fluid. Fig. 7 shows the results of three different scales similar to the previous example. For all three cases, solvation free energy obtained from the present model agrees well with the MD simulation results. In contrast, solvation free energy computed for Nm = 4 without scaling correction yields underestimated value compared with the MD simulation result, as also consistently observed in Fig. 6. Remarkably, Fig. 7 is also evidence of the limitation of the scaling correction in Eq. (15). We note that if spherical size is comparable with the fluid particle size, the solvation free energy is overestimated. This discrepancy

FIG. 7. Solvation energy of a spherical particle of various radius within bulk water simulated by the present model with N m = 1, 2, and 4 compared with the MD simulation result from Ref. 55.

In this study, we put forth a coarse grained particle-based model based on smoothed dissipative particle dynamics to resolve the interplay between nanoscale hydrodynamic fluctuations and solvent response to nanoscale interfaces stemming from the inclusion of explicit molecular scale interaction. We derive an analytical formula to impose a prescribed surface tension by incorporating an additional pairwise force between the standard SDPD particles. We show that the present model captures the capillary wave spectrum of a molecular model (SPC/E) of water at low and intermediate wave number. For high wave number, the simulated capillary wave spectrum deviates from MD results, where the point of deviation depends on the resolution of the SDPD fluid model. As the model resolution approaches the molecular scale, the capillary wave spectrum is in excellent agreement with MD through the entire spectrum. Furthermore, if we construct the instantaneous density field using fixed smooth length determined by the SDPD size, the capillary wave spectra collapse, further validating that the present model conserves the fluctuation spectrum at large length scales level while preserving the quality of fluctuations at the molecular level as we increase the model resolution. In contrast to the hydrodynamic fluctuations, we show that straightforward scaling of SDPD fluid model fails to capture hydrophobicity due to the inherent ambiguities of the finite size of the SDPD fluid particle. We note that this effect is widely neglected in most of the continuum fluid models based on LLNS. We resolve this inconsistency by relating the density fluctuations of the coarse grained SDPD system to the compressibility of the target fluid. The effective SDPD fluid particle size is determined such that the density fluctuations match with the prediction in thermodynamic limit. Using the derived effective sizes of the SDPD fluid model, the void probability of a cubic volume and the solvation free energy of nonpolar hard sphere are predicted by the present model in good agreement with MD simulation results. The present Lagrangian fluid model is not plagued with the numerical and technical issues associated with its Eulerian based counter-parts.15,30 This coarse grained model is now poised to explore a wide range of applications pertaining to nanoscale hydrodynamics. Future studies will be aimed to study nano-scale fluid transport in bulk and confined geometries such as nano-channels, where both thermal fluctuations and the solvent response to the boundary will play a dominant role. Research problems such as modeling the solid boundary with mixed hydrophobic, semi-hydrophobic, and hydrophilic properties, modeling consistent thermal fluctuation near the boundary, and extending the present fluid model to boundary with complex geometry are worth further exploration. Moreover, the present study based on SDPD is adaptive to integrate other effects such as electrostatic interaction that is important to study models of electrolytes at finite concentration under nano-confinement that are known to be important in many energy related and biological applications.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 22 May 2015 05:26:48

194504-8

Lei et al.

J. Chem. Phys. 142, 194504 (2015)

ACKNOWLEDGMENTS

G.K.S. and H.L. are supported from the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research as part of the Collaboratory on Mathematics for Mesoscopic Modeling of Materials (CM4). C.J.M. is supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences. N.K.V. acknowledges support from the National Science Foundation under Grant No. DMS 1418962. Computation of this research is conducted at supercomputer clusters supported by the PNNL Institutional Computing (PIC) program at Pacific Northwest National Laboratory. We would like to thank Nathan Baker, Xin Bian, Wenxiao Pan, Amish Patel, Rick Remsing, Alexandre Tartakovsky, and John Weeks for helpful discussion.

APPENDIX A: DERIVATION OF SURFACE TENSION FORMULATION

Following Refs. 42 and 43, the surface tension between fluids α and β is σ = Tαα + Tβ β − 2Tα β ,



where f (z) represents the pairwise force between fluid particles by distance z. Therefore, the potential energy between fluids α and β is given by  ∞ ∞ 2π⟨n⟩2 ψ(z ′)dz ′dz, (A3) 0

z

and the interfacial energy is given by  ∞ ∞ Tα β = π⟨n⟩2 ψ(z ′)dz ′dz 0 ∞ z 2 = π⟨n⟩ ψ(z)zdz 0  ∞ π z 4 f (z)dz. = ⟨n⟩2 8 0

 Pi P j ∂W * + + n2 n2j - ∂r i j j , i

since f p depends on the local particle number density.



0

z h

= z

( z) 2 r W (r)dr. 2π ⟨n⟩ 1 − r

Therefore, the total number density is given by    h ( z) 2 n(z) = ⟨n⟩ 1/2 + r W (r)dr , 2π 1 − r z

(A8)

(A9)

which derives Eq. (7d).

APPENDIX B: MD SIMULATION OF SPC/E WATER

In this study, we conduct MD simulation of water using SPC/E model. We simulate fluid volume of 20 × 20 × 16 nm3 with 213 120 water molecules at thermal temperature T = 298 K using Nose-Hoover thermostat. We run 12 identical independent MD simulations. For each simulation, we run 10 000 steps of simulation with time step 0.5 fs to reach thermal equilibrium. Then, we conduct 600 000 steps of simulation with time step of 2.0 fs. We collect sampling data by using the snapshot of the last 400 000 steps. 1K.

Lum, D. Chandler, and J. D. Weeks, J. Phys. Chem. B 103, 4570 (1999). Rein ten Wolde, S. X. Sun, and D. Chandler, Phys. Rev. E 65, 011201 (2001). 3G. Hummer, S. Garde, A. E. Garca, A. Pohorille, and L. R. Pratt, Proc. Natl. Acad. Sci. U.S.A. 93, 8951 (1996). 4G. Hummer, S. Garde, A. E. Garca, M. E. Paulaitis, and L. R. Pratt, J. Phys. Chem. B 102, 10469 (1998). 5D. Chandler, Nature 437, 640 (2005). 6B. J. Berne, J. D. Weeks, and R. Zhou, Annu. Rev. Phys. Chem. 60, 85 (2009). 7B. Davidovitch, E. Moro, and H. Stone, Phys. Rev. Lett. 95, 244505 (2005). 8M. Moseler and U. Landman, Science 289, 1165 (2000). 9L. D. Landau and E. M. Lifshitz, Fluid Mechanics, Course of Theoretical Physics Vol. 6, 2nd ed. (Butterworth-Heinemann, 1987). 10J. M. O. de Zarate and J. V. Sengers, Hydrodynamic Fluctuations in Fluids and Fluid Mixtures (Elsevier, San Diego, CA, 2006). 11M. Bixon and R. Zwanzig, Phys. Rev. 187, 267 (1969). 12R. F. Fox and G. E. Uhlenbeck, Phys. Fluids (1958-1988) 13, 2881 (1970). 13N. K. Voulgarakis and J.-W. Chu, J. Chem. Phys. 130, 134111 (2009). 14N. K. Voulgarakis, S. Satish, and J.-W. Chu, J. Chem. Phys. 131, 234115 (2009). 2P.

(A4)

Regarding the pressure force term, the contribution to the surface tension term cannot be neglected. We need to explicitly consider the contribution from pressure force term. However, we note that Eq. (A2) does not apply to the pressure force term given by fp =−

0

Therefore, the interfacial energy is given by  1 h P(n(z))dz, (A7) T (d) = − 2 0 where derives Eq. (7c). Regarding n(z), we consider the number density of a particle at z0 and above which [z0 z + z0] are vacuum. The density contribution is from two parts, the hemisphere below z0 and the segment above z + z0. Given the number density defined by the smooth kernel W (r), the contributions from the two parts are given by ⟨n⟩ /2 and  R  arccos(z/r ) seg 2 n = r dr (−d cos θ)2π ⟨n⟩ W (r)dr

(A1)

where Tα β is the interfacial energy between fluids α and β. The interfacial energy between α and β, due to the pairwise force f α β , can be derived by following Refs. 42 and 43. Consider a particle “p” of fluid α with distance z above the interfacial plane of fluid β. The potential energy between particles “p” and layer of fluid β is given by 2π ⟨n⟩ ψ(z), where ψ(z) is given by  ∞  z′ ψ(z) = f (z ′′)z ′dz ′′, (A2) z

Alternatively, we derive the pressure contribution as following: we compute the work required to separate a bulk fluid into two separate layers in a quasi-static way. Then, the required force is the pressure applied on the unit area of the fluid layer, e.g., P(z). The required energy to fully separate the fluid layer is given by  ∞ P(n(z))dz. (A6)

(A5)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 22 May 2015 05:26:48

194504-9

Lei et al.

15B. Z. Shang, N. K. Voulgarakis, and J.-W. Chu, J. Chem. Phys. 137, 044117

(2012). 16A. Donev, A. J. Nonaka, A. K. Bhattacharjee, A. L. Garcia, and J. B. Bell, Phys. Fluids 27, 037103 (2015). 17K. Balakrishnan, A. L. Garcia, A. Donev, and J. B. Bell, Phys. Rev. E 89, 013017 (2014). 18M. Serrano and P. Español, Phys. Rev. E 64, 046115 (2001). 19J. B. Bell, A. L. Garcia, and S. A. Williams, Phys. Rev. E 76, 016708 (2007). 20A. Donev, E. Vanden-Eijnden, A. L. Garcia, and J. B. Bell, Commun. Appl. Math. Comput. Sci. 5, 149 (2010). 21A. Donev, J. B. Bell, A. de la Fuente, and A. L. Garcia, Phys. Rev. Lett. 106, 204501 (2011). 22P. J. Atzberger, P. R. Kramer, and C. S. Peskin, J. Comput. Phys. 224, 1255 (2007). 23P. J. Atzberger, J. Comput. Phys. 230, 2821 (2011). 24A. J. C. Ladd, Phys. Rev. Lett. 70, 1339 (1993). 25P. Varilly, A. J. Patel, and D. Chandler, J. Chem. Phys. 134, 074109 (2011). 26S. Vaikuntanathan and P. L. Geissler, Phys. Rev. Lett. 112, 020603 (2014). 27G. Jeanmairet, M. Levesque, and D. Borgis, J. Chem. Phys. 139, 154101 (2013). 28G. Jeanmairet, M. Levesque, R. Vuilleumier, and D. Borgis, J. Phys. Chem. Lett. 4, 619 (2013). 29N. K. Voulgarakis, B. Z. Shang, and J.-W. Chu, Phys. Rev. E 88, 023305 (2013). 30S. Popinet and S. Zaleski, Int. J. Numer. Methods Fluids 30, 775 (1999). 31P. Español, M. Serrano, and H. C. Öttinger, Phys. Rev. Lett. 83, 4542 (1999). 32P. Español and M. Revenga, Phys. Rev. E 67, 026705 (2003). 33L. B. Lucy, Astron. J. 82, 1013 (1977). 34J. J. Monaghan, Rep. Prog. Phys. 68, 1703 (2005). 35J. P. Morris, P. J. Fox, and Y. Zhu, J. Comput. Phys. 136, 214 (1997).

J. Chem. Phys. 142, 194504 (2015) 36M.

Grmela and H. C. Öttinger, Phys. Rev. E 56, 6620 (1997). C. Öttinger and M. Grmela, Phys. Rev. E 56, 6633 (1997). 38X. Hu and N. Adams, J. Comput. Phys. 213, 844 (2006). 39A. Tartakovsky and P. Meakin, Phys. Rev. E 72, 026301 (2005). 40A. M. Tartakovsky and P. Meakin, Adv. Water Resour. 29, 1464 (2006). 41J. Kordilla, T. Geyer, and A. Tartakovsky, “Simulation of film and droplet flow on wide aperture fractures using smoothed particle hydrodynamics,” in Proceedings of the 7th International SPHERIC Workshop, 2012. 42 Lord Rayleigh, in On the Theory of Surface Forces (Philosophical Magazine, 1890), Vol. 3. 43A. Tartakovsky and A. Panchenko, “Pairwise force smoothed particle hydrodynamics model for multiphase flow: Surface tension and contact line dynamics,” J. Comput. Phys. (submitted). 44M. P. Allen and D. J. Tildsley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1987). 45R. Courant, K. Friedrichs, and H. Lewy, Math. Ann. 100, 32 (1928). 46C. Vega and E. de Miguel, J. Chem. Phys. 126, 154707 (2007). 47A. Vázquez-Quesada, M. Ellero, and P. Español, J. Chem. Phys. 130, 034901 (2009). 48A. P. Willard and D. Chandler, J. Phys. Chem. B 114, 1954 (2010). 49F. P. Buff, R. A. Lovett, and F. H. Stillinger, Phys. Rev. Lett. 15, 621 (1965). 50R. Evans, Adv. Phys. 28, 143 (1979). 51F. Sedlmeier, D. Horinek, and R. R. Netz, Phys. Rev. Lett. 103, 136102 (2009). 52A. J. Patel, P. Varilly, S. N. Jamadagni, H. Acharya, S. Garde, and D. Chandler, Proc. Natl. Acad. Sci. U. S. A. 108, 17678 (2011). 53D. Chandler, Phys. Rev. E 48, 2898 (1993). 54A. J. Patel, P. Varilly, and D. Chandler, J. Phys. Chem. B 114, 1632 (2010). 55D. M. Huang, P. L. Geissler, and D. Chandler, J. Phys. Chem. B 105, 6704 (2001). 37H.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 22 May 2015 05:26:48

Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.

Modeling nanoscale hydrodynamics by smoothed dissipative particle dynamics.

Thermal fluctuation and hydrophobicity are two hallmarks of fluid hydrodynamics on the nano-scale. It is a challenge to consistently couple the small ...
691KB Sizes 2 Downloads 8 Views