Molecular Simulation of Ion Transport in Silica Nanopores - The

Jan 5, 2009 - Velocity profiles of K+ ions in the z-direction inside silica nanopores of Lz0 = 2 nm (top) and Lz0 = 4 nm (bottom) calculated from NEMD...
0 downloads 0 Views 2MB Size
Download PDF
J. Phys. Chem. B 2009, 113, 1041–1047

1041

Molecular Simulation of Ion Transport in Silica Nanopores Katsuhiro Shirono,† Naoya Tatsumi,‡ and Hirofumi Daiguji*,‡ National Institute of AdVanced Industrial Science and Technology (AIST), Tsukuba 305-8564, Japan, and Institute of EnVironmental Studies, Graduate School of Frontier Sciences, The UniVersity of Tokyo, Kashiwa 277-8563, Japan ReceiVed: June 20, 2008; ReVised Manuscript ReceiVed: September 17, 2008

Ion distribution and transport of KCl aqueous solutions at the junction of hydrophobic and hydrophilic regions inside silica nanopores were studied by using two kinds of molecular simulation: grand canonical Monte Carlo (GCMC) simulations and nonequilibrium molecular dynamics (NEMD) simulations. The nanopores were 2 nm diameter silica pores in which surface functional groups, -SiOH, had been modified by hydrophobic surface functional groups, -SiCH3, within three different lengths along the pore direction (z-direction), Lz0 ) 0, 2, and 4 nm. If Lz0 is long enough, water could not enter the hydrophobic region, but for all Lz0 studied here, water entered the hydrophobic region. When an external electric field was applied along the z-direction, ions could not pass through the hydrophobic region when the external electric field was less than a threshold level, E0, whereas the ionic current increased relatively linearly with increasing electric field strength above E0. In 2 nm diameter fluidic pores, the electrical potential barrier appeared at the junction between the hydrophilic and hydrophobic regions due to the difference in dipole moment of the surface functional groups, although the overall surface charge of the pore was neutral. This junction formed an electrical potential threshold, and the ionic current could be modulated by adjusting the external electric field. 1. Introduction The molecular basis for ion transport through nanoscale channels, in particular, through nanochannels in biological membranes, has been actively researched. This has led to the discovery of membrane proteins that work as selective ion channels1-3 and, recently, to the application of ion-channel sensors that mimick this phenomenon in biological4-6 and artificial synthetic7-10 nanochannels. X-ray structural analyses of the molecular arrangement11,12 as well as molecular simulations13-16 have provided key insight into this ion selectivity. Ion transport through nanoscale channels has also been studied using inorganic nanopores for the sensing of single DNA molecules, biomolecules, and chemical substances. Pores in biological membranes are typically less than 1 nm in diameter (diameter D < 1 nm), and ions and water molecules usually cannot penetrate though the channel simultaneously. In contrast, inorganic pores used in ion transport studies to date are about 3 nm or larger in diameter. For example, Li et al.7 fabricated nanopores (D ≈ 5 nm) in a Si3N4 membrane by using an ion beam and observed a decrease in ionic current when doublestranded DNA molecules were translocated through the nanopore, suggesting that ionic current was blocked when the pore was occupied by double-stranded DNA molecules. Danelon et al.10 functionalized nanopores (D ≈ 4 nm) in a Si3N4 membrane with silicon oxide deposition and found asymmetry in ionic current-voltage curves with the sign of the applied potential. They reported that the asymmetry might be caused by only one side of the pore being coated with silicon oxide or by a slightly asymmetric radius profile along the channel axis. Fan et al.17 integrated silica nanotubes (D ≈ 50 nm) into microfluidic systems and observed the modulation of ionic current when * To whom correspondence should be addressed: e-mail [email protected], Ph +81 4 7136 4658, Fax +81 4 7136 4659. † National Institute of Advanced Industrial Science and Technology. ‡ The University of Tokyo.

