0 Applications of All-Atom Molecular Dynamics to Nanofluidics

In recent years, due to the progress in the fabrication of micro and nanodevices, nanofluidics has become an intense research field. The interest of the scientific community is evident from the large amount of published papers and by the issuing of dedicated journals. In this scenario, a better understanding of the key aspects of fluid motion in nanoscale systems is fundamental and hence computational techniques able to provide a deeper insight in the complex phenomena involved in nanoscale mass transport play a crucial role in nanofluidic research. The aim of this chapter is to introduce the reader to the application of classical all-atom molecular dynamics (MD) to nanofluidic problems. Nanofluidics is the study of the fluid motion in confined structures whose characteristic size is some nanometers, typically 1 − 100 nm (Eijkel & Berg, 2005). Confined fluids in nanoscale geometries exhibit physical behaviors that, in several cases, largely differ from macroscale dynamics. The crucial difference is that in nanoscale systems the usual mathematical description for continuum fluid dynamics often fails to reproduce the correct behavior. Here we deal with simple liquids (Hansen & McDonald, 2006) and the appropriate macroscopic model is given by the incompressible Navier-Stokes equations for mass and momentum conservation


Introduction
In recent years, due to the progress in the fabrication of micro and nanodevices, nanofluidics has become an intense research field.The interest of the scientific community is evident from the large amount of published papers and by the issuing of dedicated journals.In this scenario, a better understanding of the key aspects of fluid motion in nanoscale systems is fundamental and hence computational techniques able to provide a deeper insight in the complex phenomena involved in nanoscale mass transport play a crucial role in nanofluidic research.The aim of this chapter is to introduce the reader to the application of classical all-atom molecular dynamics (MD) to nanofluidic problems.Nanofluidics is the study of the fluid motion in confined structures whose characteristic size is some nanometers, typically 1 − 100 nm (Eijkel & Berg, 2005).Confined fluids in nanoscale geometries exhibit physical behaviors that, in several cases, largely differ from macroscale dynamics.The crucial difference is that in nanoscale systems the usual mathematical description for continuum fluid dynamics often fails to reproduce the correct behavior.Here we deal with simple liquids (Hansen & McDonald, 2006) and the appropriate macroscopic model is given by the incompressible Navier-Stokes equations for mass and momentum conservation ∇•u = 0 (1) where ρ is the constant fluid mass density, u the fluid velocity, p the pressure, ν the kinematic viscosity and f is the force per unit of mass due to external loads.Equations ( 1) and ( 2) are usually completed with the impermeability and the no-slip boundary conditions at solid walls that, taken together, can be written as where D is the solid-liquid interface.The absence of slip at a rigid wall is largely confirmed by direct observations at the macroscale and there are only few well documented cases (see Lauga et al. (2005)) where the use of the no-slip boundary condition (3) at solid walls does not reproduce the correct fluid behavior at macroscale.
Since fluids are composed by molecules it is, in principle, always possible to investigate the flow of a fluid in a nanoconfined system by simulating the motion of each single atom via all-atom molecular dynamics (MD) simulations.This approach, that in general is inefficient and, more often, inapplicable to usual fluid dynamics problems, is crucial for nanofluidics for two main reasons i) it is not based on continuum assumption -that often fails in nanoconfined geometries -and ii) it does not require assumptions on boundary conditions at the interfaces.
In recent years MD proved to be a powerful tool for the analysis of several nanofluidics problems such as, among others, the meniscus and contact line dynamics (De Coninck & Blake, 2008;Gentner et al., 2003), the role of precursor film in wetting (Chibbaro et al., 2008), the thermal exchange properties of carbon nanotube (Chiavazzo & Asinari, 2011) and the interface dynamics of a two immiscible liquid system (Orlandini et al., 2011).Stimulated by the experimental results concerning the flow rate through carbon nanotubes that, for narrow channels (a few nanometers), exceed predictions from the no-slip Poiseuille flow by up to several orders of magnitude (Majumder et al., 2005), several authors used MD simulation to investigate the liquid transport through nanopores.Both Lennard-Jones (Cannon & Hess, 2010;Chinappi et al., 2006) and more realistic models, where water and pore are modeled with the state of the art of classical force fields (Falk et al., 2010;Hummer et al., 2001;Thomas et al., 2010) have been used.Interpretations of the observed flow rate enhancement in terms of viscosity changes in the depletion region close to the wall (Myers, 2011) and change of bulk and interface properties due to confinement (Thomas et al., 2010) have been proposed and, in part, tested via MD techniques.Another nanofluidics topic that has attracted the interest of researcher is the liquid slippage on solid walls, i.e. the presence of a finite velocity at the wall and, hence, the failure of no-slip boundary condition (3).As we will see later in more details slippage is associated to both chemical and geometrical features of the solid surface.MD was largely used to explore the role of surface hydrophobicity (Chinappi & Casciola, 2010;Huang et al., 2008), of surface roughness (Zhang et al., 2011), and of the shear rate influence (Niavarani & Priezjev, 2010;Priezjev et al., 2005) on liquid slippage.
Although the large amount of applications and the increasing computational power of calculation systems, for a large number of nanofluidics phenomena the typical time and length scales accessible to MD simulations are still far from application to realistic systems.This implies that Navier Stokes equation (2) has to be used for the description of nanofluidic phenomena.Two main issues arise when using Navier-Stokes equations for nanofluidic systems, i) is the continuum assumption reasonable at the scale of the system? and ii) which boundary condition has to be applied at the wall?In this chapter we discuss these two issues for the case of liquids in nano confined geometries providing examples of how to employ MD simulations for their analysis.Concerning the former issue, we set up a model system for the estimation of the mass flux through a pore due to an external forcing and compare the results with the hydrodynamic prediction obtained via dimensional analysis of Navier-Stokes equations (2).We show that for simple-liquids the threshold for the validity of the continuum assumption is of the order of five times the molecule dimension and that, below this scale, the hydrodynamic prediction underestimates the flow rate.The second issue is addressed by estimating the slippage for a smooth solid wall using different solid surfaces.In agreement with literature results (Huang et al., 2008) we find a relationship between surface wettability and slippage.The section ends with a discussion on the rough surface case and a comparison with experimental results for a specific system of technological interest consisting in water flowing on an hydrophobic coated surface.
The chapter is structured as follows.A first brief section concerns the continuum assumption and the sub-continuum behavior for a liquid flow through a nanopore.Then we introduce the concept of liquid slippage and we present a simulation set-up for MD characterization of slippage for both smooth and rough surfaces.The final section of the chapter is dedicated to a perspective on near future applications of MD to nanofluidic problems.

