Atomistic Description of Pressure-Driven Flow of ... - ACS Publications

Article pubs.acs.org/JPCC

Atomistic Description of Pressure-Driven Flow of Aqueous Salt Solutions through Charged Silica Nanopores Neil R. Haria† and Christian D. Lorenz*,‡ †

School of Engineering, University of Warwick, Coventry CV4 7AL, United Kingdom Theory & Simulation of Condensed Matter Group, Department of Physics, King’s College London, London WC2R 2LS, United Kingdom

‡

S Supporting Information *

ABSTRACT: Lab-on-a-chip devices and nanoscale porous membranes made from silica substrates are currently of significant interest for application to water desalination, biosensing, and chemical separation technologies. In each of these applications, the interaction between ionic solutions and the silica interface is key to developing design rules. In this paper, we present the results of extensive all-atom molecular dynamics simulations of the streaming current flow of 0.5 M NaCl and 0.5 M CaCl2 ionic solutions through cylindrical charged silica nanopores with diameters of 1.5, 2.0, 2.5, and 3.0 nm. We present the results from these simulations which provide a detailed description of the ion transport through the pores. We investigate the effect of pore size on ion pairing, ion hydration, and the structure of ions within the aqueous solutions. We also present the flux of the various species of the solutions through the various pores and the average current that is carried by the ions through the pores. In the 0.5 M NaCl systems, we observe that as the pore diameters increase the average current becomes increasingly positive, while for the 0.5 M CaCl2 system, the average current becomes increasingly negative as the pore diameter increases. This difference is a result from the fact that charge inversion is present in the CaCl2 systems but not in the NaCl systems.

■

metal oxides,47 and zeolites.48,49 Additionally, the dynamics of ionic solutions that are driven by electric fields,29,31,32,34−36,44,50−52 concentration gradients,22,45 temperature gradients,27 and pressure drops21,24,25,34,38,48,49 have been investigated using molecular dynamics simulations, and there has been a recent study in which the effect of competing osmotic and electric field gradients is studied.30 In the majority of the aforementioned simulation studies, the ionic solution contains simple monovalent or divalent salt solutions; however, there have been two initial studies of the transport of biological molecules in the form of a protein38 and single-stranded DNA52 using simulations. In this study, we are specifically interested in the transport and exclusion of ions when ionic solutions are driven through charged silica nanopores. Previous simulation studies have applied electric fields across charged silica slit pores in order to measure the transport properties of NaCl,34,53 KCl,54 and CaCl234 aqueous solutions in electroosmotic flow (EOF). Additionally, the transport properties of pressure-driven flow of NaCl and CaCl2 solutions through charged silica slit pores has also been measured.34 Transport of ionic solutions through

INTRODUCTION The emergence of lab-on-a-chip type devices and nanoscale porous membrane technology has resulted in an increase in the importance of understanding the interactions between fluids and substrates at the nanoscale.1−11 These nanofluidic systems are being utilized in a diverse range of application areas including water desalination,12−15 biosensing,16,17 and chemical separations,18 where in each case having a detailed understanding of the effect of the interactions between ionic solutions and synthetic nanopores plays a key role in the design of the systems. Subsequently, research into the transport of ions and biomolecules through nanochannels of various materials has increased, with a particular interest in obtaining an atomistic and molecular scale description of how the chemical and physical properties of the materials that make up the nanopore affect the structure and dynamics of ionic solutions within the nanopores. Molecular simulations are being employed extensively in order to provide atomistic detail of the structure and dynamics of these nanofluidic systems,8−10,19,20 and therefore simulations are able to complement experimental studies, which are not sensitive enough to measure these properties to the same detail. Simulations have been applied to study ionic solutions in contact with nanopores made from a wide range of materials including carbon nanotubes,12,21−27 graphene,28−32 silica,12−46 © 2015 American Chemical Society

Received: December 29, 2014 Revised: April 11, 2015 Published: April 30, 2015 12298

DOI: 10.1021/jp5129639 J. Phys. Chem. C 2015, 119, 12298−12311

Article

The Journal of Physical Chemistry C cylindrical charged silica nanopores when applying electric fields35,36,44 or concentration gradients45 has been studied with molecular dynamics simulations. In this paper, we report on the first atomistic molecular dynamics (MD) simulation study of the effect of the pore diameter of cylindrical charged silica nanopores on the structural and transport properties of ionic solutions driven by a pressure drop (i.e., streaming current) across the pore. The systems that have been studied consist of a nanopore that connects two reservoirs of ionic solution (see Figures 1 and 2),

Figure 1. Snapshots of the simulated systems such that the pore openings for (a) d = 1.5 nm, (b) d = 2.0 nm, (c) d = 2.5 nm, and (d) d = 3.0 nm can be seen. The water molecules are represented with small blue spheres, and the silica substrate is represented by small red, yellow, and white spheres for the oxygen, silicon, and hydrogen atoms, respectively, of which it consists. The large red spheres represent the bare oxygens which are the charged species on the pore interfaces.

Figure 2. Snapshots of the simulated systems with pores with diameters of (a) 1.5, (b) 2.0, (c) 2.5, and (d) 3.0 nm. In our systems, the pores are aligned such that they are parallel to the z-axis, and there is a periodic reservoir of aqueous salt solutions. The water molecules are represented with small blue spheres, and the silica substrate is represented by small red, yellow, and white spheres for the oxygen, silicon, and hydrogen atoms, respectively, of which it consists. The large red spheres represent the bare oxygens which are the charged species on the pore interfaces, and from these figures one can see the random distribution of them along the pore.

which has allowed us to also investigate ion exclusion of both monovalent (NaCl) and divalent (CaCl2) ionic solutions. In the Model and Methodology section, we provide a detailed description of the setup of our simulated systems and the simulation parameters we use during the pressure-driven flow simulations. Then, in the Results section, we discuss the structures of the ionic solutions within and outside of the pore, and we specifically focus on how the hydration of the ions and ion pairing change as they travel through the pore. Also, we report the transport properties of the streaming current flow of the ionic solutions through the nanopores. Finally, in the Discussion section, we present a discussion of our results and compare them to other related simulations studies of ionic solutions in nanopores.

pressure gradient. In order to construct the initial configurations used in this study, the approach used here is similar to that used in our previous study.44 The amorphous silica nanopores were constructed by starting with a bulk α-quartz crystal and then carrying out a series of MD simulations in order to convert the crystalline substrate into an amorphous silica substrate with two free interfaces in the z-dimension. First, the crystalline silica slab was melted at a temperature of 5000 K over the course of 100 ps, and then the resulting liquid silica was quenched to a temperature of 300 K for a period of 100 ps. These two simulations result in a bulk amorphous silica substrate. In order to generate the two free surfaces of the silica substrate, the dimensions of the simulation box were increased in the zdirection. Then, four different diameter pores (1.5, 2.0, 2.5 and

■

MODEL AND METHODOLOGY In this study, classical atomistic molecular dynamics (MD) simulations were used to study the structural and transport properties of ionic solutions as they are driven through charged silica nanopores via a flow field, which is used to model a 12299


Article

The Journal of Physical Chemistry C 3.0 nm; see Figure 1) were cut out of the silica substrate such that they run along the long axis (the z-axis) of the silica substrate. Then, the silica interfaces were annealed at 300 K in order to remove the majority of the high energy defects (i.e., silicon atoms with one or two oxygens bonded to them), which generally only takes approximately 0.3 ps. All of the simulations used to melt, quench, and anneal the silica substrates were conducted using the BKS potential.55 In an attempt to model the experimentally observed interfacial properties of amorphous silica as closely as possible, several steps are taken when finalizing the hydroxylated interface of the silica. First, interfacial silicon−oxygen bonds are chosen at random and broken in order to produce the correct number of undercoordinated silicon and oxygen atoms so that a surface density of silanol groups can be produced on the interfaces of the cylindrical pore, and the interfaces that exist at either end of the pore are in good agreement with experimental observations (∼4.6 nm−2).56 Then hydrogens and hydroxyl groups are placed on undercoordinated oxygens and silicons, respectively. This process results in each interface having isolated, geminal, and vicinal silanol groups present, at densities of approximately 1.2, 0.60, and 2.8 OH/nm−2, respectively, which is in accordance with the Zhuravlev model.57 A more detailed description of the methods used to generate these interfaces can be found in ref 34. After producing a realistic silica interface, then hydrogen atoms are removed from the isolated silanol groups in order to create a surface charge density of ∼0.9 e−/nm2, which is consistent with the charge density found experimentally at a pH ∼7.5 and 300 K.58,59 The resulting partial charges of the atoms in the deprotonated silanol groups were assigned using the methodology presented in previous publications.34,60 Figure 2 shows snapshots of the distribution of the oxygens that have been deprotonated on the pore surfaces of the four different diameter pores. These figures show that while the deprotonated oxygens are randomly distributed throughout the pore, there are areas of higher charge than others. The final dimensions of the silica pore are 6.10284 nm × 6.16904 nm × 7.46585 nm in the x-, y-, and z-dimensions, respectively. After the pore was hydroxylated and then charged such that it represents experimentally observed silica interfaces, water molecules were placed inside the pores such that the water density is consistent with the bulk equilibrium density of water. Also, preequilibrated water baths with dimensions of ∼6.1 nm × 6.1 nm × 2.5 nm were placed in contact with both sides of the pore. Finally, both water baths outside of the pores are randomly populated with ions (i.e., Na+ and Cl− or Ca2+ and Cl−) such that the resulting salt solution has a concentration of 0.5 M in the entire system. Table 1 contains a list of the numbers of ions and water molecules in each system. The final systems with both water baths and the pore are displayed in Figure 2. Periodic boundaries are applied in all three dimensions, and as a result the two water reservoirs on either side of the pore are actually one large reservoir. These systems were thermalized by first simulating the system using the NPT ensemble for 1 ns and then the NVT ensemble for another 1 ns. Then the systems were subjected to a “gravity” field (or flow field) instead of a pressure difference, as is done in most experiments. This gravity field, which is a standard option in the LAMMPS molecular dynamics simulation package, is used to generate a characteristic Poiseuille flow profile through the pores, and in doing so, the same acceleration is imparted on each atom within the fluid.