DNA molecules were translocated through the nanotube. They reported that a transition from current decrease to current increase occurred with decreasing concentration of buffer solutions, suggesting that the current decrease was due to geometric blockage effect and the current increase to electrostatic charge effect. For molecular sensors based on geometric blockage, the channel size must be slightly smaller than the molecular size. However, asymmetry and enhancement of ionic current are due to electrostatic interactions between ions and channel walls, and thus for nanofluidic devices to utilize electrostatic interactions, the channel size is limited by the characteristic length scale of electrostatic interactions, that is, the Debye length. For monovalent electrolytes, the Debye length of the solution is 3 Å-100 nm depending on the ion concentration of 1-10-5 M. The characteristic length scales of transport phenomena in nanoscale channels depend on the changes in intermolecular forces near a solid surface. Electrostatic interaction energy (point charge-point charge interaction energy) is inversely proportional to the molecular distance r (i.e., r-1), dipole-point charge interaction energy is proportional to r-4-r-2, and dipole-dipole interaction energy is proportional to r-6-r-3. Dipole-point charge and dipole-dipole interactions are negligible when r is larger than 1 nm. But in 1 nm diameter channels, the dipole moment on the channel wall could affect the distribution and transport of ion and water. To obtain an adequate description of ion and water transport in channels whose diameter is smaller than 2 nm, the dipole-point charge and dipole-dipole interactions as well as the electrostatic interaction must be considered. However, because these interactions cannot be modeled by the continuum theory, the discreteness of molecules must be clarified in detail. In this study, molecular simulations were used to clarify the ion distribution and transport of KCl aqueous solutions at the junction of hydrophobic and hydrophilic regions inside silica

10.1021/jp805453r CCC: $40.75  2009 American Chemical Society Published on Web 01/05/2009

1042 J. Phys. Chem. B, Vol. 113, No. 4, 2009

Shirono et al.

nanopores. The nanopores were 2 nm diameter silica pores in which surface functional groups, -SiOH, had been modified by hydrophobic surface functional groups, -SiCH3, within three different lengths along the pore direction (z-direction), Lz0 ) 0, 2, and 4 nm. These nanopores are typical examples to understand dipole-point charge and dipole-dipole interactions between ion/water and pore walls because the hydrophilic surface is charge-polarized and capable of hydrogen bonding, whereas the hydrophobic surface is not. Ion and water transport through hydrophilic and hydrophobic pores whose diameter is smaller than 2 nm have been studied both experimentally and theoretically,18,19 and hydrophobic pores have usually been used to prevent water from entering. In this study, however, hydrophobic regions are so short that water could enter, but ion could not enter because an electric potential barrier appears at the junction of hydrophobic and hydrophilic regions due to the difference in dipole moment. Furthermore, the electric potential barrier could depend on Lz0; ion distribution and transport was studied for model silica nanopores with three different Lz0. This study aimed at (1) clarifying structural and dynamical properties of water and ions without an applied external electric field by using GCMC and equilibrium MD simulations and (2) clarifying transport properties of ions in the external electric field by using NEMD simulations. Furthermore, the ion distribution and transport properties calculated using these molecular simulations were compared with those calculated using continuum dynamics simulations, and the effect of the discreteness of molecules on the transport properties was clarified. 2. Methodology 2.1. Modeling of Silica Nanopores. A silica pore model based on a crystal structure of R-quartz was used in the molecular simulations in this study. The modeling entailed the following five-step procedure.20 In step 1, first, 8, 8, and 4 trigonal unit cells of R-quartz (Si3O6) were arrayed in the [100], [010], and [001] directions, respectively (see Figure 1a). Then, this R-quartz crystal was transformed from a rhombohedron to a rectangular parallelepiped (3.931 × 3.404 × 2.162 nm3) by relocating the positions of the unit cells. The x- and z-axes were selected as the direction parallel to the [100] and [001] directions, respectively (see Figure 1b). In step 2, Si atoms within a radius of 0.98 nm from the z-axis were removed. In step 3, O atoms that were covalently bonded with the two adjacent removed Si atoms were also removed. In step 4, H atoms were added to O atoms bonded with only one Si atom to prevent dangling bonds on the hydrophilic surface, whereas united atoms of CH3 were replaced with the O atoms on the hydrophobic surface. The length of the hydrophilic surface in the z-direction was 4 nm (exactly 4.324 nm), and that of hydrophobic surface was set at 0, 2, and 4 nm (exactly 0.0, 2.162, and 4.324 nm, respectively). In step 5, to minimize stress in the model pores, the canonical ensemble MD (NVT-MD) simulations were conducted with all atoms mobile at 300 K. 2.2. Potential Functions and Parameters. The model silica nanopores were composed of SiO2 and the surface groups of SiOH and SiCH3. The Tsuneyuki potential model21 was used for SiO2 (consisting of silicon (SiS) and oxygen (OS) atoms). The surface groups of SiOH (consisting of silicon (SiS*), oxygen (OS*), and hydrogen (HS*) atoms) and SiCH3 (consisting of silicon atom (SiS′) and a united atom of CH3 (CH3)) were assumed to be rigid bodies. The SPC/E water model22 was used for water (consisting of oxygen (OW) and hydrogen (HW) atoms), and Dang’s ion pair model23 was used for potassium (K+) and chloride (Cl-) ions. Table 1 shows the potential