Continuum vs single-file motion
As we discussed in the introduction, the main causes for the failure of the prediction of the standard continuum model -eq.( 1), ( 2) and (3) -at the nanoscale are the inappropriate boundary conditions and the failure of the continuum assumption.In each MD simulation aimed at reproducing a nanoscale flow these effects are concomitant and hence it is not easy to isolate the two contributions, and, in particular, it is not possible to clearly identify a threshold for the validity of the continuum assumption and to understand how non-continuum features affect the flow.Here we introduce a MD set-up that overcomes the problem strictly controlling the boundary condition at the solid wall.The comparison of MD results with continuum model prediction is then used to estimate a threshold for the validity of the continuum assumption.As a first step in this program we need to recall a hydrodynamic prediction for the flow rate through pores.
The motion of a liquid in a macroscale system is described by the incompressible Navier-Stokes equation (2).For microscale systems the non-linear convective term (u •∇u) is negligible respect to the viscous term ν∇ 2 u and eq.( 2) reduces to the Stokes equation With a reference length l 0 , to be specified later, a reference diffusive time, t 0 = l 2 0 /ν, speed u 0 = l 0 /t 0 = ν/l 0 and pressure p 0 = ν 2 ρ/l 2 0 , the typical reference force for unit of mass is f 0 = ν 2 /l 3 0 .As a result, the dimensionless parameter where f = |f|, is a natural measure of the external load.With the above choices, the dimensionless form of eq. ( 4) reads where f = f/ f and stars indicate dimensionless units.Since we are interested in molecular flow across a pore, it is natural to identify l 0 with a length that characterizes the pore diameter and that we will indicate as effective diameter d e .Let us use the previous formalism to obtain a hydrodynamic prediction for the mass flux through a pore.The mass flux across a surface of area S is given by with Φ * = S * u * dS * the dimensionless flux.Since Equation ( 6) is linear, u * (and hence Φ * )is proportional to f * and the scaling law

321
Applications of All-Atom Molecular Dynamics to Nanofluidics www.intechopen.comfor the mass flux is found.Equation ( 8) is the well-known power-four law for the mass flux of a viscous fluid in a pipe, that, in the particular case of an infinite pipe of diameter d e with no slip boundary condition at the wall, results in the Hagen-Poiseuille expression φ = πP ′ (d e /2) 4 /(8ν) with P ′ the pressure gradient.It is crucial to note that this dimensional argument is valid only if no other characteristic length scales appear in the problem and, in particular, in the boundary condition.This happens, for instance, in the two cases of no-slip (the velocity at the wall is zero) and no stress (the stress at the wall is zero) boundary conditions but not for partial slip condition where a new characteristic length (the slip length L s , see section 3 below) is present.Equation ( 8) is the continuum model prediction we will compare to MD simulation results.

System set-up
The molecular dynamics set-up presented here is similar to the one presented in Chinappi et al. (2008).For the sake of the clarity we report here the main features, while the interested reader could find the details in the above cited paper.We consider a cylindrical nanopore of height h and circular section of radius r.The nanopore connects two cylindrical reservoirs of radius R (see panel a of Fig. 1).A periodic boundary condition is applied in z-direction, the box length being L z .The liquid molecules are monoatomic and interact via a standard Lennard-Jones (LJ) potential V LJ (r)=4ǫ[(σ/r) 12 − (σ/r) 6 ] truncated at distance r cut = 3.1σ.
In the rest of the section we will use LJ units.In all the simulations the density of the liquid is ρ = 0.83 and the temperature is θ = 0.725.The solid wall is modeled as a continuum that occupies a volume S w and the interaction between the wall and a fluid atom located at r is given by V w (r)= S w n w f (|r − r w |)dr w , where n w is a suitable density having dimension of inverse volume, S w is the wall domain and f (r) is a LJ potential truncated at its minimum.A slab of width h endowed with the cylindrical pore separates the two reservoirs.The aspect ratio of the domain and the wall is kept constant in all simulations, namely, L z /r = 10, R/r = 4, h/r = 1.The origin of reference system is set at the center of the pore with the z-coordinate running along its axis.A sketch of the simulation box is reported in panel a of Fig. 1.Due to the steepness of the wall-liquid potential V w (r) the isosurface V w = k b θ is a natural candidate for the boundary of the volume accessible by liquid atoms.Hence, we define the effective diameter d e as the diameter of the narrow part of the pore delimited by the isosurface V w = k b θ.At this virtual interface, no tangential forces are exerted on molecules.As a consequence of that in the hydrodynamic description this virtual interface corresponds to a free-slip (no stress) impermeable boundary.The flux across the pore is induced by a homogeneous forcing f acting on each liquid atom along the pore axis direction.
The coupling to the heat-bath in the non-equilibrium simulations is achieved via a Berendsen thermostat (Berendsen et al., 1984) and preliminary tests provided confidence on the small sensitiveness of the presented results to changes in thermostat's characteristic time constant.All the simulations were performed with a molecular dynamics code obtained on the basis of the open-source code DL_PROTEIN-2.1 (Melchionna & Cozzini, 2001).
Simulations were performed for different values of the effective pore diameter d e with the number of atoms ranging from N = 435 for d e = 1.83 to 31680 for the largest one at d e = 9.25.Following Zhu et al. (2004), the mass flux is expressed by means of the collective variable n, defined in differential form by its increment in the time interval dt: where dz i is the displacement of the i − th molecule in a time step dt, and the sum runs over all the molecules within the channel at time t, i.e. −h/2 < z i (t) < h/2.Hence each molecule crossing the channel from left/right to right/left is associated to an increase/decrease n → n ± 1.The integer part of n(t + ∆t) − n(t) measures the number of molecules which cross the channel from left-to-right, minus the ones which cross the channel from right-to-left, in a time interval ∆t.In statistically stationary conditions, the flux of molecules is defined as , where the average is taken by sampling the system in time.For a given pore diameter d e simulation were performed at different forcing intensities f in order to verify that the system is in the linear response regime, i.e.Φ∝ f .In the following only linear response results are considered.