Table 1. Number of Ions in Solution, Number of Bare Oxygens on the Silica Pore Surface, and Number of Water Molecules in the Various Simulated Systems pore diam (nm)

system

cation

anion

bare oxygens

H2O

1.5

NaCl CaCl2 NaCl CaCl2 NaCl CaCl2 NaCl CaCl2

153 101 162 107 175 115 184 122

58 107 60 112 64 119 67 127

95 95 102 102 111 111 117 117

6710 6270 6965 6474 7405 6894 7689 7421

2.0 2.5 3.0

This results in each atom to which this gravity field is applied (in our case that is the water and ions in solution) receiving an added force which is proportional to its mass and therefore can be used to model a pressure field. The resulting flows from a flow field and a pressure field have the same parabolic velocity profile, so in this study a flow field is employed. The small differences between a flow field and a pressure-induced field are irrelevant for the questions which are being addressed in this paper. In these simulations, there is a flow field of 5 × 10−4 kcal/ (nm g), which is equivalent to a pressure drop of ∼150 atm. While MD simulations with smaller pressure fields would certainly be of interest, the present simulations allow for a very reasonable all-atom description of electrokinetic phenomena, and currently obtaining meaningful results from simulations conducted with smaller flow fields is not possible. Furthermore, the values of these external fields are in line with previous MD simulations.34 The production simulations were run for a period of 50 ns using the NVT ensemble, using the Nosé−Hoover thermostat.61 During these simulations, instead of freezing the entire silicon substrate, which can lead to undesired temperature gradients at the interfaces, the location of the silica pore is fixed by attaching three silicon atoms to their initial positions using springs with very large spring constants, and the remainder of the atoms in the pore were thermostated. The fixed silicons were selected such that they were in different parts of the silica pore and such that they were not near any of the pore interfaces. All of the simulations in this study were carried out with LAMMPS.62 The CHARMM force field was used to model the ions63 and silica,34,60 while the TIP3P model was used for water.64 This combination of force fields has been used widely in order to simulate the interactions of various ionic solutions and biological molecules in ionic solutions in contact with silica interfaces. The bulk properties of the TIP3P water model have been reported elsewhere,65−67 and we have included the plots of the number density of the oxygen atoms in the water molecules as a function of the z-dimension for each of our simulated systems in the Supporting Information. Also, the energetics of ion solvation have been studied for these force fields68 and we have the previously reported structural properties of ion hydration for the Cl− ion and compared it to that found from neutron scattering.69 The SHAKE algorithm70 was used to constrain the bond lengths and angles of the water molecules as well as all hydrogen containing bonds on the silica interface. All intermolecular interactions were cut off at 1.0 nm and the 12300


Article

The Journal of Physical Chemistry C long-range Coulombic interactions were computed using the PPPM algorithm.71

■

RESULTS The results presented in the following sections were determined by performing analysis of the final 40 ns of the simulation trajectories from the production simulations. This was done such that all measurements were done after the ions in solutions had plenty of time to counteract the surface charge of the silica pores. First, the structural properties of the ionic solutions outside and inside of the pore are presented. Then, the effect the transport has on the solvation of the ions and the pairing of ions will be presented. Finally, the transport properties of the ionic solutions through the nanopores are presented. Structural Properties. At the Pore Entrance/Exit. The average number of cations and anions per unit volume was calculated in the water reservoirs above and below (i.e., zcoordinate is less than the minimum z-coordinate of the pore or greater than the maximum z-coordinate of the pore), and then used to determine the charge density in these reservoirs. Figures 3 and 4 show the plots of the total charge density, the

Figure 4. Charge density (e−/nm3) as a function of distance from the pore entrance/exit into the ionic solvent reservoirs for the 0.5 M CaCl2 solutions in pores with diameters of 1.5, 2.0, 2.5, and 3.0 nm. The various curves represent the total charge density (black solid line), the charge density due to the Ca2+ ions (red dashed lines), and the charge density due to the Cl− ions (blue dotted lines). The orange dashed−dotted line represents a charge density of 0.0 for reference.

zero in each of the five defined regions of the system (as can be seen in the plots in the Supporting Information). In the NaCl systems, the charge densities in Figure 3 show that there is a broad peak of Na+ ions that is centered at approximately 0.4−0.5 nm from the silica substrate, where there is a peak present in each of the curves of the charge density due to the Na+ ions (red dashed lines). Then there is a layer of Cl− ions present after this layer of Na+ ions at a distance of ∼0.8 nm from the silica substrate, where again there is a peak in the anionic charge density in each plot. The charge density near the silica substrate within the ionic solution remains positive until it reaches zero at a distance of ∼1.5 nm from the surface of the silica pore. After this point where the net charge density becomes zero, the concentrations of the Na+ and Cl− ions become constant and their respective charges cancel out another, which results in an aqueous environment with a net charge of 0.0. The overall behavior of the charge densities is similar for each of the pore sizes studied. The only noticeable difference is that in some of the plots of the charge density due to the Na+ ions (specifically the d = 1.5 and 2.5 nm plots) there are two peaks in the charge density: one at ∼0.2 nm and a second at ∼0.5 nm away from the silica pore interface. In the other two plots, there is a less noticeable slope change in the charge density due to the Na+ ions at approximately 0.2 nm and a peak at ∼0.4−0.5 nm. We believe that these two different regions within the peak are representative of the separations of dehydrated and hydrated Na+ ions and the interface of the silica pores. The total charge density and the charge densities due to the contributions of the Ca2+ and Cl− ions for the 0.5 M CaCl2 systems is shown in Figure 4. In these systems, there is a peak in the charge density of the Ca2+ ions at a distance of ∼0.4 nm from the silica interface, which is the same distance at which there was a peak in the Na+ charge densities of the NaCl systems. The peaks in the CaCl2 systems are significantly larger and less broad than those observed for the charge density due to the Na+ cation in the NaCl systems. The peak number density of Ca2+ ions (the charge density divided by 2) is approximately 2.5 nm−3, which is almost twice the number

Figure 3. Charge density (e−/nm3) as a function of the distance from the pore entrance/exit into the ionic solvent reservoirs for the 0.5 M NaCl solutions in pores with diameters of 1.5, 2.0, 2.5, and 3.0 nm. The various curves represent the total charge density (black solid line), the charge density due to the Na+ ions (red dashed lines), and the charge density due to the Cl− ions (blue dotted lines). The orange dashed−dotted line represents a charge density of 0.0 for reference.

charge density due to the cations, and the charge density due to the anions in the reservoirs outside of the pores of different diameters for the 0.5 M NaCl and 0.5 M CaCl2 systems, respectively. In the plots, the average of the charge densities for the top and bottom reservoirs have been presented in order to improve upon the statistics of the results. In the Supporting Information, plots of the number of cations and anions in each reservoir, at each interface, and within the pore as a function of time are shown, and in these plots it can be seen that the system is at a steady state during the entire production simulation. The number of ions in each region remains more or less constant throughout the simulations. Also, the net charge (as defined as the charge due to cations in solution − the charge due to anions in solution − the charge due to the deprotonated silanols on the silica interface) is approximately 12301


Article