Figure 1. (a) 8 × 8 R-quartz unit cells on the (001) plane and (b) R-quartz model pore. The length of the unit vector a is 4.9133 Å. The x- and z-axes were selected as the direction parallel to the [100] and [001] directions, respectively.

functions and parameters used in this study. The interaction between the atoms/united atoms involved two contributions: Coulomb interaction and either Born-Mayer-Huggins (BMH) interaction or Lennerd-Jones (LJ) interaction. The BMH potential function was used for the interaction involved with SiS, OS, SiS*, OS*, HS*, and SiS′ except for the interactions between CH3-SiS, CH3-OS, CH3-SiS*, CH3-OS*, and CH3-HS*. These latter five interactions (CH3-SiS, CH3-OS, CH3-SiS*, CH3-OS*, and CH3-HS*) along with the interactions between OW, HW, K+, Cl-, and CH3 were calculated using the LJ potential function. The parameters of the BMH potential functions were cited from refs 20, 21, and 24, except for a small repulsion term involved with HS*. The parameters of the LJ potential function were cited from refs 22, 23, and 25. The parameter cSiS′ was determined as the pair potential between SiS′ and OS is minimized at the distance of rSiS′-OS ) 1.415 Å because the pair potential between SiS and OS is minimized at rSiS-OS ) 1.415 Å in the Tsuneyuki potential model. The distances between SiS*-OS*, SiS*-HS*, OS*-HS*, and SiS′-CH3 were fixed at 1.61, 2.16, 1.0, and 1.0 Å, respectively. 2.3. Simulation Details. Initially, four pairs of K+ and Clions were set in the hydrophilic region. There was no water molecule inside nanopores on the initial condition, and water molecules were inserted into the pores in the GCMC simulations. The number of ions was set as the ion concentration below the solubility limit of KCl aqueous solutions. K+ ions were distributed randomly in the hydrophilic region, whereas Cl- ions were located near SiOH groups. Note that Cl- ions were located between two adjacent SiOH groups stably and did not move

Ion Transport in Silica Nanopores

J. Phys. Chem. B, Vol. 113, No. 4, 2009 1043

TABLE 1: Potential Functions and Parameters Used in This Studya ID

(united) atom

1 2 3 4 5 6 7 8 9

OS, OS* SiS, SiS* SiS′ HS* CH3 K+ ClOW HW

qi

ai (Å)

bi (Å)

ci (eV1/2 Å3)

ref

εi (eV × 102)

σi (Å)

ref

-1.2 +2.4 +1.8 +0.6 0 +1.0 -1.0 -0.8476 +0.4238

2.0474 0.8688 0.8688 0.0

0.175 66 0.032 85 0.032 85 0.0

14.657 4.828 5.497 0.0

21 21 in text

0.5272

3.402

25

1.908 1.95 1.694

0.131 0.18 0.1179

0.7588 0.4336 0.4336 0.6739

3.87 3.331 4.401 3.166

25 23 23 22 22

3.1243 6.7835 5.21

24 24 20

a OS, OS*, and OW are O atoms of silica pores, SiOH groups, and water molecules, respectively. SiS, SiS*, and SiS′ are Si atoms of silica pores, SiOH groups, and SiCH3 groups, respectively. HS* and HW are H atoms of SiOH groups and water molecules, respectively. For the (united) atom pairs designated by the ID numbers 1-1, 1-2, 1-3, 1-4, 1-5, 1-6, 1-7, 1-8, 2-2, 2-3, 2-4, 2-6, 2-7, 3-3, 3-6, 3-7, 4-6, and 4-7 in the table, Coulomb, uij(r) ) qiqje2/r, and BMH potential functions, uij(rij) ) f0(bi + bj) exp[(ai + aj - rij)/(bi + bj)] - cicj/rij6, were used, whereas for the (united) atom pairs of 5-5, 5-6, 5-7, 5-8, 6-6, 6-7, 6-8, 7-7, 7-8, and 8-8, Coulomb, and LJ potential functions, uij(rij) ) 4(εiεj)1/2{[(σi + σj)/2rij]12 - [(σi + σj)/2rij]6}, were used. For the other pairs, only Coulomb potential functions were used. Here, f0 ) 1 (kcal Å-1 mol-1) ) 4.3384 × 10-2 (eV Å-1). e ) 1.602176 × 10-19.