Flow through a cylindrical nanopore
In order to test the validity of the continuum prediction eq. ( 8) we plot in panel d of Fig. 1 the quantity Φd −4 e / f as a function of the effective diameter d e .It is apparent that for high d e the curve tends to a constant value as expected from equation ( 8) that predicts a power four scaling law of the flux Φ with the diameter d e .This hydrodynamic behavior sets in at an effective diameter in between 5 and 6 van der Waals radii.This result is coherent with existing literature (see among others Koplik & Banavar (1995) and Bocquet & Charlaix (2009)) indicating that continuum approach for simple liquids is valid on length scales one order of magnitude larger than molecule size.Snapshots of the simulated systems for d e = 1.83 and d e = 9.25 are provided in Fig. 1b, where typical configurations of the molecules are reported for equilibrium (no forcing) simulations.The molecules are represented as van der Waals spheres.Those inside the pore are drawn in actual size, while the radii of those outside the pore have been arbitrarily reduced for clarity.
Decreasing d e the flux is larger than the hydrodynamic prediction and it appears to scale roughly as d 3 e .Observing the snapshot of typical configurations of the molecules for low d e we see that in the extreme case (d e = 1.83 lower panel of Fig. 1b) only a single atom can occupy the pore section.This evidence suggests that an appropriate framework for such a quasi-1D motion is the so called single-file motion, namely a sequential motion of concomitant molecules along a line, with no possibility of overtaking.Being the molecules densely packed within the nanopore, each time a molecule enters the inlet mouth, a molecule is kicked-out from the outlet (see Fig. 1c).The many-body aspects of single-file motion can be described by a single parameter, the hopping-rate (Berezhkovskii & Hummer, 2002) k, that is the inverse of the characteristic time at which molecules hop inside the inlet mouth of the nanopore, with a sufficient energy to move all the molecules inside the pore, so that the last one is expelled.As suggested by Zhu et al. (2004) in describing the flux of water through a carbon nanotube, the relevant parameter is the potential jump associated with the passage of a particle through the pore, ∆μ = fL z .At low forcing intensity, the flux is linear in the potential jump For each panel the image on the left is a projection on zy plane, while the image on the right is the projection on the xy plane.For the sake of clarity particles inside the pores (i.e.0.5h < z < 0.5h) are drawn as spheres of diameter σ while the other particles as spheres of diameter 0.1σ.The image is realized using the VMD software (Humphrey et al., 1996).c) Schematic representation of single-file motion.When a particle enters the pore (e.g.particle 1 of the upper panel), the last particle in the pore exits from the opposite side (particle 5 in lower panel).where k b is the Boltzmann constant.
In stochastic models for single-file transport (Berezhkovskii & Hummer, 2002), the hopping rate k is a phenomenological constant of the model.Roughly speaking, for slight deviations from equilibrium, the hopping rate k should depend on fluid density and temperature and, of course, on the cross section of the pore.Hence, at a given thermodynamic state, it is reasonable to argue that the hopping rate is proportional to the pore cross section, i.e. k ∝ d 2 e .Since, given f, ∆μ scales with the longitudinal dimension of the pore, the above expression indicates that the molecular flux in the single-file regime scales as d 3 e , a result in agreement with the molecular dynamics finding.The presented results could be qualitatively compared with recent MD results on mass flux through carbon nanotubes.Indeed stimulated by experimental studies of water flowing through carbon nanotubes (reporting flow rates exceeding the predictions given by the no-slip Hagen-Poiseuille flow by orders of magnitude (Holt et al., 2006;Majumder et al., 2005;Whitby et al., 2008)) several research groups studied the water flow in nanotubes via MD simulations.Falk et al. (2010) analyzed the water flux inside carbon nanotube (CNT) of different sizes and found that there is a transition in the friction coefficient when the pore diameter is smaller than a couple of nanometers.In particular they found that the water structure is affected by confinement for CNTs of radius below 1.6nm, i.e. ∼ 5 times larger that the size of water molecules, and that, below this regime, single-file motion sets in and the friction coefficient vanishes.In the system analyzed by Falk et al. (2010) both boundary (slippage) and confinement effects are present and hence their results could not be quantitatively compared to our case.However it is apparent that the threshold for the separation between the two regimes is in both cases ∼ 5 times the molecule dimension.A similar threshold was found in Thomas et al. (2010) where the authors show that a continuum model is able to reproduce the MD results if slip length and water viscosity changes with the pore size.In particular, the flow enhancement for narrow pores is interpreted in the continuum model as a decrease of the viscosity and an increase of the slippage at the boundary.In this respect it is interesting to point out that recently in Myers (2011) a continuum model based on reduced viscosity in the depletion region (i.e. the viscosity is not uniform and it decreases close to the wall) and no-slip boundary was proposed.This model correctly predicts the flow enhancement in narrow nanotubes, the enhancement factor being analogous of the one obtained using a finite slip-length.The presence in literature of different theoretical models able to interpret the same phenomena remarks the role of MD as a powerful tool for nanofluidic research since proper MD set-ups allow to isolate the different causes that contribute to the observed behaviour.