The Journal of Physical Chemistry C density of Na+ cations observed in the NaCl system (∼1.5 nm−3). The layer of Cl− ions that is observed in these systems occurs with a peak density at a distance of ∼0.6 nm from the silica interface. The peak number density of the Cl− ions in the CaCl2 systems is ∼1.5 nm−3, which is almost 4 times larger than the peak number density of the Cl− ions in the NaCl systems (∼0.4 nm−3). As a result of this distribution of ions within the ionic solution near the silica interface, the charge density of the solution is initially positive until a distance of ∼0.6 nm from the silica pore, at which point it becomes negative due to the large density of Cl− ions that are present. Then the net charge of the ionic solution remains negative until a distance of ∼1.0 nm from the silica interface, at which point it becomes 0.0, when the charge contributions due to the Ca2+ and Cl− ions become constant at values that result in the two contributions negating one another. Therefore, the effect of the ionic buildup near the charged silica substrate is observed a further 0.5 nm into the solution in the NaCl systems than in the CaCl2 systems. There is no noticeable difference in the behavior of ionic distributions outside of the pore in the CaCl2 systems for the different pore sizes. The fact that the sign of the net charge density changes as a function of distance from the silica interface in the CaCl2 systems is representative of a phenomenon called “charge inversion”. Charge inversion can be simply stated as occurring when interfacial charges attract counterions in excess of the nominal charge of the interface, and therefore an interface with a charge that has the opposite sign to the initial surface results. The underlying physics which is present within systems where charge inversion occurs is much more complex than is represented by the simple definition given here. For a more detailed description of the underlying physics that governs the phenomenon of charge inversion, we suggest reading refs 72−74. Inside the Pore. The charge density as a function of the distance from the center of the pore was calculated for the 0.5 M NaCl aqueous solution inside the four different sized pores and is presented in Figure 5. In all pores, the solution has a net

positive charge throughout the entire pore. At the pore entrance and exit, the surface charge effect of the pore extended to approximately ∼1.5 nm into the bulk solution. A similar effect is observed within the pores as well. For the pores with diameters less than 3.0 nm, Figure 5 shows that the there is no place within the pore where the charge decays to zero, and only in the very middle of the 3.0 nm diameter pore does the charge distribution within the NaCl solution approach 0.0. Therefore, there is no charge inversion observed within the pores. The Na+ ions form a layer near the silica pore surface that is approximately 0.5 nm in thickness. As the pore diameter increases, the density of Cl− ions inside the pore increases. The Cl− ions are distributed within the pore such that they have an increasing density as the center of the pore is approached. Near the silica pore surface, there is an ∼0.3 nm region in which there are no Cl− ions present in any of the pores. Figure 6 shows the charge density distribution of the 0.5 M CaCl2 solutions within the four different sized pores. Again, the

Figure 6. Charge density (e−/nm3) within the pore as a function of distance from the center of the pore (r = 0.0) for the 0.5 M CaCl2 solutions in pores with diameters of (a) 1.5, (b) 2.0, (c) 2.5, and (d) 3.0 nm. The various curves represent the total charge density (black solid line), the charge density due to the Na+ ions (red dashed lines), and the charge density due to the Cl− ions (blue dotted lines). The magenta dashed−dotted line represents a charge density of 0.0 for reference.

Ca2+ cations form a layer nearest the silica interface which is ∼0.5 nm in thickness. There is a layer of Cl− ions that forms just next to the Ca2+ layer, which is also ∼0.5 nm in size. Then, after the two layers, the charge densities due to the cations and due to the anions reach a constant value and remain constant throughout the remainder of the pore. The net charge within the 0.5 M CaCl2 solutions is largely positive nearest the silica interface, and then decreases until there is a layer of negatively charged solution which corresponds with the location of the layer of Cl− ions. Finally, the net charge increases to zero at a distance of ∼1.0 nm from the silica pore interface. This behavior is nearly identical to that observed at the pore entrance and exit. The charge inversion is clearly observed for the pores with diameters greater than or equal to 2.0 nm, when there is enough space within the pore for the double layers to form. Hydration of Water Molecules and Ions. In order to study how the hydration of the water molecules and ions changed as they were driven through the pore, we first measured the radial distribution functions (rdf’s) of the relative

Figure 5. Charge density (e−/nm3) within the pore as a function of distance from the center of the pore (r = 0.0) for the 0.5 M NaCl solutions in pores with diameters of (a) 1.5, (b) 2.0, (c) 2.5, and (d) 3.0 nm. The various curves represent the total charge density (black solid line), the charge density due to the Na+ ions (red dashed lines), and the charge density due to the Cl− ions (blue dotted lines). The magenta dashed−dotted line represents a charge density of 0.0 for reference. 12302


Article

The Journal of Physical Chemistry C

Table 2. Average Number of First Neighbor Water Molecules Surrounding Free and Bound Water Molecules, Na+ Ions, and Cl− Ions in the Different Diameter Pore Systems pore diameter (nm) location free

reservoir

entrance/exit

inside pore

bound

entrance/exit

inside pore

type H2O Na+ Cl− H2O Na+ Cl− H2O Na+ Cl− H2O Na+ Cl− H2O Na+ Cl−

1.5 5.3 5.7 8.1 5.3 5.7 8.2 5.1 5.0 − 4.1 3.3 6.7 3.7 3.1 −

2.0

± ± ± ± ± ± ± ±

0.8 0.6 1.0 0.8 0.7 1.1 1.2 1.7

± ± ± ± ±

1.1 1.2 1.8 1.2 1.2

components so that we could determine the first neighbor distance which is characterized by the first minimum in the plot of the radial distribution function. Then we used the first neighbor distance to count the number of first neighbors each atom of a certain type had and averaged over time for atoms found in three different portions of the system: in the reservoirs outside of the pore, in the transitional region of the system where the water molecules and ions are entering or exiting the pore, and inside of the pore. Also, we have differentiated water molecules or ions that are bound to the silica interface (i.e., within the first neighbor distance of the oxygen or hydrogen of a silanol) from those that are not bound, which we refer to as “free” throughout the remainder of the paper. For the 0.5 M NaCl solution, the average number of first neighbor water molecules surrounding free and bound water molecules, Na+ ions, and Cl− ions in the three different regions of the different sized pores is shown in Table 2. The average hydration number of water molecules does not change as the diameter of the pore increases. The average hydration number of free water molecules is 5.3 ± 0.8, whether the water molecule is in the water reservoirs, at the pore exit/entrance, or inside the pore. Meanwhile, the average hydration number of bound water molecules is 4.1 ± 1.1 at the silica interfaces at the pore entrance and exit, and 3.8 ± 1.2 within the pore. Figure 7 shows histograms of the number of first neighbor water molecules around a free or bound water molecule within the 0.5 M NaCl solution in the reservoirs, at the pore entrance and exit, and inside the pores for the pore diameters of 1.5, 2.0, 2.5, and 3.0 nm. When considering the free water molecules, we observe that in all cases the most probable number of first neighbor water molecules is 3. However, the percentage of water molecules with 3 first neighbor water molecules decreases as one moves from the water reservoir to the entrance/exit of the pore to inside the pore. Similarly, the number of water molecules with less than 3 first neighbor water molecules is larger inside of the pore than it is at the pore entrance/exit, which is larger than that for the water molecules in the reservoir. So while the average number of first neighbor water molecules is nearly the same throughout the system, it is clear that there is some dehydration observed of the water molecules as a result of the confinement of the pores. Additionally, there seems to be more dehydration of water molecules (as seen by

5.3 5.6 8.1 5.3 5.6 8.2 5.1 5.7 8.3 4.1 3.4 7.5 3.7 2.8 5.8

± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

2.5 0.8 0.7 1.0 0.8 0.8 1.0 1.0 0.6 1.1 1.1 1.1 1.4 1.2 1.2 1.8

5.3 5.7 8.2 5.3 5.7 8.2 5.2 5.7 8.3 4.1 3.1 7.2 3.8 2.9 7.6

± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

3.0 0.8 0.7 1.0 0.8 0.6 1.0 1.0 0.6 1.1 1.1 1.2 1.4 1.2 1.2 1.3

5.2 5.6 8.0 5.2 5.6 8.1 5.1 5.6 8.2 4.1 3.3 7.3 3.8 3.2 6.6

± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.8 0.7 1.0 0.8 0.8 1.0 0.9 0.7 1.1 1.1 1.2 1.4 1.2 1.2 1.8

Figure 7. Distribution of the number of nearest neighbor water molecules around a water molecule within the 0.5 M NaCl solutions in the reservoirs (black empty bars), at the pore entrance/exit (red bars with cross hatching), and inside the pores (blue filled bars) for pores with diameters of 1.5, 2.0, 2.5, and 3.0 nm. Separate distributions are presented for water molecules that are bound to the silica substrates and those that are free in solution.

the slightly larger percentages of water molecules with less than 3 first neighbor water molecules) in the smaller diameter pores. Meanwhile, the histograms corresponding to the hydration of the bound water molecules in Figure 7 show that bound water molecules also have a general trend of being dehydrated within the pore as compared to at the entrance/exit of the pore. In all pore diameters, the most probable number of first neighbor water molecules around a bound water molecule at the entrance/exit of the pore is 2.5. Inside of the pore the most probable number of first neighbor water molecules around a bound water is 2.0 for all pore diameters. Similarly, there is a larger percentage of bound water molecules that have less than 2.0 first neighbor water molecules inside of the pore than at the pore entrance/exit. In Table 2, the average hydration number of free Na+ cations is the same (∼5.7) outside, at the entrance/exit, and inside of the pore for the pore diameters larger than 1.5 nm. However, for the pore with a diameter of 1.5 nm, inside the pore the 12303