appreciably from the initial positions in the following simulations. Possible reasons for the immobility of Cl- ions are as follows: (1) the deprotonation reaction (SiOH f SiO- + H+) was neglected; (2) high surface number density of SiOH groups was assumed; (3) the potential energy functions of OS* and SiS* were assumed to be the same as those of OS and SiS, and the partial charge of HS* was determined to become electrically neutral in the entire silica pore models. To clarify the ion distribution and transport in more realistic silica pores, more precise modeling of the surface groups including the deprotonation reaction should be crucial. The GCMC simulations of a KCl aqueous solution were conducted for 300 K and 1 kPa. In the GCMC simulations, the trials of the insertion and the deletion were conducted only for water molecules, and the trials of the motion were conducted for water molecules, ions, and SiOH and SiCH3 groups. These GCMC simulations yielded the equilibrium structures of SiOH and SiCH3 groups as well as those of water and ions. The calculated number of water molecules at Lz0 ) 0, 2, and 4 nm was 309, 430, and 524, respectively. The equilibrium structures obtained by the GCMC simulations were then used as the initial structure for the MD simulations. Both the initial transitional and rotational momentum were set to zero for all atoms/united atoms and molecules. In the NEMD simulations, the electric field was applied in the z-direction, and the magnitude, Eex, was set at six different levels: 0.0 (equilibrium MD), 0.1, 0.5, 0.75, 1.0, and 1.5 V/nm. Temperature was regulated at 300 K by using a Berendsen thermostat with a time constant of 200 fs. To avoid biasing the velocity in the z-direction, only the velocity components normal to the flow direction, Vx and Vy, were regulated using the thermostat. Periodic boundary conditions were applied to the pore model described in section 2.1, and the electrostatic interactions were calculated using the Ewald method. Both the cutoff distance of van der Waals interaction and the real-space term of the Coulomb interaction were 12.5 Å. To minimize calculation time, all atoms except for water molecules, ions, and SiOH and SiCH3 groups were fixed at the initial position determined by the procedure described in section 2.1.Water molecules, ions, and SiOH and SiCH3 groups were allowed to equilibrate at 300 K for 200 ps. In the MD simulations, the time step was 2 fs and the total run time was 1200 ps. To clarify the dielectric properties of water confined in nanospace, the radial and axial components of dielectric constants of water, εr and εz, were calculated by the following equations: [(εr - 1)(εs + 1)]/(εr + εs) ) 2π∑R)x,y(〈MR2〉 -

〈MR〉2)/kTV and εz ) 1 + [4π(〈Mz2〉 - 〈Mz〉2)]/kTV, where εs is the external dielectric constant, MR and Mz are the total dipole moment of water molecules in the R () x, y) and z directions, and V is the pore volume.26 The external dielectric constant was assumed to be εs ) ∞, and V was calculated by assuming that the pore diameter was 1.81 nm. To understand the transport properties of K+ ions, the total potential energy profile of K+ ions from all atoms/ions (Utotal) and the decomposed potential energy profile from each atom/ ion, namely, pore atom (Upore), Cl- ion (UCl-), water (Uwater), and K+ ion (UK+), were calculated in the z-direction. Each potential energy was calculated as a summation of pair potential between K+ ion and each atom/ion. In this calculation, electrostatic energy was calculated with Ewald sum method, and the self-energy term was classified to UK+. Furthermore, Utotal was also decomposed into the van der Waals interaction C energy, UW total, and the Coulomb interaction energy, Utotal. 3. Results and Discussion 3.1. Equilibrium States of Water and Ions in Silica Nanopores. Figure 2a-d shows cross-sectional snapshots of KCl aqueous solutions in silica nanopores produced by the GCMC simulations. If Lz0 is long enough, water could not enter the hydrophobic region, but for all Lz0 studied here, water entered the hydrophobic region. Furthermore, in the hydrophobic region, water molecules were located near the center of the pore apart from the hydrophobic wall. Figure 3 shows the water density profiles in the z-direction inside the silica pores for Lz0 ) 0, 2, and 4 nm. The water density was averaged over the cross section of the pore, and the pore diameter was assumed to be 1.81 nm. In the hydrophobic regions, the average density of water in each pore was about 0.60 g/cm3, whereas in the hydrophilic regions, it was about 0.85 g/cm3. Although the density of water in the hydrophobic region was less than that in the hydrophilic region, the hydrophobic region was not long enough to interrupt the column of a KCl aqueous solution. The calculated εr and εz inside the silica nanopores of Lz0 ) 0, 2, and 4 nm were (εr, εz) ) (11.32, 46.28), (10.92, 20.31), and (9.53, 22.54), respectively. For all Lz0, εr was smaller than εz, and εz was smaller than that in the bulk SPC/E water, ε ) 67 ( 10.27 Figure 4a,b shows the radial distribution functions of a K+ ion around a Cl- ion (i.e., gK+-Cl-) and around OS* (i.e., gK+-OS*) inside the silica nanopore of Lz0 ) 0 nm. The first peaks are respectively located at rK+-Cl- ≈ 0.33 nm and rK+-OS* ≈ 0.65 nm. Because the ionic radiuses of K+ and Cl- are about 1.38 and 1.81 Å, respectively, the first peak position of gK+-Cl- is