Liquid slippage on solid walls
In the continuum framework, liquid slippage at the wall is described in terms of the Navier boundary condition, which, in the case of parallel flow over a non-moving planar wall, reads where v(z) is the velocity profile with z the wall-normal coordinate and v w is the slip velocity, i.e. the value the velocity profile attains at the wall.The parameter L s is called the slip length.
It is geometrically interpreted as the distance below the wall where the extrapolated fluid velocity vanishes, see Fig. Fig. 2. Slippage on a planar wall.Left, no-slip surface: the velocity at the wall is zero.Right: partial slip surface, Navier boundary condition: the velocity at the wall v w is proportional to the velocity gradient dv(z)/dz with z directed as the wall normal toward the liquid (eq.( 11)).
The constant of proportionality is called the slip length (L s ).
typically classified into two broad classes: intrinsic slip, sometimes called molecular slip, and apparent slip.In the first case (intrinsic slip) a non-vanishing slip velocity at the smooth wall results from the first few layers of liquid molecules sliding on the solid surface.In the second case the slip appears at scales intermediate between the characteristic size of the system and the molecular scale.A typical case of apparent slippage is the presence of gas nano-bubbles trapped in the surface roughness elements (center panel of Fig. 5 showing the so called Cassie state), as often happens both for regularly patterned (Gogte et al., 2005) and randomly rough surfaces (Govardhan et al., 2009).The capability of a surface to trap gas bubbles is often associated to the so called superhydrophobic condition where the surfaces are characterized by a low wettability (large contact angle of a water sessile drop together with a low contact angle hysteresis) and slippage (Rothstein, 2010).
In the next section we present some results concerning slippage on smooth walls and discuss the relationship between slippage and wettability.The last section is dedicated to an example of rough surface with nanoscopic defects.

Intrinsic slippage on smooth surfaces
On smooth walls the mechanism that is responsible for liquid slippage is the so called intrinsic or molecular slippage where the first few layers of liquid molecules slide on the solid surface.Intrinsic slippage has been largely studied via MD simulations.The general picture emerging is that wettability and slippage are deeply related, in particular the larger the contact angle θ the larger the slip length L s .In this respect the data presented by Huang et al. (2008) support the existence of a quasi-universal relationship between contact angle θ and slip length L s .However a recent research (Ho et al., 2011) reports MD simulation results that clearly indicate that water slippage could occur also for hydrophilic surfaces suggesting that the connection between θ and L s is purely coincidental.Here we use two simple simulation set-ups to measure the contact angle and the slip length on solid surface with different degree of hydrophobicity and discuss our results in the framework of the ongoing debate on the role of wettability in liquid slippage.
Table 1.Summary of the contact angle measurement simulations.The second column reports the value of the parameter c SL that rules the attractive part of the modified Lennard-Jones potential (12) for the solid-liquid interaction (c SL = 1 for standard LJ potential, c SL = 0 for repulsive potential).

