Theory and simulations of water flow through carbon ... - ENS-phys

Home

Search

Collections

Journals

About

Contact us

My IOPscience

Theory and simulations of water flow through carbon nanotubes: prospects and pitfalls

This article has been downloaded from IOPscience. Please scroll down to see the full text article. 2011 J. Phys.: Condens. Matter 23 184110 (http://iopscience.iop.org/0953-8984/23/18/184110) View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 82.225.22.112 The article was downloaded on 29/04/2011 at 20:04

Please note that terms and conditions apply.

IOP PUBLISHING

JOURNAL OF PHYSICS: CONDENSED MATTER

J. Phys.: Condens. Matter 23 (2011) 184110 (10pp)

doi:10.1088/0953-8984/23/18/184110

Theory and simulations of water flow through carbon nanotubes: prospects and pitfalls Douwe Jan Bonthuis1 , Klaus F Rinne1, Kerstin Falk2 , C Nadir Kaplan1,3 , Dominik Horinek1 , A Nihat Berker1,4, Lydéric Bocquet1,2 and Roland R Netz1 1 2 3 4

Physik Department, Technische Universität München, Garching 85748, Germany Laboratoire PMCN, Université Lyon 1, CNRS, UMR 5586, 69622 Villeurbanne, France Martin Fisher School of Physics, Brandeis University, Waltham, MA 02454, USA ¨ Sabancı Universitesi, Orhanlı-Tuzla, Istanbul 34956, Turkey

Received 1 June 2010, in final form 31 August 2010 Published 20 April 2011 Online at stacks.iop.org/JPhysCM/23/184110 Abstract We study water flow through carbon nanotubes using continuum theory and molecular dynamics simulations. The large slip length in carbon nanotubes greatly enhances the pumping and electrokinetic energy conversion efficiency. In the absence of mobile charges, however, the electro-osmotic flow vanishes. Uncharged nanotubes filled with pure water can therefore not be used as electric field-driven pumps, contrary to some recently ventured ideas. This is in agreement with results from a generalized hydrodynamic theory that includes the angular momentum of rotating dipolar molecules. The electro-osmotic flow observed in simulations of such carbon nanotubes is caused by an imprudent implementation of the Lennard-Jones cutoff. We also discuss the influence of other simulation parameters on the spurious electro-osmotic flow. (Some figures in this article are in colour only in the electronic version)

1. Introduction

acquire a better understanding of proteins like aquaporin and other membrane channels, the mechanisms behind single molecular flows need elucidation [10]. Because of the large surface to volume ratio, describing the fluidic properties on this scale involves a precise determination of the electrostatic and hydrodynamic boundary conditions. A carbon nanotube is a promising candidate for use as a channel in nanofluidic devices [11]. In addition, it may serve as a simple model for a biological channel [12]. However, the fact that nanotubes are uncharged limits their use as active electro-osmotic components. In recent years, many alternative ways of pumping fluids with carbon nanotubes have been suggested [13–15]. For flow through tiny capillaries, the electrostatic boundary conditions strongly influence the permeability to charged species. The electrostatics of a small channel embedded in a medium of very low dielectric constant leads to a large self-energy barrier for ions to enter the confined region [16, 17]. In many biological channels, transport of ions is facilitated by inclusion of fixed charges in the channel

The number of applications of nanofluidic technology has been rising rapidly over the past few years. Inherent to the small size of nanofluidic devices, surface characteristics dominate the bulk flow properties, and as device dimensions shrink even further, interfacial effects become ever more important. A limiting case of small scale fluid flow is found in biological systems, where membrane channels and pumps transport fluids and biological molecules on a single molecular scale [1, 2]. Apart from microscopic applications like micrometer-scale laboratories, arrays of nanofluidic channels can be used in macroscopic pumps, desalination devices and electrokinetic power plants [3–7]. Whereas efficiency is not crucial for microscopic applications, it becomes important when nanofluidic devices are used in parallel arrays that are designed to reach macroscopic sizes. For all applications, ranging from microscopic analysis devices to large-scale power plants, a profound comprehension of the behavior of liquid flowing under strong confinement is required [8, 9]. Likewise, to 0953-8984/11/184110+10$33.00

1

© 2011 IOP Publishing Ltd Printed in the UK & the USA

J. Phys.: Condens. Matter 23 (2011) 184110

D J Bonthuis et al

walls [18], coating walls with dipolar surface groups [19, 20] and screening by salt [21]. Fixed charges outside a carbon nanotube are found to affect the pressure-driven passage of water molecules through the nanotube as well [22]. Even more intriguing is the observation that carbon nanotubes filled with pure water exhibit electro-osmotic flow when either an electric field is applied [23] or point charges are fixed outside the channel [12]. Similarly, electro-osmotic flow is found in uncharged channels filled with solutions of the relatively symmetric salts sodium chloride or potassium chloride [24–26]. The appearance of a non-zero ζ -potential in the absence of free charges was speculated to be related to the electrostatics in the boundary layer: water molecules tend to orient, leading to a strong dipolar field in the first few molecular layers next to the interface. The coupling of the electric field to this dipole density was thought to induce flow via rotation of individual molecules [23, 27]. Although exciting, these results raise concern as fundamental laws of physics appear to be violated. In particular, the effects seem at odds both with Onsager’s reciprocal theorem, since an externally applied pressure drop cannot cause a steady electric current due to the absence of free charges, and with thermodynamics, since the electric field performs no work in the steady state. A growing amount of literature indicates that on small scales the hydrodynamic boundary conditions also deviate from the usual no-slip condition [28, 29]. Experimentally, the viscosity is defined as the proportionality constant between the shear force per unit area and the resulting velocity gradient. Close to a wall, this proportionality constant appears to be different from the bulk value. In the first few molecular layers next to the interface, this is caused by the changing water structure, most notably the dipolar orientation and the oscillating density profile. In addition, for hydrophobic walls there is a density depletion gap between the wall and the fluid [30–32]. Directly at the wall, the velocity gradient can therefore be very different from the bulk in a Couette shearing scenario. On a molecular scale, even ‘true slip’ of the first molecular layer along the wall can be imagined. On a simple level, the assumption that the surface stress is linearly related to the surface velocity via a friction coefficient which equates the surface stress to the viscous shear stress leads to the notion of a slip length [33], as sketched in figure 1(a). A finite slip length greatly enhances flow rates through small channels. For planar surfaces, the slip length exhibits a sensitive dependence on the microscopic properties of the surface, increasing as the contact angle grows and the surface becomes more hydrophobic, and scaling proportional to the depletion length to the fourth power [34]. Carbon nanotubes appear to have very large slip lengths, up to tens of micrometers [35]. For direct investigation of interfacial properties, simulation techniques are particularly valuable. Molecular dynamics simulations have provided detailed insight into the mechanism of surface slip [36], which had been elusive and subject to speculation for a long time. However, due to the complexity of a modern-day molecular dynamics simulation package, utmost care must be taken when interpreting the results. In this paper we focus on the use of carbon nanotubes and rectangular slits as nanofluidic channels. We first calculate

Figure 1. (a) Sketch of a slipping velocity profile next to an interface. The slip length b corresponds to the distance where the linear extrapolation of the velocity profile reaches zero. (b) Illustration of the hypothetical mechanism to induce electro-osmotic flow in a dipolar fluid. The electric field exerts a torque on the oriented dipoles, leading to rotation ω2 . The rotation is coupled to the velocity u 1 via the antisymmetric part of the stress tensor.

the energy conversion efficiency gained by using channels with a slip length of the same order of magnitude as the channel width. Secondly, we formulate the fundamental hydrodynamic equations including the presence of angular momentum associated with spinning dipolar molecules. We solve those equations in the presence of an external electric field for flow through nanofluidic slab-like channels and carbon nanotubes. Finally, we discuss the prospects and pitfalls of molecular dynamics simulations of water flowing through carbon nanochannels with specific focus on various simulation artifacts.

2. Pumping and energy conversion 2.1. Dissipation in pressure-driven flows Conservation of momentum for a fluid flowing with spacedependent and time-dependent velocity u(x, t) is expressed by the Navier–Stokes equation, which in the absence of external forces reads

ρ

du =∇·P dt

with

d ∂ = + u · ∇. dt ∂t

(1)

We assume the fluid to be incompressible, leaving the mass density ρ constant. For fluids of spherically symmetric constituents, the stress tensor P in equation (1) is given by P = − pU + %

with

%i j = η(∇ j u i + ∇i u j ). (2)

The hydrostatic pressure is denoted p , U is the unit tensor and η is the viscosity. The balance equation for the kinetic energy contained in a flowing liquid follows from a scalar product of the Navier–Stokes equation with the velocity, giving [37]

ρ

1 du 2 = ∇ · (P · u) − P : (∇u). 2 dt

(3)

In steady state ∂u/∂t = 0 and equation (3) reduces to 0 = ∇ · (−ρu2 u/2 + P · u) − P : (∇u),

(4)

where the first term on the right-hand side represents the convective energy flow ρ(u · ∇)u2 /2 = ρ∇ · (u2 u)/2, which 2


D J Bonthuis et al

need not be zero in steady state. Integration over a volume V and application of Gauss’s theorem yields

" # ! ρu2 u 0= nˆ · − + P · u da − P : (∇u) dx. 2 ∂V V !

from which it can be seen that increasing the slip length greatly reduces the power required to achieve a flux Q when b becomes of the order of d . A pressure difference across a cylindrical channel like a carbon nanotube produces a Hagen–Poiseuille flow profile. At the wall of the tube of diameter d we use the boundary condition given in equation (8) restricted to the wall at +d/2. Inside a tube of diameter d # l the flow profile is given by

(5)

Using Einstein’s summation convention, we write the integral componentwise and insert equation (2), $ % ! ρu 2i u j 0= nˆ j − − pu j + η(∇ j u i + ∇i u j )u i da 2 ∂V ! − η(∇ j u i + ∇i u j )∇ j u i dx, (6)

ui = −

(p 2 [d + 4bd − 4x 2j ]. 16ηl

(12)

The dissipation for fixed Q equals

V

where we used P · u = − pδi j u i + η(∇ j u i + ∇i u j )u i and P : (∇u) = Pi j ∇ j u i and the incompressibility of the fluid,

Pdiss =

∇i u i = 0. The surface of the volume ∂ V consists of three parts: an inlet, an outlet and a wall. At impenetrable walls nˆ j u j = 0. For flow through a translationally invariant channel also the integrals over the inlet and outlet of the first term in the surface integral vanish. The second term in the surface integral represents the input power Q(p , with the volume flux ! Q= nˆ j u j da, (7)

128 Q 2 ηl . π(d 4 + 8bd 3 )

(13)

In carbon nanotubes, the slip length b can be two orders of magnitude larger than in other hydrophobic channels and much larger than the diameter d [35], making the energy dissipation scale inversely linear with the slip length. The resulting massive reduction of the energy dissipation is a major advantage when using large arrays of parallel nanotubes for filtering or desalination.

)

2.2. Electrokinetic energy conversion efficiency

where ) is the cross section of the channel and (p is the static pressure difference between the inlet and outlet. The remaining term in the surface integral in equation (6) represents the power dissipated by friction at the walls and depends on the extent of slip. The volume integral represents the power dissipated within the fluid because of velocity gradients. In our continuum description, we take the viscosity as a constant and introduce a slip length b as a boundary condition,

− nˆ j b∇ j u i |±d/2 = u i |±d/2 ,

The calculation above shows that a large slip length reduces the energy dissipation in pressure-driven flows. However, a hydrostatic pressure difference is sometimes difficult to generate on small scales. It has long been recognized that electro-osmosis is a very simple and scalable mechanism to use in nano-scale pumps [38]. We consider an electro-osmotic pump consisting of a slit-like or cylindrical channel with a fixed surface charge or surface potential operated with a solution of monovalent ions. Under the influence of an electric field, the ions move along the surface and drag the fluid along, driving a hydrostatic load impedance. The load impedance connected to the channel is the hydrostatic equivalent of an electrical load resistance connected to an electrical power source. In this section we calculate the efficiency of the energy conversion from electrical power to a hydrostatic pressure difference as a function of the slip length. Because of the inhomogeneity of the volume and current flow densities, the standard flow densities and Onsager coefficients are integrated over the cross section ) of the channel [37]. In the linear response regime, the electrical current I and fluid volume flux Q are related to the pressure gradient ∇ p and the voltage gradient ∇ϕ via $ % $ %$ % Q ,11 ,12 −∇ p = . (14) I −∇ϕ ,21 ,22

(8)

with d the height of the channel. For single-directional fluid flow, as appropriate for an incompressible fluid flowing through a translationally invariant channel, the balance equation becomes

Q(p = −η

!

∂V

u 2i da − η b

!

(∇ j u i )2 dx,

(9)

V

where the right-hand side equates to the dissipated power Pdiss . For a given volume flux Q , we can now calculate the dissipation as a function of the slip length. Ignoring edge effects, a pressure-driven fluid flow in a slit-like rectangular channel forms a Poiseuille flow profile,

ui = −

(p 2 [d + 4bd − 4x 2j ], 8ηl

(10)

The hydrostatic output power per unit length of the channel equals Pout = Q · ∇ p , while the electrical input power equals Pin = −I · ∇ϕ . For electrokinetic power generation by pressure-driven flow, input and output power are defined the other way around, but the efficiency is calculated the same way. The electrokinetic energy conversion efficiency is given by the ratio E = Pout /Pin , which can be expressed in terms of the electrokinetic response functions via equation (14).

with u i the fluid velocity lengthways, x j the spatial coordinate along the height of the channel and l the length of the channel. For a given volume flux Q , the dissipation from equation (9) in a channel of width w becomes

Pdiss =

12 Q 2 ηl , w(d 3 + 6d 2 b)

(11) 3


D J Bonthuis et al

In the channel, the fluid velocity and electrostatic potential (apart from the constant voltage gradient driving the pump) depend only on the coordinate perpendicular to the channel wall, throughout the following calculation denoted x j . The velocity, pressure gradient and driving voltage gradient are only non-zero in the direction along the channel wall, denoted x i . Using the Navier–Stokes equation (equation (1)) and the Poisson equation, ∇ 2j ψ = −e(n + − n − )/εε0 with ψ(x j ) the electrostatic potential, e the elementary charge, n ± the number densities of positive and negative monovalent ions and εε0 the dielectric constant, the Onsager coefficients can be expressed as ! ( p) u i (x j ) ,11 = − da, ∇i p ) ! εε0 ,12 = ,21 = − [ψ(d/2) − ψ(x j )] da, η ) (15) !

,22 = e

− e

)