Article

The Journal of Physical Chemistry C average hydration number (∼5.0) is significantly smaller than in the reservoir or at the pore entrance/exit, which demonstrates some dehydration of the Na+ ion. This dehydration of the free Na+ ion in the pore with a 1.5 nm diameter is reinforced by the histogram of the number of first neighbor water molecules around a Na+ ion within the 0.5 M NaCl solutions shown in Figure 8.

Figure 9. Distribution of the number of water molecules in the second hydration shell of a Na+ ion within the 0.5 M NaCl solutions in the reservoirs (black empty bars), at the pore entrance/exit (red bars with cross hatching), and inside the pores (blue filled bars) for pores with diameters of 1.5, 2.0, 2.5, and 3.0 nm. Separate distributions are presented for Na+ ions that are bound to the silica substrates and those that are free in solution.

Also, as was the case with the number of first neighbor water molecules, the pore size effect is noticeable for the distribution of the number of water molecules within the first two hydration shells for pore diameters less than 3.0 nm. The distributions of the number of first neighbor water molecules and the number of water molecules in the first two hydration shells around a Cl− ion within the 0.5 M NaCl solutions are shown in Figures 10 and 11. For the free Cl− ions, there is no change in the number of first neighbor water molecules as they go through the pore observed in any of the diameter pores and only a slight decrease in the number of water molecules in the first two hydration shells. There is a significant decrease in the number of water molecules in the first and first two solvation shells of bound Cl− ions in the d =

Figure 8. Distribution of the number of nearest neighbor water molecules around a Na+ ion within the 0.5 M NaCl solutions in the reservoirs (black empty bars), at the pore entrance/exit (red bars with cross hatching), and inside the pores (blue filled bars) for pores with diameters of 1.5, 2.0, 2.5, and 3.0 nm. Separate distributions are presented for Na+ ions that are bound to the silica substrates and those that are free in solution.

The bound Na+ ions, which are directly interacting with the silica substrate, have approximately only half the number of first neighbor water molecules as compared to the free cations (see Table 2). Again, a dehydration effect is observed for the Na+ ions when comparing the hydration of bound cations at the entrance/exit of the pore to those that are bound to the silica inside of the pore. From the histograms shown in Figure 8, the distribution of the number of first neighbor water molecules around a bound Na+ ion is different for those water molecules at the pore entrance/exit from those within the pore for pores with diameters up to 2.5 nm in size, and then in the 3.0 nm pore the distributions are almost identical. So there is a more significant dependence on the pore size on the hydration of bound Na+ ions than free Na+ ions in our systems. We also measured the number of water molecules that were within the first two hydration shells of the Na+ ions in the three different regions of our simulation systems for the four different pore sizes, and the distributions are shown in Figure 9. We did this by calculating the number of water molecules that are first neighbors of those water molecules that are first neighbors of the Na+ ions. For the free Na+ ions, there is a clear trend that the number of water molecules within the first two hydration shells is smaller at the pore entrance/exit than in the reservoir, and smaller yet inside of the pore. As was the case when considering the number of first neighbor water molecules around a Na+ ion, the most significant effect is observed in the 1.5 nm diameter pore, while the distributions look quite similar for the 2.0, 2.5, and 3.0 nm diameter pores. For the bound Na+ ions, again there is a noticeable decrease in the number of water molecules within the first two hydration shells of cations within the pore when compared to those at the pore entrance/exit.

Figure 10. Distribution of the number of nearest neighbor water molecules around a Cl− ion within the 0.5 M NaCl solutions in the reservoirs (black empty bars), at the pore entrance/exit (red bars with cross hatching), and inside the pores (blue filled bars) for pores with diameters of 1.5, 2.0, 2.5, and 3.0 nm. Separate distributions are presented for Cl− ions that are bound to the silica substrates and those that are free in solution. 12304


Article


Figures 12 and 13, respectively. For the pores with diameter less than 2.5 nm, the largest percentage of free Ca2+ ions have

Figure 11. Distribution of the number of water molecules in the second hydration shell of a Cl− ion within the 0.5 M NaCl solutions in the reservoirs (black empty bars), at the pore entrance/exit (red bars with cross hatching), and inside the pores (blue filled bars) for pores with diameters of 1.5, 2.0, 2.5, and 3.0 nm. Separate distributions are presented for Cl− ions that are bound to the silica substrates and those that are free in solution.

Figure 12. Distribution of the number of nearest neighbor water molecules around a Ca2+ ion within the 0.5 M CaCl2 solution in the reservoirs (black empty bars), at the pore entrance/exit (red bars with cross hatching), and inside the pores (blue filled bars) for pores with diameters of 1.5, 2.0, 2.5, and 3.0 nm. Separate distributions are presented for Ca2+ ions that are bound to the silica substrates and those that are free in solution.

2.0 and 3.0 nm pores, while there is no significant difference observed in the d = 2.5 nm pore. For the 0.5 M CaCl2 solutions, Table 3 shows the average number of first neighbor water molecules surrounding the free and bound water molecules, Ca2+ ions, and Cl− ions in the four different pore sizes. The values in Table 3 show that there is no evidence of dehydration of water as it travels through the pore, and this was reinforced by the distributions of the number of first neighbor water molecules around a water molecule throughout the system (not shown). Also, we observe that in the d = 1.5 and 2.0 nm pores the average number of first neighbor water molecules around both ions actually increases within the pore as compared to those ions in the reservoir or at the pore entrance/exit. The distributions of the number of nearest neighbor water molecules and the number of water molecules within the first two hydration shells for the 0.5 M CaCl2 systems are shown in

more than 7 water molecules in the first hydration shell, but at the entrance and exit to the pore and in the reservoirs on either side of the pore, the largest percentage of free Ca2+ ions have 4 water molecules within their first hydration shells. Then in the pores with diameters of 2.5 and 3.0 nm, the distributions inside and outside of the pore are nearly identical. The bound cations also show a similar trend where the bound cations inside the pore have more water molecules within their first hydration shell, as compared to those bound to the pore entrance/exit. Also, as the pore gets larger again the distributions become more similar. While there is an increase in the number of water molecules in the first hydration shells of free Ca2+ ions inside the pore as compared to those not in the pore, there is a significant

Table 3. Average Number of First Neighbor Water Molecules Surrounding Free and Bound Water Molecules, Ca2+ Ions, and Cl− Ions in the Different Diameter Pore Systems pore diameter (nm) location free

reservoir

entrance/exit

inside pore

bound

entrance/exit

inside pore

type H2O Ca2+ Cl− H2O Ca2+ Cl− H2O Ca2+ Cl− H2O Ca2+ Cl− H2O Ca2+ Cl−

1.5 4.6 5.0 8.2 4.6 5.2 8.3 4.5 7.0 8.8 3.5 4.6 7.9 3.3 5.3 −

± ± ± ± ± ± ± ± ± ± ± ± ± ±

2.0 0.8 1.4 1.0 0.8 1.4 1.1 1.1 0.6 1.1 1.0 1.0 1.3 1.2 0.8

12305

4.6 5.3 8.1 4.5 5.2 8.3 4.4 6.4 8.5 3.5 4.6 8.1 3.3 4.5 8.1

± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

2.5 0.8 1.5 1.0 0.8 1.4 1.1 0.9 1.3 1.1 1.0 0.9 1.2 1.1 1.3 1.3

4.5 5.2 8.1 4.5 5.3 8.3 4.4 5.5 8.2 3.5 4.5 7.8 3.3 4.8 8.1

± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

3.0 0.8 1.4 1.0 0.8 1.5 1.1 0.9 1.5 1.0 1.0 1.0 1.3 1.1 1.0 1.2

4.5 5.1 8.1 4.5 5.1 8.3 4.5 5.5 8.3 3.5 4.5 7.8 3.3 4.9 7.7

± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.8 1.4 1.0 0.8 1.4 1.1 0.9 1.5 1.1 1.0 1.1 1.3 1.0 1.0 1.2


Article


Table 4. Distribution of the Number of Cl− Anions Bound to a Na+ Cation within the 0.5 M NaCl Systems in the Pores of Different Diameters free Na+ pore diam (nm)

location

0

1

0

1

1.5