1044 J. Phys. Chem. B, Vol. 113, No. 4, 2009

Shirono et al.

Figure 2. Snapshots of KCl aqueous solutions confined in 1.96 nm diameter silica pores modified by hydrophobic groups within three different fixed lengths along the pore axis (z-axis), Lz0. Cross section of a model pore with (a) Lz0 ) 0 nm vertical to z-axis and (b) Lz0 ) 0 nm, (c) Lz0 ) 2 nm, and (d) Lz0 ) 4 nm along the z-direction. In (b)-(d), hydrophilic regions are at both ends, and the hydrophobic region is in the central part along the z-axis. Each cross section is 0.6 nm thick, and the number of water molecules in each pore is Nwater ) 309, 430, and 524, respectively.

Figure 3. Water density profiles in the z-direction inside silica pores of three different Lz0. The water density was averaged over the cross section of the pore, and the pore diameter was assumed 1.81 nm. Dashed lines represent the boundary between hydrophilic and hydrophobic regions.

almost the same as the sum of the ionic radiuses of K+ and Cl-, i.e., 3.19 Å. This result suggests that K+ and Cl- ions could form a directly linked structure. In contrast to gK+-Cl-, gK+-OS* showed a broad distribution. Because water was hydrated around OS* (the first peak position of gOS*-OW is about 0.275 nm20), K+ could approach the hydration water around OS* but could not form a directly linked structure to OS*. 3.2. Ion Transport in Silica Nanopores. In all the NEMD simulations studied here, K+ ions were transported through the silica nanopores, whereas Cl- ions hardly moved from their initial positions. In the following sections, only the transport properties of K+ ions are therefore discussed. Figure 5a shows the number of K+ ions passing through the junction between

Figure 4. Radial distribution functions of a K+ ion around (a) Clion, gK+-Cl-, and (b) OS*, gK+-OS*, inside a silica nanopore of Lz0 ) 0 nm at Eex ) 0 V/nm. Radial distribution function was normalized using an average density of K+ ions inside the nanopore. Dashed lines represent the first peak positions of rK+-Cl- and rK+-OS*.

the hydrophilic and hydrophobic regions during 1 ns calculation as a function of the applied external electric field (NK+* - Eex) for three different lengths of the hydrophobic region, Lz0 ) 0, 2, and 4 nm. For Lz0 ) 2 and 4 nm, NK+* ) 0 at Eex e 0.75 V/nm and NK+* increases with increasing Eex at Eex g 0.75 V/nm, suggesting that a threshold exists for the external electric field, E0, for a K+ ion passing through the junction and that the magnitude of E0 does not depend on Lz0. Figure 5b shows the average velocity of K+ ions as a function of Eex (VK+ - Eex). Also shown in this figure are the curves for VK+ ) eEex/ξ (solid line) and VK+ ) e(Eex - E0)/ξ (dashed line), where ξ is the friction coefficient of a K+ ion and eEex/ξ represents the velocity of K+ ions in the bulk water in an applied external electric field

Ion Transport in Silica Nanopores

Figure 5. (a) Number of K+ ions passing through the z ) 0 plane, NK+*, vs external electric field, Eex, and (b) average velocity of K+ ions, VK+, vs external electric field, Eex, inside silica nanopores of different Lz0. Also shown are curves for VK+ ) eEex/ξ (solid line) and VK+ ) e(Eex - E0)/ξ (dashed line). VK+ was calculated from VK+ ) NK+*/ FA, where F is the number density of K+ ions, A is the sectional area of the pore, and t is calculation time. Dashed line represents the threshold of external electric field, E0 ) 0.75 V/nm.