!

Figure 2. Efficiency E of electrokinetic energy conversion as a function of slip length b in (a) a rectangular channel and (b) a cylindrical channel, calculated from equation (16). The electrolyte used is 0.1 mM monovalent ions with the electrophoretic mobilities of KCl. For both geometries, the surface potential is fixed at −240 mV. The corresponding surface charge density equals −13 mC m−2 in the d = 0.5 nm cylindrical channel, −24 mC m−2 in the d = 0.5 nm rectangular channel and ∼−60 mC m−2 for all other channels.

[m + n + (x j ) − m − n − (x j )] da,

)

efficiency with respect to the channel diameter d is nonmonotonic and different for different values of the slip length. For a slip length of ∼2 nm, which is a reasonable number for hydrophobic surfaces [36], the efficiency in the 10 nm high rectangular channel reaches ∼20%. In an extension to the work of [42] we plot the efficiency for negative slip lengths. The measured efficiency of 3% in a 75 nm high channel and 1% in a 490 nm high channel [41] with a surface charge density of −60 mC m−2 , which is the measured surface charge density of silicon oxide [43], is consistent with a slip length of b = −0.3 nm, corresponding to one static molecular layer, in good agreement with predictions by [36] for an OH-covered surface. In figure 2(b), the efficiency is shown for a cylindrical channel. The increase in efficiency is similar to the planar case, but in carbon nanotubes the slip length can be as much as ∼50 µm [35]. In all nanotubes considered, the efficiency for slip lengths of that order is over 90%.