reservoir entrance/exit inside pore reservoir entrance/exit inside pore reservoir entrance/exit inside pore reservoir entrance/exit inside pore

0.88 0.90 1.00 0.88 0.91 0.97 0.88 0.90 0.93 0.87 0.89 0.91

0.12 0.10 0.00 0.12 0.09 0.03 0.12 0.10 0.07 0.13 0.11 0.09

− 0.95 1.00 − 0.93 0.99 − 0.88 0.98 − 0.93 0.87

− 0.05 0.00 − 0.07 0.01 − 0.12 0.02 − 0.07 0.13

2.0

2.5

3.0

Figure 13. Distribution of the number of water molecules in the second hydration shell of a Ca2+ ion within the 0.5 M CaCl2 solution in the reservoirs (black empty bars), at the pore entrance/exit (red bars with cross hatching), and inside the pores (blue filled bars) for pores with diameters of 1.5, 2.0, 2.5, and 3.0 nm. Separate distributions are presented for Ca2+ ions that are bound to the silica substrates and those that are free in solution.

bound Na+

Table 5. Distribution of the Number of Cl− Anions Bound to a Ca2+ Cation within the 0.5 M CaCl2 Systems in the Pores of Different Diameters free Ca2+ pore diam (nm)

decrease in the number of water molecules within the first two hydration shells of the cations inside the pore as compared to those not inside of the pore for pores with diameters of 1.5 and 2.0 nm where the whole distribution is shifted to a smaller number of water molecules. In the pores with diameter 2.5 and 3.0 nm, the distributions of the number of water molecules in the first two hydration shells become more broad at the pore entrance/exit and inside of the pore (ranging from 12.5 to 27.5 water molecules) whereas in the bulk it only ranges from 20 to 27.5 water molecules. For the bound Ca2+ ions, at the pore entrance/exit the number of water molecules in the first two hydration shells is consistently between 12.5 and 22.5 water molecules for all pore diameters. Inside the pore, the distribution ranges from 5 to 20 water molecules for pore diameters of 1.5, 2.0, and 2.5 nm. Only in the 3.0 nm diameter pore are the distributions of the number of water molecules in the first two hydration shells for the bound cations at the pore entrance/exit and inside of the pore nearly identical. We have measured the distribution of the number of water molecules in the first and the first and second hydration shells of the Cl− ions in the 0.5 M CaCl2 solutions. These distributions (not shown) show no dehydration effect of the free or bound anions in the various regions of the pore system or with pore diameter. Ion Pairing. We have also determined the amount of ion pairing in the various systems we have simulated. Table 4 shows the fraction of free and bound Na+ ions that have 0 or 1 Cl− ions within a first neighbor distance. The values show that for pores with diameters of 1.5 and 2.0 nm the ion pairing of both free and bound Na+ ions is disrupted by the confinement within the pores, as the fraction of cations paired with 1 Cl− ion decreases as the ion moves from the reservoir to the pore entrance/exit and then further decreases when the ions are inside the pore. There is not as drastic of an effect for the 2.5 and 3.0 nm diameter pores. Table 5 shows the fraction of free and bound Ca2+ ions that have 0, 1, or 2 Cl− ions within a first neighbor distance of it. The values show that for pores with diameters of 1.5, 2.0, and 2.5 nm the ion pairing of the free Ca2+ ions is disrupted by the

1.5

2.0

2.5

3.0

bound Ca2+

location

0

1

2

0

1

2

reservoir entrance/ exit inside pore reservoir entrance/ exit inside pore reservoir entrance/ exit inside pore reservoir entrance/ exit inside pore

0.24 0.27

0.21 0.22

0.55 0.51

− 0.58

− 0.40

− 0.02

1.00 0.31 0.27

0.00 0.20 0.22

0.00 0.49 0.51

0.90 − 0.48

0.10 − 0.50

0.00 − 0.02

0.67 0.30 0.33

0.17 0.21 0.21

0.16 0.49 0.46

0.78 − 0.51

0.22 − 0.47

0.00 − 0.02

0.43 0.24 0.26

0.08 0.22 0.24

0.49 0.54 0.50

0.83 − 0.54

0.17 − 0.46

0.00 − 0.00

0.38

0.19

0.43

0.78

0.22

0.00

confinement within the pores. The bound cations have their ion pairing significantly affected within the pores as compared to the cations at the pore entrance and exit. This is probably due to the anions in the diffuse layer being driven through the nanopore and therefore disrupting the pairing. Also, the bound ions, which are acting to counter the surface charge, generally have one less anion that pairs with it as a result of its interaction with the silica interface. Transport Properties. The average velocity in the zdimension of the water molecules was calculated as a function of the distance from the center of the pore. The distance was divided into 0.1 nm increments, and the z-component of the velocity of the water molecules in a given bin is then averaged over the number of water molecules present in the bin and over time. The resulting velocity profiles for the water molecules in the 0.5 M NaCl systems are shown in Figure 14. As shown in Figure 14, the peak velocity in the central region of the pore increases as the pore diameter increases. The shape of the velocity profile changes from a more plateau-like curve for the pore with diameter of 1.5 nm to a deformed parabola at d = 2.0 12306


Article


Figure 14. Velocity profiles of the water molecules in the 0.5 M NaCl solutions within pores with diameters of 1.5 nm (black ○), 2.0 nm (red □), 2.5 nm (green ◇), and 3.0 nm (blue ×). The center of the pore is designated as r = 0 nm.

Figure 15. Velocity profiles of the water molecules in the 0.5 M CaCl2 solutions within pores with diameters of 1.5 nm (black ○), 2.0 nm (red □), 2.5 nm (green ◇), and 3.0 nm (blue ×). The center of the pore is designated as r = 0 nm.

nm and then finally a more well-defined parabola for d = 2.5 and 3.0 nm. Plots of the velocity profiles for the Na+ and Cl− ions (which are presented in the Supporting Information) show the same general behavior, in that the maximum velocities are observed in the center of the pore with very little motion of ions within ∼0.8 nm of the pore interface. Also, it is worth noting that in all cases the ions are observed to move in the same direction as the water, which in turn is the same direction as the applied flow field. The same parabola-shaped velocity profile of the water molecules, Na+ ions, and Cl− ions was observed in simulations of streaming current flow of 0.4 M NaCl solutions inside charged silica slit pores.34 We have also calculated the flux values for the water molecules, Na+ ions, and Cl− ions for the 0.5 M NaCl systems, and present those values in Table 6. The values show the same

in the NaCl system. However, for the 2.5 and 3.0 nm diameter pores, the peak velocities have approximately the same values. In both the NaCl and CaCl2 systems, the velocities of the water molecules are near zero for molecules within ∼0.4 nm of the pore wall, which is due to the tightly bound cations to the negatively charged silica interface. This is consistent with the width of the first peaks in the charge density as shown in Figures 5 and 6, for the NaCl and CaCl2 systems, respectively. Then, in the CaCl2 systems, there is a region in the velocity profiles of the water molecules that is ∼0.4−0.6 nm in thickness where the velocities are slightly larger but not as large as those found in the NaCl systems within the same pore size. This is due to the fact that in the CaCl2 system there is charge inversion, while there is no such effect in the NaCl systems, and therefore a significant diffuse layer of Cl−, whose thickness is consistent with the size of this region in the velocity profiles as can be seen in Figures 6 and 15. The velocity profiles for the cations and anions in the 0.5 M CaCl2 systems share the same general shape as is observed in the velocity profiles of the water molecules in that they have their maximum velocity in the center of the pore and have a generally parabolic shape (see Supporting Information). Also, as described above, the velocity profiles of the Ca2+ ions demonstrate very small velocities within 1.0 nm of the pore walls, which corresponds to the thickness of the region over which we observe charge inversion in our systems. Meanwhile, the Cl− ions show very small velocities when found within 0.7 nm of the pore wall, and then their velocities begin to increase, and therefore there is further evidence that in the diffuse layer, where the Cl− ions bind to the layer of Ca2+ ions that are tightly bound to the negatively charged silica interface, the anions are able to be transported through the pore. Table 7 shows the flux of the water molecules, Ca2+ ions, and − Cl ions that have been determined from our simulations as a function of pore diameter. As was the case in the NaCl systems, we observe the same trend as would be expected from the velocity profiles, which is that as pore diameter increases so does the flux of each component of our ionic solution. In the case of the CaCl2 solutions, we observe that as the pore diameter increases the flux of the Cl− anions increases more rapidly than does the flux of the Ca2+ cations, which results in an increasingly negative flux in the net charge. When comparing the fluxes of the individual components of the ionic solutions

Table 6. Flux Values for the Water Molecules, Na+ Ions, and Cl− Ions through the Pores with Diameter d = 1.5, 2.0, 2.5, and 3.0 nm pore diam (nm)

water

Na+

Cl−

net charge

1.5 2.0 2.5 3.0