Eex. For all the pores, the gradients of VK+ with respect to Eex, ∂VK+/∂Eex, are almost the same as the gradient of e/ξ at Eex g E0. Furthermore, the velocities of water molecules were calculated to examine the effect of electroosmotic flow. The velocity of water molecules is about 0.1-7.5 m/s, and it is only about 5% of the velocity of K+ ion. 3.3. Energetic Properties. The transport properties of K+ ions were studied here by evaluating the potential energy profiles of K+ ions in the z-direction. Figure 6 shows the profiles of the total potential energy of K+ ions, Utotal, in the z-direction. Because the thermal oscillation energy (kT ≈ 0.026 eV) is much smaller than the potential energy, E0 can be determined by the maximum gradient of the total potential energy with respect to z, (∂Utotal/∂z)max. The calculated (∂Utotal/∂z)max is about 0.75 eV/ nm at the junction between the hydrophilic and hydrophobic regions, and this value is the same for Lz0 ) 2 and 4 nm when Eex ) 0.75 eV/nm. The value of (∂Utotal/∂z)max ) 0.75 eV/nm agrees well with the value of E0, namely, 0.75 V/nm (see Figure 5). E0 is independent of Lz0, suggesting that (∂Utotal/∂z)max is determined by the local structure at the junction between the hydrophilic and hydrophobic regions. Furthermore, ∂VK+/∂Eex is almost the same in all pores at Eex > E0 (Figure 5b), suggesting that VK+ is determined by the electric field at the junction between the hydrophilic and hydrophobic regions. When Eex is applied, K+ ions accumulate at the junction and the potential barrier weakens, although the K+ ions cannot overcome the potential barrier at Eex < E0. With increasing Eex above E0, the K+ ions can overcome the potential barrier, and the effective electric field for K+ ions becomes E ) Eex - E0. Therefore, ∂VK+/∂Eex ≈ e/ξ, that is, VK+ ≈ e(Eex - E0)/ξ (Figure 5b). The slight nonlinearity of VK+-Eex curves is possibly due to the change in ξ with respect to Eex. Although the overall charge of the pore surface was neutral, this junction could form an electrical potential threshold, and thus the ionic current could be modulated by adjusting Eex.

J. Phys. Chem. B, Vol. 113, No. 4, 2009 1045

Figure 6. Total potential energy profiles of K+ ions, Utotal, in the z-direction inside silica nanopores of different Lz with an applied external electric field of Eex ) 0.75, 1.0, and 1.5 V/nm. Dashed lines represent the boundary between hydrophilic and hydrophobic regions, and the arrow shows the direction of Eex.

If Utotal is decomposed into the van der Waals interaction C energy, UW total, and the Coulomb interaction energy, Utotal, then C W Utotal is almost the same as Utotal because Utotal is only about 0.01 eV at any position along the z-axis. Therefore, the potential C . energy barrier at the junction mainly depends on Utotal Furthermore, the Coulomb interaction energy can be categorized into three types of interaction energy: (i) ion-dipole moment interaction energy, i.e., Upore, (ii) ion-ion interaction energy, i.e., Uion ) UK+ + UCl-, and (iii) ion-solvent interaction energy, i.e., Uwater. Figure 7 shows Utotal, Upore, UCl-, Uwater, and UK+ in the z-direction inside the silica nanopores of Lz0 ) 2 and 4 nm. Upore is relatively independent of Uwater, but Uion has negative correlation with Uwater, implying that the dielectric effect of water does not act on the ion-dipole moment interaction but does act on the ion-ion interaction. Possible reasons for this weak dielectric effect of water on the ion-dipole moment interaction are as follows: (a) motion of water molecules near the pore wall is constrained, (b) only one or two water molecules can be positioned between K+ ions and pore atoms, and (c) ion-dipole moment interaction is intrinsically weaker than ion-ion interaction. The weak dielectric effect of water causes a strong interaction between K+ ions and pore atoms. As a result, the local structure of K+ ions around pore atoms as well as the distribution of K+ ions in the pore direction affects the formation of potential barrier at the junction between the hydrophilic and hydrophobic regions. 3.4. Kinetic Features of Ion Transport. Figure 8 shows the velocity profiles of K+ ions in the z-direction inside silica nanopores of Lz0 ) 2 and 4 nm obtained from the NEMD simulations(symbols)andone-dimensionalPoisson-Nernst-Plank (PNP) equations (solid lines). In both results, the velocity of K+ ions increases in the hydrophobic region and decreases in the hydrophilic region and is minimum at the junction. In the calculation of one-dimensional PNP equations, the external electric potential, U0, was expressed as U0 ) -Eexz + Upore + UCl- + Uwater calculated with NEMD in the z-direction, where Eexz is the external electric field in the z-direction and Upore,

1046 J. Phys. Chem. B, Vol. 113, No. 4, 2009