εε0 [ψ(d/2) − ψ(x j )][n + (x j ) − n − (x j )] da, η

where the monovalent ions have electrophoretic mobilities m ± and number densities n ± = n exp(∓eψ/kB T ) with n the bulk density. The electrostatic potential is calculated numerically from the Poisson–Boltzmann equation, which holds sufficiently well on small scales [39]. The surface of a carbon nanotube is uncharged, but a finite surface potential can be achieved for nanotubes by application of an external potential. Surface potentials also arise naturally due to unequal surface adsorption of cations and anions. Note that in an experimental setup the driving electric field will be at least partially screened by electron transport through the nanotube, which will reduce the electric field strength ( p) inside the nanotube. The pressure-driven flow velocity u i is calculated from the Navier–Stokes equation as in section 2.1. The relation between the pressure gradient in the channel and the volume flux depends on the hydrostatic load connected to the channel. The hydrostatic load impedance Z L defines the proportionality between the two, Q = ∇ p/Z L . The efficiency is maximized for a specific value of Z L ,11 , which is the hydrostatic load impedance divided by the hydrostatic channel −1 . The value of the maximum efficiency has a impedance ,11 very simple form in terms of the variable α = ,212 /,11 ,22 , given by [40] α . E= √ (16) 2+2 1−α−α

3. Rotation–vorticity coupling 3.1. Conservation of linear and angular momentum The calculations up to this point are performed for charged systems using standard hydrodynamic theory. However, the simulations by [12, 23] show electro-osmotic flow also for uncharged channels filled with pure water. These results raise concern because according to Onsager’s symmetry the streaming conductivity ,21 equals the electro-osmotic coefficient ,12 . This means that when an electric field induces an electro-osmotic flow, the conjugate pressure gradient field should cause an electric current, which seems impossible in a system without free charges. In addition, the energy supplied by a static, space-invariant electric field is only non-zero for a non-zero current density. Nevertheless, an intuitively appealing mechanism to pump a dipolar fluid involves field induced rotation of the dipoles. The hypothetical mechanism is sketched in figure 1(b): an electric field gives rise to a spin field ω2 , which is coupled to the velocity field u 1 . To analyze the interaction of an external electric field with a dipolar fluid we