Contact angle measurement
The system we consider is formed by a solid Lennard-Jones (LJ) wall and by a LJ liquid droplet (Fig. 4a).Each atom interacts with the others via a modified LJ potential with c ij the parameter used to change the affinity between atoms, indeed c ij = 0 corresponds to a completely repulsive interaction while for c ij = 1 the usual attractive tail of LJ potential is recovered.For liquid-liquid interaction we used standard LJ potential, i.e. ǫ LL = 1, σ LL = 1 and c LL = 1.The solid is more self attracting than the liquid, ǫ SS = 10, σ SS = 1 and c SS = 1 moreover solid atoms are constrained to a face centered cubic (FCC) lattice by a harmonic spring.Concerning solid-liquid interaction we have ǫ SL = 1, σ SL = 1 while c SL is varied from 0.1 to 1.In the initial configuration the solid atoms forms a FCC slab of dimension L x , L y and h in x,y, and z direction respectively, while the liquid atoms are arranged as a spherical cut, the center of the sphere being at a distance b from the last layer of solid atoms.Periodic boundary conditions are applied in the three directions being L x , L y and L z the dimensions of the periodic box.Both solid and liquid atoms' initial positions are on a FCC lattice of density ρ w = ρ l = 0.83.Initial velocities are assigned to give to the system an initial kinetic energy corresponding to a temperature T = 0.75.During the equilibration phase a thermostat is applied to both wall and liquid atoms.The system temperature T = 0.75 (smaller than LJ fluid critical temperature) and the volume available to fluid atoms correspond to a two phase liquid-vapor system on the LJ phase diagram (Hansen & McDonald, 2006).During equilibration the initial FCC lattice structure used to assign the initial position to the liquid atoms disappears and some atoms leave the droplet surface and enter in the vapor phase.Moreover the droplet rearranges until reaching a steady state where the contact angle θ does not change.The system is considered at equilibrium when the time evolution of the number of fluid particles in vapor phase and the contact angle do not appreciably change.At this point the thermostat is turned off and a NVE equilibrium simulation is run.In all the cases analyzed we do not noticed further changes of the droplet contact angle during the NVE run.and e of Fig. 3) to a weak hydrophobicity for c SL = 0.7 (panels b,d and f of Fig. 3).We observe (data not shown) that the final state is independent from the initial position of the droplet and in particular from the value of b that rules the contact angle for the initial configuration.However the time needed to reach equilibrium dramatically increase if the initial guess of contact angle is far from the equilibrium value, this effect is particularly relevant if the initial contact angle is much larger than the equilibrium one.Equilibrium contact angle is estimated as in Werder et al. (2003) and Chinappi & Casciola (2010), i.e. calculating the density profile, defining the droplet surface as the set of points for which ρ is the half of the liquid density inside the droplet and fitting the surface points to a sphere.It is known that at nanometric scale, the contact angle of a drop may significantly differ from its macroscopic value due to line-tension.For instance in the case of a liquid water droplet on a hydrophobically coated surface (Chinappi & Casciola, 2010), the equilibrium contact angle is θ 120 • for a droplet of radius r 34Å and θ 112 • for a droplet of radius r 60Å leading to a macroscopic contact angle θ macro 105 • , obtained fitting the modified Young equation, cos θ = cos θ macro − τ/(r b γ lv ) with τ the line tension at the three-phase line, r b its curvature radius and γ lv the surface tension at the liquid-vapor interface.For the purpose of the present section, since we are interested in the correlation between contact angle θ and slippage and not in the measurement the precise values of θ, we do not consider this systematic correction due to the line tension τ.In Table 1 the value of c SL for the cases considered in this chapter and the measured contact angles θ are reported.It is apparent that θ increases when the attraction between solid and liquid decreases and, in particular, for c SL = 0.1 (case E) the thermal agitation is able to detach the drop from the solid wall.

The Couette flow MD set-up and the slip length measurement
A usual way to measure intrinsic slippage by MD simulation is to induce a shear in the fluid and to estimate the slip length L s from the extrapolation of the bulk velocity profile.To this purpose the Couette flow is a natural choice since the shear is homogeneous and, hence, the bulk velocity profile in stationary state is linear.In order to induce a Couette flow we prepared a system formed by a channel where a liquid is confined by two solid walls.Periodic boundary conditions are implemented in wall parallel directions (x and y) being L x and L y the box dimensions.As in the case of droplet simulations the lower wall is constrained in a FCC lattice by harmonic potential.The upper wall is constrained only in z direction and a constant force is applied to each wall atom in the x direction resulting, after a transient, in a stationary Couette flow.The atoms of the lower wall are coupled to a thermostat in order to dissipate the heat produced by viscous drag in the liquid.

Results, effect of liquid-solid interaction
In panel b of Fig. 4 the intrinsic slip length L s is plotted as a function of the equilibrium contact angle θ.It is apparent that the for the hydrophilic surface there is no slippage.Moreover the slippage is found to increase with the contact angle.The figure reports also the data from a previous study on water slippage on hydrophobic coatings (Chinappi & Casciola, 2010).These results qualitatively confirm the picture of Huang et al. (2008) that propose a quasi universal relationship between the contact angle θ and intrinsic slip length L s .In particular we do not observe slippage for hydrophilic walls.As we pointed out in section 3.1 a recent research (Ho et al., 2011) reports a positive slip length also for hydrophilic surfaces.In that paper the slippage is observed for hydrophilic surface only in the case the wall lattice spacing is  significantly smaller that liquid molecule size.Here, as in Huang et al. (2008), the dimension of the liquid molecules is similar to solid lattice spacing, hence it is not surprising that our results confirm the picture that associates hydrophobicity and slippage.On the other hand this fact suggests a natural way to further investigate the issue.Indeed it is easy to systematically vary the solid lattice dimension and repeat both contact angle and slippage measurements.Such analysis, performed with a minimal model such as the LJ system, could reveal if the slippage on hydrophilic surfaces is a general phenomenon to be ascribed only to the ratio between solid and liquid molecule sizes or is due to specific choices for the force field implemented in Ho et al. (2011).

Slippage on rough surfaces: the example of OTS coatings
In order to address the role of surface nanoscale defects on slippage, we consider the case of an Octedecyltrichlorosilane (OTS) coated surface.The OTS molecule is formed by a linear alkyl of 17 carbon atoms with a methyl group on one end and a SiCl 3 (silane) group on the other end.OTS molecules are able to spontaneously form layers (sometimes a monolayer) where the molecules are assembled in hexagonal cells with the silane group attached to the solid surface and the methyl group exposed toward the fluid.Due to their ability to form compact layers exposing the methyl (hydrophobic) group, OTS coatings represent a promising technology for surface functionalization.Several groups quantified the slip length L s of liquid water on smooth OTS coated surfaces (Bouzigues et al., 2008;Cottin-Bizonne et al., 2008 Fig. 4. a) Sketch of the MD set-up for slip length measurement.The lower wall is constrained and thermalized while a constant forcing is applied to the upper wall in wall parallel direction, moreover the z position of the upper wall is constrained.Periodic boundary conditions are applied in x and y. b) Contact angle θ vs slip length L s for the simulations described in the section 3.1 (filled squares).The empty circle corresponds to a result of a previous study on wettability and slippage of an Octadecyltichlorosilane (OTS) coated surface Chinappi & Casciola (2010).In order to map this result in LJ units we use the value 3.3Å as van der Waals radius of the water molecule.