5.52 11.81 19.09 31.51

0.01 0.05 0.12 0.32

0.00 0.02 0.08 0.21

0.01 0.03 0.04 0.11

trends as would be expected from the relative comparison of the velocity profiles. As the pore diameter increases, the flux of each component increases. The flux of cations increases more rapidly than the flux of anions, and therefore as the pore size increases the net flux of charge across the pore also increases. Figure 15 shows the velocity profiles in the z-direction of the water molecules in the 0.5 M CaCl2 systems as a function of the distance from the center of the pore. As in the NaCl systems, the peak velocity is generally found in the center of the pore and the value of the peak velocity increases as the pore diameter increases. In the CaCl2 systems, a parabolic velocity profile is observed for all of the pores except the 1.5 nm diameter pore. In this smallest pore, the maximum water velocity is observed nearest the pore walls and in the center the water has an average velocity in the z-dimension of ∼0. In the d = 1.5 and 2.0 nm pores, the peak velocity values of the water molecules in the CaCl2 system are smaller than those measured 12307



■

Table 7. Flux Values for the Water Molecules, Ca2+ Ions, and Cl− Ions through the Pores with Diameter d = 1.5, 2.0, 2.5, and 3.0 nm pore diam (nm)

water

Ca2+

Cl−

net charge

1.5 2.0 2.5 3.0

2.97 9.10 11.68 22.36

0.00 0.01 0.04 0.09

0.00 0.01 0.14 0.25

0.00 0.01 −0.06 −0.07

1 ΔtLz

N

∑ qi(zi(t + Δt ) − zi(t )) i=0

DISCUSSION

Ion transport through cylindrical charged silica nanopores has been studied using all-atom molecular dynamics simulations with a flow field applied as a mimic for a pressure-driving force. The transport of both 0.5 M solutions of NaCl and CaCl2 through cylindrical pores of diameter 1.5, 2.0, 2.5, and 3.0 nm has been studied. The interfacial properties of the pores were created such that they represent the silica interface observed in experimental systems, using a method that had been utilized in several previous studies.34,44,53 The simulations that were carried out produced trajectories that were used to study the structural and transport properties of the ionic solutions as they travel through the different diameter pores. We measured the distribution of ions within the aqueous solutions as a function of distance from the silica interfaces at the pore entrance/exit and within the pores. From these measurements, we observe a small amount of charge inversion for all simulated systems of 0.5 M CaCl2 solutions. There is no such charge inversion observed for any of the simulated systems of the 0.5 M NaCl solutions. These findings are consistent with previous simulation studies which have observed charge inversion for 0.2,34 0.5,44 and 1.0 M34 CaCl2 solutions, but no charge inversion for 0.434 and 0.5 M44 NaCl solutions. Experimental results have previously shown that charge inversion does not occur at similar concentrations of NaCl but that, for concentrations larger than 0.4 M, CaCl2 solutions do show charge inversion.59 The effect of the charged silica interface is observed to affect the structure of the NaCl solution up to ∼1.5 nm from the interface, while for the CaCl2 solution it is observed to have an effect up to ∼1.0 nm from the interface. This is consistent with the findings of previous simulation studies of ionic solutions near charged silica interfaces.33,34,44 We have measured the hydration of free and bound water molecules, cations (Na+ or Ca2+), and the Cl− anions when they are in a reservoir, within the interfacial region of the pore entrance/exit, and inside the pore. For the free species, we observed slight dehydration of the water molecules in the two smallest diameter pores that we studied (1.5 and 2.0 nm). In the 0.5 M NaCl systems, we observed significant dehydration of first neighbor water molecules around the Na+ ions in the 1.5 nm diameter pore, and very little effect of the first neighbor hydration shell of the Na+ ions in the larger pores. When investigating the effect to the first two hydration shells, we observe that there is significant dehydration of the first two hydration shells around Na+ ions in all pore diameters. There is no noticeable dehydration of first neighbor water molecules around Cl− ions and only minor dehydration of the water molecules that make up the first two hydration shells of the Cl− ions. In the 0.5 M CaCl2 solutions, we observe that for the smallest pores (1.5 and 2.0 nm diameter) the free Ca2+ ions that enter have a larger number of water molecules in the first hydration but a significantly smaller number of water molecules within the first two hydration shells than found in the 2.5 and 3.0 nm diameter pores. In the 2.5 and 3.0 nm diameter pores, there is no dehydration effect to the first hydration shell but there is a noticeable dehydration of the water molecules in the first two hydration shells. Again, there is no noticeable dehydration of the free Cl− ions in the CaCl2 solutions either. We also measured the amount of ion pairing between cations and anions in the reservoirs, the interfacial regions, and the

between the NaCl and the CaCl2 solutions, we see that the net fluxes of water molecules and anions are almost identical, whereas the net flux of the cations is approximately 4 times larger in the NaCl solutions as compared to that in the CaCl2 solutions. Finally, we have calculated the instantaneous ionic current I(t) using the expression I (t ) =

Article

(1)

where Lz is the length of the pore in the z-dimension, Δt is the difference in time between snapshots in the trajectory, qi is the charge of the ion, and zi(t) and zi(t + Δt) are the z-coordinates of the ion at times t and t + Δt, respectively. The average ionic current is then calculated by averaging the value determined using eq 1 over the time of the production simulation. Figure 16 shows a plot of the average net ionic current and the average current that is due to the flow of cations and anions

Figure 16. Plot of average ionic current as a function of the pore diameter. The total average current (black □), as well as the average ionic current contributed by cations (red ○) and anions (blue ◇), is plotted for the (a) NaCl aqueous solution and (b) CaCl2 aqueous solution systems.

in the NaCl and CaCl2 solutions as a function of pore diameter. In the 0.5 M NaCl systems, the ionic current generally increases as the pore diameter increases, and this is due to the fact that the magnitude of the current generated by the flow of Na+ ions increases more rapidly than the magnitude of the current generated by the flow of Cl− ions. In the 0.5 M CaCl2 systems, we observe a trend of increasingly negative average net current as the pore diameter increases. This is due to the fact that in these systems the magnitude of the average current due to the flow of the Cl− ions increases more rapidly with pore diameter than does the average current due to the flow of the Ca2+ ions. 12308


Article


have used to study ion transport through silica nanopores as well.12,34,44 Currently, we are carrying out simulations of model silica nanopores functionalized with amino acids, which are commonly found in ion channels, to understand ion transport through these biological nanopores from a fundamental level. In doing so, we plan to investigate the origin of the different behavior that is observed in simulation studies of ion conduction through solid state and biological nanopores (one example we discussed in our previous paper44).

pore for the different pore diameters. In the 0.5 M NaCl and CaCl2 systems, there is a disruption to the ion pairing within all of the different diameter pores as compared to the bulk, and the disruption of the ion pairing seems to begin as it transitions from the bulk into the interfacial region. As the pore diameter increases, the amount of ion pairing within the pores increases. In the case of the CaCl2 systems, the number of Ca2+ ions that is paired to one Cl− ion is nearly the same in the pore and in the bulk for the largest diameter pore (3.0 nm), but the number of Ca2+ ions that are paired with two Cl− ions is still 20% less than in the bulk. In the NaCl system, the number of Na+ ions paired with a Cl− ion is still 30% less inside the 3.0 nm diameter pore than is found in the reservoirs. The velocity profiles of the water molecules and ions in the NaCl and CaCl2 solutions were measured within the different diameter pores. In the 0.5 M NaCl systems, the velocity profiles of the water molecules transition from a pluglike flow to the expected parabolic shape as the pore diameter increases. The magnitude of the velocities of each components increases with the pore diameter, and the maximum velocity of each component in the center of the pores increases with pore diameter. In the 0.5 M CaCl2 systems, the velocity profiles of the water molecules show the expected parabolic shape for all pore diameters except for the smallest pore diameter (d = 1.5 nm). An increase in the magnitudes and the maximum velocity of the various components is also observed in the 0.5 M CaCl2 systems as the pore diameter increases. The maximum velocities in the NaCl and CaCl2 systems are fairly consistent, which is different from what is observed in electroosmotic flow where the external field causes ions to move in opposite directions34,44 but is consistent with what was observed in charged silica slit pores previously.34 Finally, we have investigated the transport properties of the ionic solutions through the different diameter pores. The transport of ions in the two different ionic solutions, 0.5 M NaCl and 0.5 M CaCl2, is significantly different. In the 0.5 M NaCl systems, we observed that larger numbers of Na+ ions are transported through the pores than Cl− ions, and the difference increases as the pore diameter increases. Therefore, as the pore diameter increases the resulting charge flux through the pore is increasingly positive. In the 0.5 M CaCl2 systems, the Cl− ions are transported through the pores in larger numbers than the Ca2+ ions. The difference between the number of anions and the number of cations transported through the pores increases as the diameter of the pore increases, and therefore the net charge flux through the pore becomes increasingly negative as the pore diameter increases. Similar behavior was observed in a previous study of transport of NaCl and CaCl2 through charged silica slit pores with a diameter of 7.5 nm.34 Therefore, the pressure-driven flow of these ionic solutions may not completely isolate either of the cationic or anionic species; there is potential to dilute the concentrations of Cl− ions in NaCl solutions and Ca2+ ions in the CaCl2 solutions. Similar properties are of importance when considering the mechanisms by which biological ion channels function. Recent review articles have summarized the various methods that have been used and the contributions that have been made applying computer simulations to studying a variety of these ion channels.75,76 Generally, those groups that use molecular dynamics simulations to study the conduction of ions through biological ion channels either use free-energy methods like umbrella sampling in order to drive the ions through the ion channel, or they will use a small external electric field, which we