Using a no-slip boundary condition, the predicted maximum efficiency is only 7%, while the measured efficiency in a silicon oxide channel of d = 75 nm reaches no more than 3% [41]. To investigate the effect of slip, we use the hydrodynamic boundary condition of equation (8). The efficiency from equation (16) is shown in figure 2 as a function of slip length for different channel dimensions. We use the ion mobilities of KCl and a bulk concentration of 0.1 mM, giving a Debye screening length of 30 nm. Clearly, the efficiency increases drastically with increasing slip length. The variation of the

4


D J Bonthuis et al

now extend the Navier–Stokes equations to include the effect of anisotropic and rotating molecules. The generalized hydrodynamic theory that includes the effect of an electric field on dipolar fluid particles involves not only the velocity field u(x, t), but also the spin field ω(x, t), which measures the local mean angular particle velocity. In this case, the antisymmetric part of the stress tensor, which vanishes for monoatomic fluids, is non-zero in general [44]. Instead of equation (2), we now split the viscous part of the stress tensor % into a symmetric part %s and an antisymmetric part %a , P = − pU + %

% = %s + %a .

and

(17)

The fluid flow is determined by the conservation laws for linear and angular momentum. The angular momentum density is the sum of the vorticity contribution ρx × u and the spin contribution ρ I · ω. The moment of inertia per unit mass I is a tensor that in general depends on the moments of inertia of a single particle and the local particle orientational distribution. Because of the approximate mass isotropy of a water molecule, we approximate I to be diagonal and proportional to I . The balance equation for the spin angular momentum reads [45, 46]

ρI

dω = ρ" + ∇ · C + 2#a dt

Figure 3. (a) Sketch of the cylindrical geometry. (b)–(d) Snapshots of the three different carbon nanotubes used in the simulations: (10, 0) with diameter d = 0.782 nm, (16, 0) with d = 1.25 nm and (16, 16) with d = 2.17 nm, respectively.

where c1 is an integration constant to be fixed by boundary conditions. Combining equations (21) and (22) yields the differential equation

(18)

where #a is the antisymmetric part of the shear stress tensor in axial vector representation, " is the external torque density vector and the spin stress tensor is denoted as C. Including the antisymmetric part of the stress tensor, the Navier–Stokes equation for incompressible fluids in the absence of a body force reads

ρ