2010).
Here our aim is to assess the effect of nanoscale roughness of the coating on the water slippage on a OTS layer.

Defected OTS coating: MD set-up
The typical defect considered here consists of a hole of diameter D. In the greater part of the simulations the defect is obtained removing the OTS molecules in a circle of diameter D hence exposing the LJ (hydrophilic) uncoated surface.For the largest system considered here, for reason that will be clear later, the defect is obtained in a slightly different way; using alkyl molecules of different lengths, namely 11 carbon atoms alkyl chain for the defect and 29 carbon atoms alkyl chain for other molecules.The two different length chains are represented as dark and light blue molecules in the left panel of Fig. 6.Concerning the water molecules we used the TIP3P model (Jorgensen et al., 1983).This model fails to reproduce some of the water properties, in particular the viscosity (one third of the actual water viscosity) and the surface tension (slightly smaller).More accurate models exist such as the TIP4P/2005 (Abascal & Vega, 2005; Alejandre & Chapela, 2010) but they are more computationally demanding and, for the specific case of slippage on smooth OTS surfaces, they did not lead to results qualitatively different from the TIP3P model (Chinappi et al., 2011).A summary of the simulations, performed with the NAMD software (Phillips et al., 2005), is reported in Table 2.The box dimensions L x and L y are reported in columns 2 and 3 while column 4 reports the diameter of the defect.Equilibration was performed at 300 K and 1 atm.Systems were equilibrated following the same procedure described in Chinappi & Casciola (2010) for water slippage on OTS ideal (not defected) coatings and here briefly reported for Liquid Wall Ls

Gas
Cassie state Wenzel state Wall Liquid v(z) Fig. 5. Wenzel and Cassie states.In Wenzel state (left) the liquid fills all the asperities and wets the whole surface, while in the Cassie state (center) gas (or vapor) pockets are trapped in the surface grooves resulting in a patterned interface (liquid-solid zones alternated with liquid-gas zones).The right panel illustrates the concept of apparent slip on a Cassie state (superhydrophobic) surface.The negligible stresses at the liquid-vapor interface results in a bulk velocity profile that could be reproduced by a continuum model with a Navier boundary condition with a uniform slip length L s on an effective flat surface.
the sake of clarity.The wall is modeled as a Lennard-Jones (LJ) solid with the parameters σ and ǫ selected to have a melting point at 1 atm well above the simulation temperature.The solid-liquid interface is parallel to the xy plane and corresponds to a 111 plane of the LJ FCC structure.The alkyl chain head group binds the wall and it is treated in hybrid manner, as a LJ atoms of the wall as concerning no-bonded interactions and as a carbon atom of the alkyl chains for the interactions with bonded atoms along the same chain.The equilibration phase was performed in a triply-periodic box where, by periodicity, the coated wall forms a unique bunch of solid with the upper wall.During equilibration a Langevin piston is applied in the wall-normal direction (z) in order to relax the system to the desired thermodynamic state.After equilibration the two walls were separated by inserting a void region (larger than the LJ cut-off radius of 12Å) before starting the Couette simulation.Concerning the coated wall, the atoms of its lower plane were keep fixed, while the other LJ atoms were coupled to a Langevin thermostat.Concerning the uncoated wall, the atoms of its upper plane were constrained in the wall-normal direction by a harmonic spring and, as in the case of smooth walls described in the previous section, a constant force parallel to the solid-liquid interface is applied to all the upper wall atoms.