Shirono et al. shielding effect of water to K+ ions has a negligible effect on the entire potential energy profile of Uwater. The detailed calculation method to solve PNP equations is shown in the Appendix. When the average density of K+ ions, njK+, was calculated by assuming the diameter of the pore to be 1.0 nm and DK+ was assumed to be the bulk value of 1.96 × 10-9 m2/ s28 and εz was assumed to be 20 (near the calculated values of 20.31 and 22.54 for Lz0 ) 2 and 4 nm, respectively), the velocity profiles of K+ ions obtained from one-dimensional PNP equations agree well with those obtained from NEMD simulations. When the external potential energy profile of K+ ions, U0, was determined by the atomistic interactions between K-pore atoms, K+-Cl-, and K+-water, and the parameters of DK+, εz, and njK+ were assumed to be the values not in the bulk but in silica nanopores, the transport property of ions can be obtained from one-dimensional PNP equations. 4. Conclusions

Figure 7. Potential energy profiles of K+ ions passing through silica nanopores of Lz0 ) 2 (top) and 4 nm (bottom) with an applied external electric field Eex ) 1.5 V/nm. Also shown are total potential energy profile (Utotal) and decomposed potential energy profiles from each atom/ ion: pore atom (Upore), Cl- ion (UCl-), water (Uwater), and K+ ion (UK+). Dashed lines represent the boundary between hydrophilic and hydrophobic regions, and the arrow shows the direction of Eex.

Figure 8. Velocity profiles of K+ ions in the z-direction inside silica nanopores of Lz0 ) 2 nm (top) and Lz0 ) 4 nm (bottom) calculated from NEMD simulations (symbols) and one-dimensional PNP equations (lines) with an applied external electric field Eex ) 1.0 V/nm (blue) and Eex ) 1.5 V/nm (red). Dashed lines represent the boundary between hydrophilic and hydrophobic regions, and the arrow shows the direction of Eex.

UCl-, and Uwater are the electric potential energies between K+-pore atom, K+-Cl-, and K+-water, respectively (see Figure 7). The electric potential of Uwater includes the shielding effect of water to K+ ions. In the calculation of one-dimensional PNP equations, the shielding effect was considered in ε as well as U0, and thus the effect was double-counted. However, in this study, Uwater was used without any correction because the

Ion distribution and transport of KCl aqueous solutions at the junction of hydrophobic and hydrophilic regions inside silica nanopores were studied by using molecular simulations. The nanopores were 2 nm diameter silica pores in which surface functional groups, -SiOH, had been modified by hydrophobic surface functional groups, -SiCH3, within three different lengths along the pore direction (z-direction), Lz0 ) 0, 2, and 4 nm. GCMC simulations and NEMD simulations of the ion distribution and transport were conducted for six different external electric fields, Eex ) 0.0, 0.1, 0.5, 0.75, 1.0, and 1.5 V/nm, and for a temperature of 300 K. The results yielded the following conclusions: 1. In the GCMC simulations, water entered the hydrophobic region as well as the hydrophilic region for all Lz0 studied here. In the hydrophobic region, water molecules were located near the center of the pore apart from the hydrophobic wall. Clions are located near SiOH groups stably and do not move appreciably from the initial positions in the simulations. In contrast, K+ ions diffuse inside the hydrophilic region and strongly interact with the surface molecules but do not enter the hydrophobic region unless an external electric field Eex is applied. 2. When Eex is applied in the pore direction, ions cannot pass through the hydrophobic region when Eex is less than a threshold E0 ) 0.75 V/nm, whereas when Eex > E0 ) 0.75 V/nm, the ionic current increases relatively linear with increasing Eex. The values of E0 and ∂VK+/∂Eex are independent of Lz0 because the local structure of the junction between hydrophilic and hydrophobic region dominates the electrical potential barrier for K+ ions. 3. The dielectric property of water affects the ion-ion interaction but does not affect the ion-dipole moment interaction, resulting in a strong interaction between K+ ions and pore atoms. The local structure of K+ ion around pore atoms as well as the distribution of K+ ions in the pore direction affects the formation of potential barrier at the junction between the hydrophilic and hydrophobic regions. 4. One-dimensional PNP equations were solved by assuming the external potential energy of K+ ions, U0, that was obtained from NEMD simulation and the parameters of njK+, DK+, and εz that were estimated from MD simulation. Velocity profiles of K+ ions calculated using the continuum theory agree well with those calculated using NEMD simulations. Appendix. Calculation Method for Poisson-Nernst-Plank (PNP) Equations The Poisson-Nernst-Plank equations are expressed as follows:

Ion Transport in Silica Nanopores

d2(U(z) - U0) dz2

)-

J. Phys. Chem. B, Vol. 113, No. 4, 2009 1047

enK+(z) ε0εz

(1)

dJK+ )0 dz

(2)