■

ASSOCIATED CONTENT

S Supporting Information *

Plots of the (a) velocity profiles of the cations and anions for each pore diameter in the 0.5 M NaCl and 0.5 CaCl2 systems; (b) number density of the oxygen in the water molecules as a function of the z-dimension for each system; (c) number of cations and anions in the bottom reservoir (reservoir 1), at the bottom interface (interface 1), within the pore, at the top interface (interface 2), and within the top reservoir (reservoir 2) as a function of time; and (d) net charge (defined as the charge of the cations in solution − the charge of the anions in solution − the charge of the deprotonated silanols at the interface of the silica pore) as a function of time in the bottom reservoir (reservoir 1), at the bottom interface (interface 1), within the pore, at the top interface (interface 2), and within the top reservoir (reservoir 2). The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/jp5129639.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS C.D.L. recognizes the support and computational time provided by the Barcelona Supercomputing CentreCentro Nacional de Supercomputation. Via our membership in the United Kingdom’s HPC Materials Chemistry Consortium, which is funded by the Office of Science and Technology through EPSRC’s High End Computing Programme, Grant EP/L000202, we utilized the facilities of ARCHER, the United Kingdom’s national high performance computing service, for aspects of this work. C.D.L. and N.H. recognize the support of the Engineering and Physical Sciences Research Council for funding N.H.’s studentship during which this work was carried out.

■

REFERENCES

(1) Plecis, A.; Schoch, R. B.; Renaud, P. Ionic Transport Phenomena in Nanofluidics: Experimental and Theoretical Study of the ExclusionEnrichment Effect on a Chip. Nano Lett. 2005, 5, 1147−1155. (2) Noy, A.; Park, H. G.; Fornasiero, F.; Holt, J. K.; Grigoropoulos, C. P.; Bakajin, O. Nanofluidics in Carbon Nanotubes. Nano Today 2007, 2, 22−29. (3) Schoch, R. B.; Han, J. Y.; Renaud, P. Transport Phenomena in Nanofluidics. Rev. Mod. Phys. 2008, 80, 839−883. (4) Mortensen, N. A.; Xiao, S. S.; Pedersen, J. Liquid-Infiltrated Photonic Crystals: Enhanced Light-Matter Interactions for Lab-On-AChip Applications. Microfluid. Nanofluid. 2008, 4, 117−127. (5) Craighead, H. Future Lab-On-A-Chip Technologies for Interrogating Individual Molecules. Nature 2006, 442, 387−393.

12309


Article

The Journal of Physical Chemistry C (6) Abgrall, P.; Nguyen, N. T. Nanofluidic Devices and Their Applications. Anal. Chem. 2008, 80, 2326−2341. (7) Napoli, M.; Eijkel, J. C. T.; Pennathur, S. Nanofluidic Technology for Biomolecule Applications: A Critical Review. Lab Chip 2010, 10, 957−985. (8) Li, Y.; Xu, J.; Li, D. Molecular dynamics simulation of nanoscale liquid flows. Microfluid. Nanofluid. 2010, 9, 1011−1031. (9) Daiguji, H. Ion Transport in Nanofluidic Channels. Chem. Soc. Rev. 2010, 39, 901−911. (10) Modi, N.; Winterhalter, M.; Kleinekathofer, U. Computational Modeling of Ion Transport through Nanopores. Nanoscale 2012, 4, 6166−6180. (11) Zharov, I.; Khabibullin, A. Nanoporous Films and Membranes with Controlled Ionic and Molecular Transport. Acc. Chem. Res. 2014, 47, 440−449. (12) Leung, K.; Rempe, S. B.; Lorenz, C. D. Salt Permeation and Exclusion in Hydroxylated and Functionalized Silica Pores. Phys. Rev. Lett. 2006, 96, 095504. (13) Corry, B. Designing Carbon Nanotube Membranes for Efficient Water Desalination. J. Phys. Chem. B 2007, 112, 1427−1434. (14) Fornasiero, F.; Park, H. G.; Holt, J. K.; Stadermann, M.; Grigoopoulos, C. P.; Noy, A.; Bakajin, O. Ion Exclusion by Sub-2-nm Carbon Nanotubes Pores. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 17250−17255. (15) Elma, M.; Yacou, C.; Wang, D. K.; Smart, S.; Diniz da Costa, J. C. Microporous Silica Based Membranes for Desalination. Water 2012, 4, 629−649. (16) Goraj, S.; Hubbard, W.; Reina, A.; Kong, J.; Branton, D.; Golovchenko, J. A. Graphene as a Subnanometer Trans-Electrode Membrane. Nature 2010, 467, 190−193. (17) Makra, I.; Gyurscányi, R. E. Electrochemical Sensing with Nanopores: A Mini Review. Electrochem. Commun. 2014, 43, 55−59. (18) Gin, D. L.; Noble, R. D. Designing the Next Generation of Chemical Separation Membranes. Science 2011, 332, 674−676. (19) Faraudo, J. The Missing Link Between the Hydration Force and Interfacial Water: Evidence from Computer Simulations. Curr. Opin. Colloid Interface Sci. 2011, 16, 557−560. (20) Thomas, M.; Corry, B.; Hilder, T. A. What Have We Learnt About the Mechanisms of Rapid Water Transport, Ion Rejection and Selectivity in Nanopores from Molecular Simulation? Small 2014, 10, 1453−1465. (21) Park, J. H.; Sinnott, S. B.; Aluru, N. R. Ion Separation Using a YJunction Carbon Nanotube. Nanotechnology 2006, 17, 895−900. (22) Zhao, K.; Wu, H. Size Effects of Pore Density and Solute Size on Water Osmosis through Nanoporous Membrane. J. Phys. Chem. B 2012, 116, 13459−13466. (23) Hao, L.; Su, J.; Guo, H. Water Permeation Through a Charged Channel. J. Phys. Chem. B 2013, 117, 7685−7694. (24) Yang, D.; Liu, Q.; Li, H.; Gao, C. Molecular Simulation of Carbon Nanotube Membrane for Li+ and Mg2+ Separation. J. Membr. Sci. 2013, 444, 327−331. (25) Wang, L.; Dumont, R. S.; Dickson, J. M. Nonequilibrium molecular dynamics simulation of pressure-driven water transport through modified CNT membranes. J. Chem. Phys. 2013, 138, No. 124701. (26) Ohba, T.; Kanoh, H. Energetic Contribution to Hydration Shells in One-Dimensional Aqueous Electrolyte Solution by Anomalous Hydrogen Bonds. Phys. Chem. Chem. Phys. 2013, 15, 5658−5663. (27) Tang, D.; Kim, D. Temperature Effect on Ion Selectivity of Potassium and Sodium Ions in Solution. Chem. Phys. 2014, 428, 14− 18. (28) Suk, M. E.; Aluru, N. R. Water Transport Through Ultrathin Graphene. J. Phys. Chem. Lett. 2010, 1, 1590−1594. (29) Zhao, S.; Xue, J.; Kang, W. Ion Selection of Charge-Modified Large Nanopores in a Graphene Sheet. J. Chem. Phys. 2013, 139, 114702.