Results
Water molecules did not fill the hole during the equilibration phase for all the cases we simulated with the exception of the largest hole we considered (case E in Table 2) where in the first stage of the equilibration process a great number of water molecules enter the defects and get trapped at the LJ hydrophilic surface resulting in a complete wetting (Wenzel) state (panel a of Fig. 5).In order to avoid this effect we introduced the slightly different system we discussed in the previous section where the defect is obtained by using alkyl chain of different lengths resulting in a hole where also the bottom surface is hydrophobic (see the left panel of Fig. 6).This system shows a very stable Cassie state.For all the simulations the slip length L s is measured from the velocity profile (see right panel of Fig. 6).For the five cases considered L s is reported in the sixth column of the Table 2.It is apparent that L s increases with the hole diameter and with the system size.Moreover L s is larger with respect to the smooth case (L s 5 − 7Å (Chinappi & Casciola, 2010)) as expected from the vanishing friction at the liquid-vapor interface.
Following the definition introduced in section 3 the observed slippage has to be classified as apparent slippage.In this context it is interesting to compare our MD data with continuum model results that are available in literature.The simplest way to study the apparent slippage with a continuum model is to consider a patterned surfaces where the interface is composed by solid areas where the local slippage either vanishes or, alternatively, conforms to the small intrinsic slip at solid-liquid interface (L in ) and by gaseous areas (corresponding to the nanobubbles) where no shear stress is acting on the liquid phase (Ng & Wang, 2010;Ybert et al., 2007).Defining the solid fraction Φ s as the ratio between the area of the liquid-solid interface and the projected area, i.e., for the present case Φ s = 1 − πD 2 /(4L x L y ) reported in column 5 of Table 2), Ng & Wang (2010) found that, if no intrinsic slip is assumed (L in = 0at the solid-liquid interface), a continuum model based on Stokes equation leads to the relation where A is the cell side length of a squared periodic lattice.The suffix "0" in expression ( 13) indicates that no-slip (L in = 0) is assumed at the liquid-solid portion of the interface.For a partially slipping solid surface (L in > 0) embedding the free-slip hole, the apparent slip length increases in proportion to the intrinsic slip length L in and in inverse proportion to the solid fraction Φ s (Ybert et al., 2007), with order one proportionality constant.Eqs. ( 13) and ( 14) allow to predict the apparent slip length for the case considered in our MD simulations, the prediction is reported in column 6 of Table 2 for three values of the intrinsic slip length, namely L in = 0, 5, 10Å.The comparison between MD and continuum model indicates that the continuum description of the apparent slip on patterned surface is valid also at the nanoscale.Moreover the value of intrinsic slippage for L in that have to be used in equation ( 14) to obtain a quantitative agreement is in between 5 and 10Å that is consistent with the intrinsic slip measured for smooth OTS coated surfaces (Chinappi & Casciola, 2010).A similar quantitative agreement was evidenced by the comparison of MD results for LJ fluid slippage on superhydrophobic surface by Cottin-Bizonne et al. (2003) and lattice Boltzmann simulation by Benzi et al. (2006).
As in the case of simple liquid analyzed in section 2 the assessment of the capability of a continuum model to reproduce (or not) a nanoscale fluid dynamics system is a precious contribution that MD is able to provide to nanofluidic research.
Moreover the presented results concerning apparent slippage on defected OTS coatings and, in particular, the stability of the superhydrophobic (Cassie) state we observed for all the considered systems, suggest a possible explanation for an interesting issue pointed out by a careful analysis of both experimental and numerical results on water slippage on hydrophobic surfaces presented by Bocquet & Charlaix (2009).The issue is the following: while for where A =(L x L y ) 0.5 is used as characteristic length.The three values correspond to intrinsic slip on solid-liquid interface L in = 0 (no-slip at solid-liquid interface), and L in = 5, 10 Å.The symbol −− for case C with no slip condition on solid surface is due to the fact that, as pointed out by Ng & Wang (2010) expression ( 13) provides a good fit of their numerical data only in the range of Φ s ∈ (0.22, 0.75).
smooth surfaces an increase in slip length with hydrophobicity was found for both MD and experiments, for a given contact angle θ the value of the experimental estimated L s is larger than MD results of about one order of magnitude.In particular in the case of the OTS coatings the most credited experimental data (obtained with different techniques (Bouzigues et al., 2008;Cottin-Bizonne et al., 2008;Li & Yoda, 2010)) indicates a slip length of ∼ 20nm while the MD value is in the range 0.5 − 1.5nm.Since the MD simulations were performed by different authors (Chinappi & Casciola, 2010;Chinappi et al., 2011;Huang et al., 2008;Sendner et al., 2009) using different force fields, computational codes and numerical set-ups, this discrepancy could hardly be ascribed to modeling inaccuracies affecting the MD simulations.
Hence we turn our attention on the putative presence of a non-negligible amount of wall roughness in the smooth surface analyzed in the cited experiments.Indeed several studies on the structure of OTS coatings suggest that perfectly smooth coatings as those analyzed in ideal MD simulations do not exist in practice.Cottin-Bizonne et al. (2008) and Joseph & Tabeling (2005) indicate the peak to peak distance between asperities in their OTS coatings in a few nanometers.Additionally, spectroscopic variable angle ellipsometry, neutron reflection and atomic force microscopy by Styrkas et al. (1999) suggest that OTS films often consist of nearby multilayer domains of different thickness, ranging from 1 to 3 and even 4 times the monolayers thickness.Our conjecture is that apparent slip effects occurring on nanopatterned surfaces can, in principle, explain the difference observed between MD simulations and experiments.By being able to trap nanobubbles, such defects locally induce slippage and increase the apparent slip length.As we show, the resulting apparent L s can be estimated with continuum model and hence, a value for the typical defect size needed to obtain the observed slippage can be estimated from equation ( 13) and ( 14).In order to validate or reject the proposed scenario, the following issues need to be addressed with care: i) Presence of suitable nanoscale surface defects able to explain the experimentally measured L s .ii) Stability of the superhydrophobic (Cassie) state for those nanoscale defects.The stable Cassie state we observed for all the discussed MD simulations provides a first step in the assessment of the second issue, however further analysis are needed.In particular the defect size needed to reconcile MD simulation with experiment is larger that the one analyzed here.Moreover another crucial point to be clarified in the future is the robustness of the picture to changes in the defects geometry and the effect of meniscus curvature.All these issues could be analyzed with appropriate MD simulations.