where U(z) is the electric potential, U0 is external electric potential, nK+(z) is number density of K+ ions, e is electronic charge, ε0 is permittivity, ε is dielectric constant, and JK+ is mass flux of K+ ions. JK+ is given by JK+ ) -DK+{∂nK+(z)/∂z + [enK+(z)/kT][∂U(z)/∂z]}, where DK+ is diffusivity of K+ ions, k is Boltzmann constant, and T is temperature. The boundary conditions are as follows:

dU dU ) , dz z)0 dz z)Lz dnK+ dnK+ , z)0 ) dz dz z)Lz

U(0) ) U(Lz) + EexLz

njK+ )

∫0L nK (z) dz/Lz z

+

(3)

(4)

where njK+ is the average number density of K+ ions. Equations 1 and 2 were solved under boundary conditions (3) and (4) using a finite difference algorithm, yielding U(z) and nK+(z). The velocity profiles of K+ ions were calculated as VK+(z) ) JK+/ nK+(z). References and Notes (1) Hodgkin, A. L.; Huxley, A. F.; Katz, B. J. Physiol. 1952, 116, 424–448. (2) Neher, E.; Sakmann, B. Nature (London) 1976, 260, 799–802. (3) MacKinnon, R. Nature (London) 1991, 350, 232–235. (4) Movileanu, L.; Howorka, S.; Braha, O.; Bayley, H. Nat. Biotechnol. 2000, 18, 1091–1095.

(5) Mathe, J.; Aksimentiev, A.; Nelson, D. R.; Schulten, K.; Meller, A. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 12377–12382. (6) Alcaraz, A.; Ramirez, P.; Gimenez, E. G.; Lopez, M. L.; Andrio, A.; Aguilella, V. M. J. Phys. Chem. B 2006, 110, 21205–21209. (7) Li, J.; Stein, D.; McMullan, C.; Branton, D.; Aziz, M. J.; Golovchenko, J. A. Nature (London) 2001, 412, 166–169. (8) Li, J.; Gershow, M.; Stein, D.; Brandin, E.; Golovchenko, J. A. Nat. Mater. 2003, 2, 611–615. (9) Chen, P.; Gu, J.; Brandin, E.; Kim, Y.-R.; Wang, Q.; Branton, D. Nano Lett. 2004, 4, 2293–2298. (10) Danelon, C.; Santschi, C.; Brugger, J.; Vogel, H. Langmuir 2006, 22, 10711–10715. (11) Brejc, K.; Dijk, W. J.; Klaassen, R. V.; Schuurmans, M.; Oost, J.; Smit, A. B.; Sixma, T. K. Nature (London) 2001, 411, 269–276. (12) Zhou, Y.; Morais-Cabral, J. H.; Kaufman, A.; MacKinnon, R. Nature (London) 2001, 414, 43–48. (13) Berneche, S.; Roux, B. Nature (London) 2001, 414, 73–77. (14) Compoint, M.; Carloni, P.; Ramseyer, C.; Girardet, C. Biochim. Biophys. Acta 2004, 1661, 26–39. (15) Liu, Z.; Xu, Y.; Tang, P. J. Phys. Chem. B 2006, 110, 12789– 12795. (16) Roux, B.; Karplus, M. Biophys. J. 1991, 59, 961–981. (17) Fan, R.; Karnik, R.; Yue, M.; Li, D.; Majumdar, A.; Yang, P. Nano Lett. 2005, 5, 1633–1637. (18) Steinle, E D.; Mitchell, D T.; Wirtz, M.; Lee, S B.; Young, V Y.; Martin, C R. Anal. Chem. 2002, 74, 2416–2422. (19) Beckstein, O.; Biggin, P. C.; Sansom, M S. P. J. Phys. Chem. B 2001, 105, 12902–12905. (20) Shirono, K.; Daiguji, H. J. Phys. Chem. C 2007, 111, 7938–7946. (21) Tsuneyuki, S.; Tsukada, M.; Aoki, H.; Matsui, Y. Phys. ReV. Lett. 1988, 61, 869–872. (22) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269–6271. (23) Dang, L. X. J. Am. Chem. Soc. 1995, 117, 6954–6960. (24) Kinoshita, T.; Mashimo, T.; Kawamura, K. J. Phys.: Condens. Matter 2005, 17, 1027–1035. (25) Hussain, H.; Titiloye, J. O. Microporous Mesoporous Mater. 2005, 85, 143–156. (26) Ballenegger, V.; Hansen, J.-P. J. Chem. Phys. 2005, 122, 114711. (27) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757–3766. (28) Lobo, V. M. M.; Ribeiro, A. C. F.; Verissimo, L. M. P. J. Mol. Liq. 1998, 78, 139–149.

JP805453R