(30) Picallo, C. B.; Gravelle, S.; Joly, L.; Charlaix, E.; Bocquet, L. Nanofluidic Osmotic Diodies: Theory & Molecular Dynamics Simulations. Phys. Rev. Lett. 2013, 111, 244501. (31) He, Z.; Zhou, J.; Lu, X.; Corry, B. Bioinspired Graphene Nanopores with Voltage-Tunable Ion Selectivity for Na+ and K+. ACS Nano 2013, 7, 10148−10157. (32) Chen, H.; Ruckenstein, E. Nanomembrane Containing a Nanopore in an Electrolyte Solution: A Molecular Dynamics Approach. J. Phys. Chem. Lett. 2014, 5, 2979−2982. (33) Lorenz, C. D.; Travesset, A. Charge Inversion of Divalent Ionic Solutions in Silica Channels. Phys. Rev. E 2007, 75, 061202. (34) Lorenz, C. D.; Crozier, P. S.; Anderson, J. A.; Travesset, A. Molecular Dynamics of Ionic Transport and Electrokinetic Effects in Realistic Silica Channels. J. Phys. Chem. C 2008, 112, 10222−10232. (35) Cruz-Chu, E. R.; Aksimentiev, A.; Schulten, K. Ionic Current Rectification through Silica Nanopores. J. Phys. Chem. C 2009, 113, 1850−1862. (36) Shirono, K.; Tatsumi, N.; Daiguji, H. Molecular Simulation of Ion Transport in Silica Nanopores. J. Phys. Chem. B 2009, 113, 1041− 1047. (37) Argyris, D.; Cole, D. R.; Striolo, A. Ion-Specific Effects under Confinement: The Role of Interfacial Water. ACS Nano 2010, 4, 2035−2042. (38) Carr, R.; Comer, J.; Ginsberg, M. D.; Aksimentiev, A. Modeling Presure-Driven Transport of Proteins through a Nanochannel. IEEE Trans. Nanotechnol. 2011, 10, 75−82. (39) Videla, P. E.; Sala, J.; Marti, J.; Guárdia, E.; Laria, D. Aqueous Electrolytes Confined within Functionalized Silica Nanopores. J. Chem. Phys. 2011, 135, 104503. (40) Bourg, I. C.; Steefel, C. I. Molecular Dynamics Simulations of Water Structure and Diffusion in Silica Nanopores. J. Phys. Chem. C 2012, 116, 11556−11564. (41) Zhu, H.; Ghoufi, A.; Szymczyk, A.; Balannec, B.; Morineau, D. Computation of the Hindrance Factor for the Diffusion for Nanoconfined Ions: Molecular Dynamics Simulations versus Continuum-Based Models. Mol. Phys. 2012, 40, 1107−1114. (42) Bonnaud, P. A.; Coasne, B.; Pellenq, R. J.-M. Solvated Calcium Ions in Charged Silica Nanopores. J. Chem. Phys. 2012, 137, No. 064706. (43) Ho, T. A.; Argyris, D.; Cole, D. R.; Striolo, A. Aqueous NaCl and CsCl Solutions Confined in Crystalline Slit-Shaped Silica Nanopores of Varying Degree of Protonation. Langmuir 2012, 28, 1256−1266. (44) Haria, N. R.; Lorenz, C. D. Ion Exclusion and Electrokinetic Effects Resulting from Electro-Osmotic Flow of Salt Solutions in Charged Silica Nanopores. Phys. Chem. Chem. Phys. 2012, 14, 5935− 5944. (45) Rizzi, F.; Jones, R. E.; Debusschere, B. J.; Knio, O. M. Uncertainity Quantification in MD Simulations of Concentration Driven Ionic Flow through a Silica Nanopore. I. Sensitivity to Physical Parameters of the Pore. J. Chem. Phys. 2013, 138, 194104. (46) Renou, R.; Ghoufi, A.; Szymczyk, A.; Zhu, H.; Neyt, J.-C.; Malfreyt, P. Nanoconfined Electrolyte Solutions in Porous Hydrophilic Silica Membranes. J. Phys. Chem. C 2013, 117, 11017−11027. (47) Phan, A.; Ho, T. A.; Cole, D. R.; Striolo, A. Molecular Structure and Dynamics in Thin Water Films at Metal Oxide Surfaces: Magnesium, Aluminum, and Silicon Oxide Surfaces. J. Phys. Chem. C 2012, 116, 15962−15973. (48) Turgman-Cohen, S.; Araque, J. C.; Hoek, E. M. V.; Escobedo, F. A. Molecular Dynamics of Equilibrium and Pressure-Driven Transport Properties of Water through LTA-Type Zeolites. Langmuir 2013, 29, 12389−12399. (49) Liu, L.; Chen, X. Fast Ion Transport and Phase Separation in a Mechanically Driven Flow of Electrolytes through Tortuous SubNanometer Nanochannels. ChemPhysChem 2013, 14, 2413−2418. (50) Liu, J.; Wang, M.; Chen, S.; Robbins, M. O. Molecular Simulations of Electroosmotic Flows in Rough Nanochannels. J. Comput. Phys. 2010, 229, 7834−7847. 12310


Article

The Journal of Physical Chemistry C (51) Jin, X.; Aluru, N. R. Gated Transport in Nanofluidic Devices. Microfluid. Nanofluid. 2011, 11, 297−306. (52) De Biase, P. M.; Markosyan, S.; Noskov, S. Microsecond Simulations of DNA and Ion Transport in Nanopores with Novel IonIon and Ion-Nucleotides Effective Potentials. J. Comput. Chem. 2014, 35, 711−721. (53) Zhang, H.; Hassanali, A. A.; Shin, Y.-K.; Knight, C.; Singer, S. J. The Water-Amorphous Silica Interface: Analysis of the Stern Layer and Surface Conduction. J. Chem. Phys. 2011, 134, 024705. (54) Qiao, R.; Aluru, N. R. Atomistic Simulation of KCl Transport in Charged Silicon Nanochannels: Interfacial effects. Colloids Surf., A 2005, 267, 103−109. (55) van Beest, B. W. H.; Kramer, G. J.; van Santen, R. A. Force Fields for Silicas and Aluminophosphates Based on Ab Initio Calculations. Phys. Rev. Lett. 1990, 64, 1955−1958. (56) Iler, R. K. The Chemistry of Silica; Wiley: New York, 1979. (57) Zhuravlev, L. T. The Surface Chemistry of Amorphous Silica. Zhuravlev Model. Colloids Surf., A 2000, 173, 1−38. (58) Tadros, T. F.; Lyklema, J. Adsorption of Potenial-Determining Ions at Silica-Aqueous Electrolyte Interface and Role of Some Cations. J. Electroanal. Chem. 1968, 17, 267. (59) van der Heyden, F. H. J.; Stein, D.; Besteman, K.; Lemay, S. G.; Dekker, C. Charge Inversion at High Ionic Strength Studied by Streaming Currents. Phys. Rev. Lett. 2006, 96, 224502. (60) Cruz-Chu, E. R.; Aksimentiev, A.; Schulten, K. Water-Silica Force Field for Simulating Nanodevices. J. Phys. Chem. B 2006, 110, 21497−21508. (61) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (62) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (63) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (64) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (65) Mahoney, M. W.; Jorgensen, W. L. A Five-Site Model for Liquid Water and the Reproduction of the Density Anomaly by Rigid, Nonpolarizable Potential Functions. J. Chem. Phys. 2000, 112, 8910. (66) Vega, C.; Abascal, J. L. F. Simulating Water with Rigid NonPolarizable Models: A General Perspective. Phys. Chem. Chem. Phys. 2011, 13, 19663−19688. (67) van der Spoel, D.; van Maaren, P. J.; Berendsen, H. J. C. A Systematic Study of Water Models for Molecular Simulation: Derivation of Water Models Optimized for Use with a Reaction Field. J. Chem. Phys. 1998, 108, 10220−10230. (68) Grossfield, A.; Ren, P.; Ponder, J. W. Ion Solvation Thermodynamics from Simulation with a Polarizable Force Field. J. Am. Chem. Soc. 2003, 125, 15671−15682. (69) Busch, S.; Pardo, L. C.; O’Dell, W. B.; Bruce, C. D.; Lorenz, C. D.; McLain, S. E. On the structure of water and chloride ion interactions with a peptide backbone in solution. Phys. Chem. Chem. Phys. 2013, 15, 21023−21033. (70) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesion Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (71) Hockney, R. W.; Eastwood, J. W. Computer Simulation Using Particles; McGraw-Hill: New York, 1981. (72) Greberg, H.; Kjellander, R. Charge Inversion in Electric Double Layers and Effects of Different Sizes for Counterions and Coions. J. Chem. Phys. 1998, 108, 2940−2953. (73) Besteman, K.; Zevenbergen, M. A. G.; Lemay, S. G. Charge Inversion by Multivalent Ions: Dependence on Dielectric Constant and Surface-Charge Density. Phys. Rev. E 2005, 72, 061501.

(74) Faraudo, J.; Travesset, A. The Many Origins of Charge Inversion in Electrolyte Solutions: Effects of Discrete Interfacial Charges. J. Phys. Chem. C 2007, 111, 987−994. (75) Roux, B.; Allen, T.; Bernéche, S.; Im, W. Theoretical and Computational Models of Biological Ion Channels. Q. Rev. Biophys. 2004, 37, 15−103. (76) Maffeo, C.; Bhattacharya, S.; Yoo, J.; Wells, D.; Aksimentiev, A. Modelling and Simulation of Ion Channels. Chem. Rev. 2012, 112, 6250−6284.

12311


Atomistic Description of Pressure-Driven Flow of ... - ACS Publications

Recommend Documents