Conclusion and perspectives
In this chapter we introduced the reader to some applications of classical all-atom molecular dynamics to nanofluidic problems and in particular to two main crucial issues i) the validity continuum assumption at nanoscale for simple liquids and the sub-continuum behavior and ii) the role of hydrophobicity and surface roughness on the liquid slippage at solid walls.For both problems dedicated MD set-ups have been prepared and some results were presented and compared with existing literature.Concerning the former issue we showed, in agreement with previous findings, that continuum model is appropriate for simple liquids when the characteristic size of the system is 5 − 6 times larger than the liquid molecule dimension.Below this threshold the molecular flow across a pore appears to be best described by a single-file model, that predicts a more efficient power-three scaling law for the mass flux.This result, achieved with a minimal model (Lennard-Jones liquid) provides a further ingredient in the interpretation of the high flow rate measured by recent experiment for water flux through carbon nanotubes.Concerning the second issue we discussed an ongoing debate on the connection of hydrophobicity with slippage.Our results suggest that the picture proposed by Huang et al. (2008) is valid when the sizes of liquid and solid molecules are similar however further investigation are needed to assess if the slippage for hydrophilic surface recently observed by Ho et al. (2011) is a general feature that emerges when solid molecules are significantly smaller that liquid ones or it is associated to the specific choices the authors did for liquid and solid molecules.Concerning slippage we also presented results for the role of surface roughness discussed here for a specific system of technological relevance, namely, the OTS coatings, providing a further example of the usage of all-atom MD for nanofluidics research.MD not only allowed to measure the apparent slippage for the specific shape of the defect but, more interestingly, was here used to propose a different interpretation of the experimentally measured slippage for hydrophobic coatings on smooth surface.Future MD studies will provide (or not) further elements to the proposed conjecture eventually stimulating experimental validations.
Beside the results for the specific problems presented here, a secondary aim of the chapter was to use examples taken from ongoing research to give a picture of the capabilities (and the limits) of all-atom Molecular Dynamics simulation for nanofluidic applications.In general the use of MD for nanofluidics problems could divided in two broad classed.From one hand one could try to use the state of the art of the force fields to make a quantitative predictions for a specific issue.Examples are the water flow through carbon nanotubes cited in section 2 and the water slippage on a given surface (such as the OTS coating discussed in section 3.2).
On the other hand MD could be used to perform in silico experiments on simple systems employing a minimal model for the interatomic interaction (e.g.Lennard-Jones liquids) aimed at isolating and analyzing new phenomena or to address specific questions such as the validity of continuum assumption 2, the role of hydrophobicity in wetting (section 3.1 presented here and many other issues discussed in the literature).Both usages will be crucial for nanofluidics research and the increase of computational power in the next years will, for sure, allow the researcher to tackle problems on length and time scales that are out of the range that can be presently explored via MD, such as macromolecules behavior in liquid flows (currently addressed via coarse grained model) and large multiphase systems.
Fig. 1. a) Sketch of the system geometry.The isosurface V w (r)=k b θ (dot-dashed line) represents a natural choice for the boundary of the volume available to the liquid atoms, hence, the effective diameter d e is defined as the diameter of the narrow part of the pore delimited by V w = k b θ.At this virtual interface, no tangential forces are exerted on molecules.b) Snapshots of equilibrium simulations having two different pore sizes: namely, lower panel d e = 1.83, upper panel d e = 9.25.For each panel the image on the left is a projection on zy plane, while the image on the right is the projection on the xy plane.For the sake of clarity particles inside the pores (i.e.0.5h < z < 0.5h) are drawn as spheres of diameter σ while the other particles as spheres of diameter 0.1σ.The image is realized using the VMD software(Humphrey et al., 1996).c) Schematic representation of single-file motion.When a particle enters the pore (e.g.particle 1 of the upper panel), the last particle in the pore exits from the opposite side (particle 5 in lower panel).d) Particle flux Φ divided by the forcing intensity f and the hydrodynamic scaling factor d 4 e .The dashed line represents the single-file scaling Φ∝d 3 e .
Fig. 1. a) Sketch of the system geometry.The isosurface V w (r)=k b θ (dot-dashed line) represents a natural choice for the boundary of the volume available to the liquid atoms, hence, the effective diameter d e is defined as the diameter of the narrow part of the pore delimited by V w = k b θ.At this virtual interface, no tangential forces are exerted on molecules.b) Snapshots of equilibrium simulations having two different pore sizes: namely, lower panel d e = 1.83, upper panel d e = 9.25.For each panel the image on the left is a projection on zy plane, while the image on the right is the projection on the xy plane.For the sake of clarity particles inside the pores (i.e.0.5h < z < 0.5h) are drawn as spheres of diameter σ while the other particles as spheres of diameter 0.1σ.The image is realized using the VMD software(Humphrey et al., 1996).c) Schematic representation of single-file motion.When a particle enters the pore (e.g.particle 1 of the upper panel), the last particle in the pore exits from the opposite side (particle 5 in lower panel).d) Particle flux Φ divided by the forcing intensity f and the hydrodynamic scaling factor d 4 e .The dashed line represents the single-file scaling Φ∝d 3 e .
2. Wall slippage was observed for both smooth and patterned surfaces by several groups, for a review see, among others,Lauga et al. (2005).Slippage is

Fig. 3
reports snapshots of the equilibration phase of a droplet from initial configuration to equilibrium state for the cases c SL = 0.3 and c SL = 0.7, starting from an initial configuration where the center of the sphere (liquid phase) is located at a distance b = 2 from the surface.It is apparent that the equilibrium contact angle is strongly affected by the liquid-solid attraction parameter c SL resulting in an almost complete dewetting condition for c SL = 0.3 (panels a,c 327 Applications of All-Atom Molecular Dynamics to Nanofluidics www.intechopen.com

Fig. 3 .
Fig. 3. Time evolution of a droplet from initial condition to steady state for a highly hydrophobic casec SL = 0.3 (panels a,c,e) -and a slightly hydrophobic one, c SL = 0.7 (panels b,d,f).

Fig. 6 .
Fig. 6.Left.Snapshot of the defected OTS coated surface for simulation E of Table 2.The LJ (hydrophilic, in green) surface is coated by alkyl chains of different lengths, 11 carbon atoms (dark blue) in the circular defect and 29 carbon atoms chains otherwise.Right.Velocity profile obtained for the defect in the left panel, the observed slip length is ∼ 25Å.

Table 2 .
Summary of the molecular dynamics simulations of the Couette flow on defected OTS-SAM coatings.Columns 2 and 3 report the periodic cell dimensions, x and L y respectively, column 4 the diameter D of the circular defect and column 5 the solid fraction Φ s = 1 − πD 2 /(4L x L y ).The apparent slip length L s obtained from MD simulations is in the 6 th column.The last column provides the value obtained combining expression 13 and 14