du = −∇ p + η(u − ∇ × #a , dt

0 = ρ12 /ν + (∇32 − κ 2 )ω2 + 2c1 ηr /(νx 3 ),

for arbitrary torque density 12 and the spin screening length κ −1 defined by κ 2 = 4ηηr /[ν(η + ηr )]. To solve equation (23) for the situation of a constant and static external electric field E 10 in x 1 direction, we calculate the torque " = µ × E . We model the dipolar ordering of the interfacial water layer by an additional x 3 -dependent electric µ field E 3 (x 3 ) in the x 3 -direction. This field accounts for the interaction between neighboring water molecules and plays the role of the crystal field that is used in ordinary mean-field theory for magnetic systems. The total electric field in the µ system becomes E = eˆ3 E 3 (x 3 ) + eˆ1 E 10 . The linear evolution equation for the polarization density µ = µ(x 3 ) is given by

(19)

with p the hydrostatic pressure and η the shear viscosity, and includes via #a a coupling to the spin field ω(x, t). 3.2. Material equations To leading order the material equations for #a and C read

#a = ηr (∇ × u − 2ω)

and

C = ν∇ω,

1 dµ = (αE − µ) + ω × µ, dt τ

(20)

(21)

µ1 =

for i = 1, 2, 3 and with ∇ 2 the Laplacian in cylindrical coordinates. From equation (21) it follows that the spin fields ω1 and ω3 are zero if the torque components 11 and 13 vanish. By integrating equation (19) once over x 3 we obtain

(η + ηr )∇3 u 1 = 2ηr ω2 + c1 (η + ηr )/x 3 ,

(24)

with τ the relaxation time of the orientation and α the polarizability per unit mass. Because of translational invariance and directional symmetry in the x 2 -direction, the dipole moment in the x 2 -direction vanishes: µ2 = 0. Since also E 2 = 0, it follows that the torque 11 = 13 = 0, and therefore ω1 = ω3 = 0 by virtue of equation (21). Equation (24) can be solved for µ in steady state,

with ηr the vortex viscosity and ν the spin viscosity. For a cylindrical geometry as sketched in figure 3, we have u = eˆ 1 u 1 (x 3 ) and ω = ω(x 3 ), with x 1 the direction along and x 3 the direction perpendicular to the channel wall. The only nonzero component of the vorticity ∇ × u is the x 2 -component. In steady state, equations (18) and (20) reduce to 0 = ρ1i + ν∇32 ωi + 2ηr (δi 2 ∇3 u 1 − 2ωi )

(23)

α E 10 + τ ω2 α E 3µ 1 + τ 2 ω22

µ3 =

and

α E 3µ − τ ω2 α E 10 . 1 + τ 2 ω22 (25)

For the torque density it follows that 2

12 = −

(22) 5

2

τ ω2 α(E 10 + E 3µ ) . 1 + τ 2 ω22

(26)


D J Bonthuis et al

Inserting equation (26) into equation (23) yields 2

0=−

µ2

2c1 ηr ρτ ω2 α(E 10 + E 3 ) . (27) + (∇32 − κ 2 )ω2 + νx 3 ν(1 + τ 2 ω22 )

Strikingly, equation (27) is quadratic in the external field strength E 10 , which means that switching the sign of E 10 leaves the equation invariant. Also, the sign of ω2 is undetermined, meaning that a non-zero value of ω2 would require a spontaneous symmetry breaking, which seems unphysical. In fact, it can be shown that in hydrodynamics, the stable solution corresponds to the solution of minimal dissipation or minimal entropy production, which is obviously the solution corresponding to ω2 = 0 [37]. Therefore, the only physical and also stable solution is the one corresponding to a vanishing spin field, ω2 = 0, and thus the integration constant c1 also vanishes. If ω2 = 0, the vorticity also vanishes according to equation (22), and no flow results. From the Helmholtz– Smoluchowski equation,

ζ =−

ηu i , εε0 E i

Figure 4. (a) Lennard-Jones force for the two different truncation schemes. (b)–(d) The ζ -potential of a (10, 0), (16, 0) and (16, 16) carbon nanotube, respectively, as a function of the Lennard-Jones cutoff length. We use an electric field strength of 0.1 V nm−1 ((b)–(c)) and 0.1 − 0.4 V nm−1 (d), where ζ is calculated from a linear fit to the velocity as a function of the field strength. The nanotubes are connected directly to their periodic image, without reservoir.

(28)

it is directly evident that this means that neutral solutes in a purely dipolar fluid without mobile charges have zero ζ potential. Nevertheless, the coupling between spin and vorticity provided by the antisymmetric part of the stress tensor can be exploited to pump a fluid through a channel by using a suitably chosen rotating electric field [47].

length of 0.142 nm and the carbon atoms are modeled by the GROMOS 96 force field. The nanotubes are equilibrated in a large bath of solvent using a semi-isotropic Berendsen barostat at 1 bar. After equilibration, the water outside the nanotube is removed. We use periodic boundary conditions in all directions and particle mesh Ewald summation for the electrostatics. In the axial direction, the simulation box is set to the same size as the nanotube, so the nanotube is connected directly to its periodic image. In the other two directions, the box size is 5 nm × 5 nm. We show in figures 4(b)–(d) the ζ -potential of the tubes, calculated via equation (28) from the number density per unit length inside the tube and a linear fit to the cumulative flux. Surprisingly, the observed ζ -potential for the simple cutoff depends crucially on the Lennard-Jones cutoff length. Whereas a dependence within numerical accuracy may be expected, a critical dependence like the one observed is unphysical. Identical to [23], we also simulate a (16, 0) carbon nanotube of length 9.8 nm connecting two reservoirs as shown in figure 5(a). We use a Nosé–Hoover thermostat and an anisotropic Parrinello–Rahman barostat. For each scheme and various rc , we simulate for 5 ns and collect the last 3 ns for analysis while for rc = 1.0 nm we extend the simulation time by 40 ns. The cumulative flux is shown as a function of time for both cutoff schemes in figure 5(b), and the corresponding ζ -potentials are shown in figure 5(c). Although the ζ -potential is smaller than in the absence of reservoirs, which has to do with the friction a water molecule encounters when entering or exiting the nanotube, again the results show a striking yet unphysical difference between the cutoff schemes: for the simple cutoff at rc = 1.0 nm we find an average flux of 22 ± 8 ns−1 , comparable to 34.0 ± 5.2 ns−1 from [23], while the shifted cutoff exhibits vanishing flux. In the limit rc → ∞,

4. Molecular dynamics simulations 4.1. Homogeneous electric fields in carbon nanotubes In order to keep a molecular dynamics simulation computationally feasible, several approximations are employed. First of all, the interaction potentials are truncated at a finite distance. For the long-ranged electrostatic Coulomb interaction, the error introduced by truncating the potential is compensated by Ewald summation, and the exact value of the real-space Coulomb interaction cutoff length has no influence on the dynamics. The electrostatic Lennard-Jones interaction, however, is considered short-ranged and no compensation is typically used. To clarify the above-mentioned simulations of pure water exhibiting the surprising electro-osmotic flow [12, 23], we simulate several different nanotubes for different values of the Lennard-Jones cutoff length and two different truncation schemes: simple cutoff, which leaves the force unchanged out to rc beyond which it is set to zero, creating a discontinuity at r = rc [23], and shifted cutoff, which adds a nonlinear function to the force over the whole range, ensuring a smooth decay to zero [48]. Both truncation schemes are shown in figure 4(a). We simulate three nanotubes with chiralities (10, 0), (16, 0) and (16, 16), having diameters of 0.782 nm, 1.25 nm and 2.17 nm respectively, as shown in figures 3(b)–(d). The tubes have lengths of 10, 10 and 5 nm, containing 37, 202 and 358 SPC/E water molecules, respectively. The C–C bond has a 6


D J Bonthuis et al

Figure 6. Cumulative flux through a (10, 0) carbon nanotube with (red, marked with a square, field strength 1 V nm−1 ) and without (blue, marked with a circle) external electric field. The cutoff schemes used are (a) shifted and (b) simple, both at rc = 1.0 nm. The flux is calculated according to the unconditional definition (see section 4.4 in the main text). Including reservoirs, the system size is 3.0 × 2.7 × 13.7 nm3 , containing 1008 SPC/E water molecules. Figure 5. (a) Snapshot of the (16, 0) carbon nanotube connecting two reservoirs. The system size is 4.0 × 3.6 × 11.8 nm3 containing 1057 SPC/E water molecules. (b) Cumulative flux through the tube as a function of time for an electric field strength of 1 V nm−1 for the two different Lennard-Jones cutoff schemes simulated with GROMACS. (c) The ζ -potential of the nanotube as a function of the Lennard-Jones cutoff length, simulated with GROMACS (circles and squares), and simulated with LAMMPS (cross and plus), for both cutoff schemes. For each cutoff length, the ζ -potential is calculated from the derivative of the cumulative flux as shown in (b), the line density inside the tube and equation (28).

Figure 7. (a) Snapshot of the (10, 0) nanotube with fixed negative charges. The charges are located at 1.92, 2.64 and 2.88 nm distance from the right end of the tube and the distance between the charges and the nanotube is δ = 0.093 nm. Charges are −0.5, −0.5 and −1e from left to right, with charges of opposite sign at the bottom of the box for compensation. The system size is 8.8 × 8.4 × 8.8 nm3 , the tube length is 4.8 nm and the system contains 8389 TIP3P water molecules. The neighbor list is updated every ten timesteps of 2 fs each. (b) Dependence of the average axial velocity u 1 inside the nanotube on the cutoff length rc for the two different cutoff schemes.

the spurious difference between the cutoff schemes disappears and the water flux vanishes, in accordance with the generalized hydrodynamic theory presented in the previous sections. We simulate exactly the same system of a (16, 0) carbon nanotube between reservoirs using the alternative simulation package LAMMPS [49]. In LAMMPS, both the simple cutoff and the shifted cutoff as used in GROMACS are implemented. As shown in figure 5, the flux vanishes, even for small rc , regardless of the cutoff scheme. We conclude that electro-osmosis, i.e. the electric field induced steady flow, of pure water in a carbon nanotube is a simulation artifact and related to the implementation of the cutoff scheme in GROMACS [48]. As a simple test of the effect of the electric field, we perform long simulations of a 9.8 nm long (10, 0) nanotube between reservoirs with and without electric field. The results are shown in figure 6 for zero and non-zero external electric field. Most strikingly, the flux is non-zero for the simple cutoff scheme, shown in figure 6(b), regardless of whether an external electric field is applied or not. The same effect is in fact displayed in figure 1 of [23]. Again, in the shifted scheme this unphysical effect is absent. Since it is self-evident that in the absence of an electric field no flux should result, this proves our point most forcefully that something is fundamentally wrong with the simple cutoff scheme as implemented in GROMACS, leading to erroneous coupling between orientation and flux. The residual flux fluctuations for the shifted cutoff scheme shown in figure 6(a) suggest that the simulations in the (10, 0) nanotube have not converged on the timescale of 40 ns, which simply makes the (10, 0) nanotube unsuited for studying equilibrium properties in general.

4.2. Point charges at carbon nanotubes We now turn to water-filled carbon nanotubes in the presence of point charges. For a set of fixed point charges the electric field is spatially inhomogeneous, but constant in time. We simulate the system shown in figure 7(a) for 5 ns, collect the last 3 ns for analysis and calculate the flux as a function of the cutoff length for the two different cutoff schemes. We use a Berendsen thermostat and update the neighbor list every ten timesteps. We define the three negative charges next to the nanotube to be in one charge group and the three compensating positive charges in another (see below for further explanation). The differences between our system and the system used in [12] are that the charges next to our tube are negative, external charges and tube atoms are frozen, our tube is longer and our reservoir is larger. From the flux and the average number of particles per unit length we calculate the average velocity inside the tube, shown in figure 7(b) as a function of the cutoff length for the two different truncation schemes. Clearly, the velocity depends on the cutoff length like in the case of a homogeneous electric field and simple cutoff. However, the velocity resulting from an inhomogeneous electric field also shows a cutoff length 7


D J Bonthuis et al

Figure 8. Dependence of the ζ -potential of a hydrophobic diamond (contact angle of 111◦ ) on (a) the number of steps between neighbor list updates (nstlist) using a timestep of 1 fs and (b) the integration time step using nstlist = 50. The electric field strength is 0.4 V nm−1 and the system size is 3.2 nm × 3.2 nm × 4.7 nm containing 928 SPC/E water molecules and one diamond. We use a simple Lennard-Jones cutoff at 0.8 nm.

Figure 9. Spurious ζ -potential of (a) a diamond surface as a function of contact angle and (b) an artificial hydrophobic surface (θ = 140◦ ) with 1 M NaCl as a function of the molecular water density N/V . The dashed lines are drawn as a guide to the eye. An electric field strength of 0.4 V nm−1 , timestep of 1 fs, neighbor list update frequency of 1/50 timesteps and simple Lennard-Jones cutoff at rc = 0.8 nm (a) and rc = 1.0 nm (b) are used. The system size in (a) is identical to the one used in figure 8. The system size in (b) is 3.9 nm × 3.9 nm × 16.4 nm, containing 2052 SPC/E water, 40 Na+ and 40 Cl− molecules confined between two slabs.

dependence for a shifted cutoff scheme. This shows that using a shifted cutoff scheme does not necessarily prevent the flow. Indeed, since there is nothing unphysical about the simple cutoff per se, we have no reason to believe that a shifted cutoff should perform better in any case; it is the numerical implementation of the cutoff scheme in GROMACS that produces questionable results. Following a recent publication [50], we also investigate the effect of the use of charge groups. For a charge group, Coulomb interactions are calculated for all individual charges in the group, but the position of each individual charge is approximated by the geometrical mean of the constituent particles [51]. In the simulations by [12] and the simulations shown in figure 7 the fixed charges close to the nanotube are grouped within a single charge group. The artificial pumping disappears for both cutoff schemes even for a short cutoff length of 0.7 nm if we do not use a charge group for the fixed external charges. This confirms the findings of [50]. Finally, using the alternative simulation package NAMD2, the directional flux vanishes altogether [52].

diamond is shown as a function of the integration time step for rc = 0.8 nm. Time steps larger than 2 fs are not possible because of the fast dynamics of the hydrogen atoms. The flow dynamics do not depend significantly on the integration time step, showing that the problem is not related to integration accuracy. Like any flow, the magnitude of the spurious flow depends strongly on the particularities of the system, most notably the contact angle of the substrate and the pressure in the system. In figure 9(a) we show the ζ -potential of the diamond surface in contact with pure water. The contact angle of the diamond is modified by adjusting the Lennard-Jones interaction strength 5 between the carbon atoms and the water molecules. For each simulation system, the contact angle is determined from the virial according to the method employed by [53]. All simulations are done using a simple cutoff truncation scheme and a truncation length of rc = 0.8 nm. Using shifted truncation, the electro-osmotic flow is zero for any contact angle. We show in figure 9(b) the ζ -potential of an artificial single face-centered cubic lattice (Lennard-Jones parameters 5 = 0.686 kJ mol−1 , σ = 0.337 nm, θ = 140◦ ) in a solution of 1 M NaCl in water as a function of the average number density of water N/V in the volume between the plates. Calculating the pressure from the force on the center of mass of each plate divided by the plate area gives a pressure of 1 bar at an average density of 29 nm−3 . We use a simple cutoff scheme for the Lennard-Jones interaction with a truncation length of 1.0 nm. The apparent ζ -potential decreases with increasing density. We note that the dependence of the ζ -potential on the surface contact angle and the water density shown in figure 9 is spurious, in the sense that the electro-osmotic flow should vanish altogether in this case. This is also true for the salt solution, since the cations and anions of NaCl show almost the same adsorption affinities on a hydrophobic surface. Figure 9 does show, however, that the magnitude of the artifact varies appreciably depending on the particularities of the system, and that the spurious flow is most pronounced for hydrophobic

4.3. Simulation details for planar surfaces For inhomogeneous electric fields resulting from a set of point charges, there are a few additional simulation settings that may influence the resulting water dynamics [50]. To speed up the simulation, each molecule carries a list of molecules with which it has significant interaction, the socalled neighbor list. Reducing the update frequency of this list speeds up the calculation, while increasing the risk of inaccuracies. In figure 8(a) we show the ζ -potential of a planar hydrophobic diamond (double face-centered cubic lattice of Catoms, contact angle θ = 111◦ ) in a homogeneous electric field as a function of the number of steps between neighbor list updates (nstlist). The Lennard-Jones simple cutoff is set to rc = 0.8 nm. A value of nstlist = 1 means the list is updated every step. From figure 8(a), we conclude that for a spatially constant electric field the neighbor list update frequency does not have a significant influence on the ζ -potential of the solute. In the case of numerical artifacts, a dependence on the size of the integration time step may be expected as well. In figure 8(b), the spurious ζ -potential of a hydrophobic 8


D J Bonthuis et al

Figure 11. Radial velocity inside a (16, 16) carbon nanotube without reservoir. (a) Electro-osmotic velocity calculated with binning at t (dashed red line) and the correct binning at t − δt/2 (solid blue line). (b) Graph adapted from [24] showing a pressure-driven radial velocity in the same (16, 16) carbon nanotube.

Figure 10. Axial flux ( x1 direction) through a (10, 0) nanotube (a) and a (16, 0) nanotube (b) at E 1 = 1.0 V nm−1 using a simple cutoff at rc = 1 nm. Shown in red and marked with a square is the flux calculated using the straightforward counting of molecules, shown in blue and marked with a circle is the conditional flux, where molecules exiting the tube are counted only if they previously entered from the other side. The simulation boxes are identical to the ones described in figures 5 and 6. We use a time step of 1 fs and neighbor list update every 50 steps.

behavior. From the unconditional flux, the bursting is far less pronounced than in the filling curves from Hummer et al [11].

systems with a low average density of water molecules, a class of systems to which carbon nanotubes belong as well.

4.5. Radial flow in compressible fluids In a related publication, a non-zero radial velocity component is observed in pressure-driven flow inside a (16, 16) carbon nanotube [24], see figure 11(b), in stark contrast with continuum Poiseuille flow predictions. The non-zero radial flow is claimed to be related to the non-homogeneous density profile. Let us analyze the effect of compressibility and a space-dependent density on the radial velocity profile. The continuity equation in its most general form is given by

4.4. Definition of flux In the first paper on water in carbon nanotubes, Hummer et al found that the filling of nanotubes goes in bursts [11]. This bursting behavior seems to be manifested also in the simulation trajectories of [23], particularly in the (10, 0) nanotube. Let us discuss the way in which bursts in the cumulative flux emerge in the original analysis of [23]. Due to conservation of mass, the average flux through a given cross section of the carbon nanotube does not depend on the position of the cross section. Therefore, without restricting the generality of our arguments, we define the flux as the number of water molecules passing the edge of the simulation box. This we call the cumulative flux. In [23], on the other hand, the flux is defined as the number of particles exiting the tube at one end that have previously entered at the other end. This definition we denote as the conditional cumulative flux and it implicitly contains a tube length dependence. We perform simulations of 40 ns with a 5 ns equilibration period using a 1.0 nm simple cutoff scheme and an electric field of 1 V nm−1 . In figures 10(a) and (b) we show the results for two tube diameters, (10, 0) and (16, 0) respectively, both connecting two reservoirs. The results have been analyzed using both the straightforward counting of molecules passing the edge of the box (denoted as flux), and using the definition from [23] (denoted as conditional flux). For long times, both definitions of the flux converge to the same value, but for shorter times, there are distinct differences. The first, obvious effect of the conditional definition of the flux is that the first particles exiting the tube, which are the particles that were in the tube from the start, do not contribute to the conditional flux since they did not previously enter the tube. This explains the initial 4–8 ns waiting time and suggests that the conditional definition of the flux is not very efficient in terms of computer time. The second effect is a suppression of noise: particles moving back and forth are only counted if their motion persists for the full tube length. This explains the striking lack of noise in the results shown in [23], as well as a large part of the alleged bursting

∂ρ + ∇ · (ρu) = 0. ∂t

(29)

When we take ρ(x 3 ) and u(x 3 ) depending on the radial coordinate, the radial component of equation (29) reduces in steady state to

u 3 (x 3 )∇3 ρ(x 3 ) + ρ(x 3 )

1 ∇3 (x 3 u 3 (x 3 )) = 0. x3

(30)

Multiplying by x 3 and integrating once over x 3 yields ! ! x 3 u 3 (x 3 )∇3 ρ(x 3 ) dx 3 + ρ(x 3 )∇3 (x 3 u 3 (x 3 )) dx 3 = c,

(31) with c an integration constant. Integration by parts of the first term on the left-hand side gives

x 3 u 3 (x 3 )ρ(x 3 ) = c.

(32)

Equation (32) fixes u 3 (x 3 ) as a function of ρ(x 3 ). When the density ρ(x 3 ) vanishes at the wall of the tube, the velocity should remain finite, showing that the only physical solution is c = 0 and therefore u 3 (x 3 ) = 0, also for compressible fluids of variable density. To find out why the simulations of [24] exhibit a non-zero radial flow profile, we calculate the radial velocity u 3 (x 3 ) in simulations of a (16, 16) nanotube from radially binned water molecule positions x 3 (t) according to u 3 (x 3 [t − τ ]) = (x 3 [t] − x 3 [t − δt])/δt with δt the integration time step. The correct positioning of the velocity requires choosing τ = δt/2, and in this case u 3 vanishes, as shown in figure 11(a) (solid line), 9


D J Bonthuis et al

in agreement with standard hydrodynamic theory. Incorrectly using τ = 0 gives a spurious non-vanishing u 3 profile (broken line), similar to profiles reported in [24]. The reason for this is that a water molecule that at time t happens to be close to the surface is, due to the impenetrability of the surface, likely to originate from a position further away from the surface at time t − δt .

[11] Hummer G, Rasaiah J C and Noworyta J P 2001 Nature 414 188 [12] Gong X et al 2007 Nature Nanotechnol. 2 709 [13] Král P and Tománek D 1999 Phys. Rev. Lett. 82 5373 [14] Insepov Z, Wolf D and Hassanein A 2006 Nano Lett. 6 1893 [15] Longhurst M J and Quirke N 2007 Nano Lett. 7 3324 [16] Parsegian A 1969 Nature 221 844 [17] Bonthuis D J et al 2006 Phys. Rev. Lett. 97 128104 [18] Zhang J, Kamenev A and Shklovskii B I 2005 Phys. Rev. Lett. 95 148101 [19] Doyle D A et al 1998 Science 280 69 [20] Kienker P K, Degrado W F and Lear J D 1994 Proc. Natl Acad. Sci. USA 91 4859 [21] Kamenev A, Zhang J, Larkin A I and Shklovskii B I 2006 Physica A 359 129 [22] Li J 2007 Proc. Natl Acad. Sci. USA 104 3687 [23] Joseph S and Aluru N R 2008 Phys. Rev. Lett. 101 064502 [24] Joseph S and Aluru N R 2008 Nano Lett. 8 452 [25] Kim D and Darve E 2009 J. Colloid Interface Sci. 330 194 [26] Joseph S and Aluru N R 2006 Langmuir 22 9041 [27] Knecht V et al 2008 J. Colloid Interface Sci. 318 477 [28] Barrat J and Bocquet L 1999 Phys. Rev. Lett. 82 4671 [29] Joly L, Ybert C, Trizac E and Bocquet L 2004 Phys. Rev. Lett. 93 257805 [30] Vinogradova O I 1995 Langmuir 11 2213 [31] Mamatkulov S I, Khabibullaev P K and Netz R R 2004 Langmuir 20 4756 [32] Chandler D 2007 Nature 445 831 [33] De Gennes P G 2002 Langmuir 18 3413 [34] Huang D M, Sendner C, Horinek D, Netz R R and Bocquet L 2008 Phys. Rev. Lett. 101 226101 [35] Majumder M, Chopra N, Andrews R and Hinds B J 2005 Nature 438 44 [36] Sendner C et al 2009 Langmuir 25 10768 [37] de Groot S R and Mazur P 1962 Non-Equilibrium Thermodynamics (Amsterdam: North-Holland) [38] Lyklema J 2003 Colloids Surf. A 222 5 [39] Joly L, Ybert C, Trizac E and Bocquet L 2006 J. Chem. Phys. 125 204716 [40] van der Heyden F H J, Bonthuis D J, Stein D, Meyer C and Dekker C 2006 Nano Lett. 6 2232 [41] van der Heyden F H J, Bonthuis D J, Stein D, Meyer C and Dekker C 2007 Nano Lett. 7 1022 [42] Ren Y and Stein D 2008 Nanotechnology 19 195707 [43] Stein D, Kruithof M and Dekker C 2004 Phys. Rev. Lett. 93 035901 [44] Born M 1920 Z. Phys. 1 221 [45] Condiff D W and Dahler J S 1964 Phys. Fluids 7 842 [46] Bonthuis D J et al 2010 Langmuir 26 12614 [47] Bonthuis D J et al 2009 Phys. Rev. Lett. 103 144503 [48] van der Spoel D et al 2005 J. Comput. Chem. 26 1701 [49] Plimpton S 1995 J. Comput. Phys. 117 1 [50] Wong-ekkabut J, Miettinen M S, Dias C and Karttunen M 2010 Nature Nanotechnol. 5 555 [51] Hess B, Kutzner C, van der Spoel D and Lindahl E 2008 J. Chem. Theory Comput. 4 435 [52] Zuo G, Shen R, Ma S and Guo W 2010 ACS Nano 4 205 [53] Janeˇcek J and Netz R R 2007 Langmuir 23 8417

5. Conclusions We study flow through carbon nanotubes, focusing on the implications of the modified electrostatic and hydrodynamic boundary conditions inherent to the small scale of nanochannels. The extremely large slip lengths measured in carbon nanotubes greatly reduce the fluidic resistance to pressuredriven flow and at the same time enhance the efficiency of electrokinetic energy conversion. Obtaining the non-zero ion distribution necessary for electrokinetic pumping may not be easy, however, due to the charge neutrality of the nanotube as well as the self-energy barrier keeping ions from entering small capillaries. Fixed charges and external electric fields may enhance or reduce the translocation rate of ions, but we show they cannot be used to pump pure water. Simulation techniques have proven to be useful tools in researching the properties of nanofluidic systems. In particular cases, however, molecular dynamics simulations produce spurious kinetics, most notably when using a combination of electric fields and pure water at interfaces in the simulation package GROMACS.

Acknowledgment We thank the Deutsche Forschungsgemeinschaft for funding via SPP 1164.

References [1] Catterall W A 1995 Annu. Rev. Biochem. 64 193 [2] Terry L J, Shows E B and Wente S R 2007 Science 318 1412 [3] Yao S, Herzog D E, Zeng S, Mikkelsen J C and Santiago J G 2003 J. Colloid Interface Sci. 268 143 [4] Yang J, Lu F, Kostiuk L W and Kwok D Y 2003 J. Micromech. Microeng. 13 963 [5] Daiguji H, Yang P, Szeri A J and Majumder A 2004 Nano Lett. 4 2315 [6] Olthuis W, Schippers B, Eijkel J and van den Berg A 2005 Sensors Actuators B 111 385 [7] Service R F 2006 Science 313 1088 [8] Stone H A, Stroock A D and Ajdari A 2004 Annu. Rev. Fluid Mech. 36 381 [9] Squires T M and Quake S R 2005 Rev. Mod. Phys. 77 977 [10] Kuyucak S, Andersen O S and Chung S-H 2001 Rep. Prog. Phys. 64 1427

10