Charge Asymmetry Effect in Ion Transport through Angstrom-Scale

Jan 3, 2019 - We believe that these results would shed light on future design and ... which we believe should be responsible for the charge asymmetry ...
0 downloads 0 Views 2MB Size
Article Cite This: J. Phys. Chem. C XXXX, XXX, XXX−XXX

pubs.acs.org/JPCC

Charge Asymmetry Effect in Ion Transport through Angstrom-Scale Channels YanZi Yu,† JingCun Fan,† Ali Esfandiar,‡ YinBo Zhu,† HengAn Wu,† and FengChao Wang*,† †

CAS Key Laboratory of Mechanical Behavior and Design of Materials, Department of Modern Mechanics, CAS Center for Excellence in Complex System Mechanics, University of Science and Technology of China, Hefei, Anhui 230027, China ‡ Department of Physics, Sharif University of Technology, Tehran 11155-9161, Iran

J. Phys. Chem. C Downloaded from pubs.acs.org by UNIV DE BARCELONA on 01/09/19. For personal use only.

S Supporting Information *

ABSTRACT: Structural and dynamic properties of ions confined in nanoslits are crucial to understand the fundamental mechanism underlying a wide range of chemical and biological phenomena. K+ and Cl− show similar ion mobilities in a bulk aqueous solution, whereas they exhibit a remarkable difference when transporting through an angstrom-scale channel. Our molecular dynamics simulations uncover that such discrepancy originates from the subtle differences in their hydration structures and preferable locations across the channel. Opposite charge causes different water dipolar orientations around ions, mediating the distance and tribological interactions between hydrated ions and channel’s walls. Hydrated Cl− ions experience a remarkable larger friction force inside the channel and consequently a smaller mobility compared with K+ ions. This charge asymmetry mechanism becomes dominant when the channel size approaches the molecular scale, which takes responsibility for the delicate distinction between dynamic transport behaviors of K+ and Cl− under strong confinement. We believe that these results would shed light on future design and optimization of new membranes for filtration and desalination applications.



pected.25−27 More recently, angstrom-scale channels were fabricated using 2-dimensional materials assembled into devices with controllable size, atomically smooth inner walls, and little surface charge.28 It is interesting and somewhat surprising that experiments implemented on those channels show that ion mobilities of K+ and Cl− exhibit a threefold difference,28 although K+ and Cl− show very similar hydrated diameters14 and ion mobilities in bulk aqueous solution.29 Bulk water is often modeled as a dielectric continuum, which is macroscopic isotropy.30 On the contrary, confined water in a channel of height 6.8 Å manifests a distinct bilayer structure.31 In the presence of ions, consequences of competing interactions among ions, water molecules, and walls would result in greater complexity of the local configuration, a compromise between water layers and hydration shells.32,33 Moreover, water reorientation and hydrogen-bond exchange around opposite charges must be taken into account.34,35 Apart from the ion mobility, analogous charge asymmetry effect on hydration energy has been reported.36,37 Despite significant progress from various experiments in this field, theoretical modeling on a molecular level is still needed to complete our understanding on this topic. In the present study, we have aimed to essentially answer the following

INTRODUCTION Ion sieving and filtration have gained massive attention worldwide from pure research to industrial and medical applications, such as water desalination,1,2 energy harvesting and storage,3−5 transmembrane ion transport,6 and hemodialysis.7 One of the basic principles involved in those processes is size exclusion, in which the ion transport is regulated by the size competition between channels and hydrated ions.8−10 When ions are dissolved in water, they interact with surrounding water molecules and form structural hydration shells.11 The effective hydrated diameter,12 DH, is widely adopted to characterize the size of hydrated ions. Although its value may vary slightly depending on the specific method, DH is generally in the range between several angstroms to about 1 nm for common ions.13,14 When hydrated ions try to pass a channel with size smaller than DH, the size exclusion effect leads to the complete rejection or, in a broad sense, causes larger energy barriers and changes the selectivity of certain ions.15−17 The full picture of hydrated ions passing through nanochannels is rather complex and has not been fully elucidated, especially when dimensions of channel approach the size of small ions and individual water molecules.18−22 Under strong confinement, the hydration shells undergo rearrangement, including distortion or partial dehydration, to fit themselves inside a channel.23,24 Influence of structural adjustment of hydrated ions on their dynamic behaviors can be ex© XXXX American Chemical Society

Received: October 5, 2018 Revised: December 17, 2018

A

DOI: 10.1021/acs.jpcc.8b09742 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C question: how hydrated K+ and Cl− ions modulate their structures and then diffuse through Å-scale channels? By using molecular dynamics (MD) simulations, we demonstrate that across a channel that accommodates only two layers of water, K+ and Cl− ions show different preferable locations, that is, inter-layer and intra-layer relative to the bilayer water, respectively. Our results further reveal the distinguishing frictional behaviors between hydrated ions and confining walls, which we believe should be responsible for the charge asymmetry phenomenon of ion mobilities under angstromscale confinement.



COMPUTATIONAL DETAILS Classical MD Simulations. In this work, we focus on the ion transport through an angstrom-scale slit with the effective channel height of h = 6.8 Å, which is corresponding to the fabricated nanocapillaries with two layers of graphene as the spacer.28 To calculate the current−voltage (I−V) relations, we used a simulation system analogous to that employed in the recent experiments,28 as shown in Figure 1a. The simulation box has a size of 68.0 Å × 34.4 Å × 80.0 Å. A graphite slab composed of 20 graphene layers separates the simulation box into 3 parts. Two layers in the middle were deleted; thus, it leaves us the angstrom-scale capillary (h = 6.8 Å) with a length of 46.2 Å for the ion transport. As a comparison reference, we also studied the ion transport through an h = 61.2 Å with 18 layers of graphene sheets removed. Although flexible graphene sheets seem closer to the practical situation, the carbon atoms in our simulations were fixed to their initial positions and were modeled as simple particles.38 Unlike the much thinker top layer (∼100 nm) used in experiments,28 which leads to a rapidly increasing bending rigidity with its thickness,39 here rigid channel walls would prevent the possible sagging or even closure of the channel. Apart from the region containing the graphene sheets, the simulation box was filled with the aqueous solution of KCl with a concentration of 1 M. The simple point charge extended model was adopted for water molecules.40 Forcefield parameters for K+ and Cl− ions come from ref 41. LennardJones potentials are used to describe interactions between different particles, like water−graphene.42 Parameters not explicitly specified are obtained using the Lorentz−Berthelot combination rules. Classical MD simulations were performed using LAMMPS.43 Periodic boundary conditions were applied to all three directions. The short-range interactions were truncated at a cutoff of 12 Å, and the long-range Coulomb interactions were computed by utilizing the particle−particle particle−mesh algorithm.44 After the initial energy minimization, the system was first equilibrated in the isothermal− isobaric (NPT) ensembles for 10 ns with a timestep of 1 fs at 300 K and 1 atm, whereas the barostat was applied only along the Z direction parallel to the channel. Then, simulations were performed in the canonical (NVT) ensembles using a Nosé− Hoover thermostat. The external electric field was applied in the Z direction with various values ranging from 0.001 to 0.12 V/Å, corresponding to the voltage drop ranging from 0.08 to 9.6 V. We have examined that within this range of E, the movement of ions inside the capillaries is still at a constant speed. Electric Current Calculation. In this work, the ionic currents through capillaries were calculated using the following equation,45 which actually originates from the definition of the

Figure 1. MD simulations on the ion transport through channel with a height of 6.8 Å. (a) Perspective view of a snapshot of simulation systems. Graphene layers are shown in gray. Water molecules are represented by red (oxygen) and white (hydrogen) spheres. K+ and Cl− ions are represented by purple and green spheres, respectively. (b) Calculated current−voltage (I−V) relations. (c) Individual contribution of Cl− to the total current. Red circles, blue diamonds, and green squares are simulation results on h = 6.8 Å, h = 61.2 Å, and bulk solution, respectively. Error bars indicating the standard deviation are displayed only if they are larger than the symbols. Shadow area in (c) shows the results from a reference simulation on the bulk solution.54

electric current, i.e., the electric charge Q transferred through a surface over time t. n

I (t ) =

1 ∑ q [zi(t + Δt ) − zi(t )] Δt ·LZ i = 1 i

(1)

where Lz is the length of the simulation box along the Z direction, and zi and qi are the z coordinate and the charge of the ion i, respectively. For each individual E, the z coordinates of all ions were recorded at regular time intervals, Δt = 1 ps, and the data from a simulation run of 5.0 ns were processed B

DOI: 10.1021/acs.jpcc.8b09742 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C and analyzed to calculate the electric current I. We also used eq 1 to calculate the individual contributions from cations or anions to the total current. Then, the summation on subscripts i applies only to K+ or Cl− ions accordingly. Ab Initio MD Simulations. We performed ab initio MD (AIMD) simulations to verify our classical MD results on the location distribution of K+ and Cl− ions and water inside the h = 6.8 Å channel along the transverse direction perpendicular to graphene walls. Because AIMD simulations are time consuming, here we constructed a much smaller system than the one used in classical MD simulations. The same ion concentration was maintained. As shown in Figure S1, it consists of an ion (K+ or Cl−), 30 water molecules, and two graphene sheets (120 carbon atoms), which were placed in 12.781 Å × 12.297 Å × 40.0 Å simulation box. AIMD simulations were carried out using CP2K,46 which is based on the Gaussian and plane wave method for ab initio Born− Oppenheimer MD within the Kohn−Sham framework of density functional theory. Apart from the force evaluation, other procedures were substantially similar to those of classical MD simulations. The Perdew−Burke−Ernzerhof and the van der Waals potential with density functional dispersion correction of type 3 (DFT-D3) were utilized as the exchange and correlation function.47,48 The molecularly optimized basis functions (MOLOPT) Gaussian basis (with an energy cutoff of 280 Ry) was selected and the Goedecker, Teter, and Hutter pseudopotentials were used to treat the core electrons.49−51 Periodic boundary conditions were applied in all three directions. The geometry of graphene sheets was fixed to their initial positions. In each AIMD run, the temperature was controlled at 300 K, and the time step was 1.0 fs. First, the systems were equilibrated in the canonical ensemble (NVT) using the Nosé−Hoover thermostat for 5.0 ps. Then the simulations were continued with energy-conserving dynamics in the microcanonical ensemble (NVE) for 15 ps to produce data for subsequent analysis. Friction Force Calculation. To investigate the tribological interactions between hydrated ions with channel’s walls, we carried out classical MD simulations on an ion−water cluster, which was dragged by a spring. The other end of this spring moves at a constant speed. Thus, the dragging force applied on the ion−water cluster is variable and its magnitude is the product of the elastic coefficient and the elongation of the spring. This is a typical procedure widely used to investigate the atomic-scale friction behavior.52,53 The dragging force increases from zero until it exceeds the maximum static friction that the channel walls can provide. After a sudden slip, the length of the spring shrinks, and the spring force drops below the critical value. Thus, the ion−water cluster sticks at the adjacent equilibrium position. Before the ion−water cluster moves, the spring force is balanced by the friction force. The reaction force on the channel walls were computed and converted to the magnitude of the friction force applied on the ion−water cluster.

indicates that a large amount of ion passage is blocked owing to the size effect. Nevertheless, the I−V curves in Figure 1b only reflect the sum of K+ and Cl− ion transport behavior. It is instructive to separate and compare these two contributions. MD simulations have this advantage that the individual contributions from anions and cations can be calculated directly using eq 1. The percentage of Cl− contribution to the currents were calculated as a ratio of ICl−/Itotal and provided under various applied voltages in Figure 1c. Although the K+ and Cl− ions move in opposite directions under the action of the electric field, they both provide positive contributions to the total current, which can also be derived from eq 1. The ion mobility is defined as the ratio of the drift velocity vd to the magnitude of the applied electric field: μ = vd/E. Thus, we can connect the ionic current to the ion mobilities by I = μcνEFS, where μ is the ion mobility, c is the concentration, ν is the valence, E is the applied electric field, F is the Faraday constant, and S is the cross-sectional area of the channel. For monovalent ions K+ and Cl− that we investigated here, the concentration is assumed to be identical. Thus, theoretically, the individual contribution is proportional to the ion mobility. Experimental measurements on infinite dilution demonstrate that ion mobilities for K+ and Cl− are 7.57 × 10−8 and 7.86 × 10−8 m2·V−1·s−1, respectively.29 Consequently, the individual contribution of either K+ or Cl− should be around 50%. This conjecture is supported by our simulation results obtained from the bulk solution. As shown in Figure 1c, our results processed using eq 1 manifest that the contribution of Cl− to the total current remains in a range of 45−49%. Moreover, there is no obvious dependency of the contribution of Cl− on the applied voltage. On the other hand, the ion mobility can be calculated using the Einstein relation,54,55 μ = Dq/kBT, in which q is the ionic charge, kB is the Boltzmann constant, and D is the ionic diffusion coefficient. D is proportional to the slope of mean-squared displacement versus time, which can also be obtained from MD simulations. In this way, the ion mobilities of K+ and Cl− in bulk solution from our simulations give the estimation of this value to be 46.2 ± 0.15% (see Supporting Information), which also agrees well with experiments29 and MD simulations from others.54 In contrast, under strong confinement, we observed anomalous behavior. At low voltage, it is clear that the contribution of Cl− to the total current is only within 30−40%, a noticeable deviation from ∼50%, as expected from the ion mobilities in bulk solution. Only when the applied voltage is even larger (>2 V in our cases), the contribution of Cl− approaches that calculated from bulk solution. This can be explained by the fact that a large driving force would weaken the discrepancy caused by ion mobility. Here we should note that one cannot use I = μcνEFS directly to calculate the ionic currents in a similar system, because ion mobilities inside nanocapillaries are not known beforehand. The confinement effect would significantly change the magnitude of the ion mobility compared with its bulk value. Although the walls in our simulations are strictly neutral, the results shown in Figure 1c qualitatively agree with the experimental observations on the charge asymmetry effect of ion mobilites reported in ref 28. Regardless of their opposite charge, K+ and Cl− show very similar hydrated diameter14 and nearly identical ion mobilities in bulk aqueous solutions.29,54 However, as presented in Figure 1c, K+ and Cl− exhibit notable asymmetry transport behaviors under strong confinement. Why is the transport of negative



RESULTS AND DISCUSSION The current−voltage (I−V) relations for the aqueous solution of KCl with a concentration of 1 mol/L are plotted in Figure 1b. To give a more intuitive comparison, we also provided I−V curves of bulk solution of the same concentration, as well as simulation results from a much larger channel, h = 61.2 Å. From Figure 1b, we can find that the currents through angstrom-scale slits were significantly suppressed, which C

DOI: 10.1021/acs.jpcc.8b09742 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Cl− ions through nanocapillaries suppressed more severely? This question motivated us to explore the underlying mechanism. To do this, we first examine the atomic structure of hydrated ions. It has been well established that dissolved ions in aqueous solution are hydrated with surrounding water molecules,11 as illustrated in Figure 2a. MD simulation results

between the hydration structure of ions and channel’s walls, which would be discussed later. Intermolecular forces are dominant on this atomic scale.30 MD simulations confirm that there is a distinct bilayer structure for confined water in a channel of height 6.8 Å.31 This bilayer structure may play a significant role in the structure of hydration shells. With these considerations in mind, we have to take into account the interactions among ions, water molecules, and walls. To test this postulation, we carried out additional simulations on location distribution of ions inside channels along the transverse direction (perpendicular to graphene walls). For the h = 6.8 Å slits, the center of carbon atoms is located at ±5.1 Å. The distribution profiles of O and H clearly confirm a bilayer structure of the confined water. Additionally, we found that the effective thickness occupied by the confined water is about 6.0 Å or even smaller. Here we should notice that h = 6.8 Å is just an empirical value, which is obtained by subtracting the effective thickness of one graphene layer 3.4 Å from the distance between the atomic centers of two innermost graphene sheets, 10.2 Å. Our results in Figure 3a demonstrate that K+ ions have a concentrated distribution between these two water layers (inter-layer), whereas Cl− ions prefer to stay in either one (intra-layer), closer to the graphene walls. For the distribution of H with the presence of K+, the long tails toward graphene walls also verify

Figure 2. Ionic hydration structures of K+ and Cl−. (a) Schematics of hydrated K+ and Cl− ions in bulk solutions. The first and second hydration shells are represented with pink and green shadow, respectively. (b) The RDF of hydrogen (solid lines) and oxygen (dash lines) atoms in water molecules around K+ (red lines) and Cl− (black lines) in the bulk aqueous solution.

of the radial distribution function (RDF) confirm such hypothesis. We can find at least two hydration shells along the radial direction (Figure 2b). Opposite charge results in different polarization of water molecules in hydration shells. The negative (oxygen) side of a water molecule is attracted by K+ with OH bonds pointing outward, whereas the positive (hydrogen) ends of water molecules are attracted to Cl−, as shown in Figure 2a. Such structural difference is also reflected in the RDF curves. The first peak of K+−O RDF can be found at 2.75 Å, smaller than that of K+−H RDF curve, 3.35 Å. Similar results can also be found in Cl−−O and Cl−−H RDF curves, 3.15 Å versus 2.25 Å, respectively. If we take the first minimum of K+−H and Cl−−O RDF curves to define their effective size, we may find that the diameter of the first/second hydration shells for K+ and Cl− ions are 8.3/13.5 and 7.7/12.9 Å, respectively. Compared with the effective channel size 6.8 Å, we can deduce that the hydration shells cannot maintain spherical structures as they do in bulk solutions. This will lead to two possible consequences. (i) There is not enough space to accommodate the entire second hydration shell, indicating a partial dehydration. Dehydration has been related to the ion selectivity of nanopores,15 which also contributes to the mechanism of charge asymmetry effect on ion mobilities. (ii) In addition, even for the first hydration shell, it is very likely to exhibit extrusion and deformation, causing intricate interplay

Figure 3. Location distribution of K+ and Cl− ions and water inside the 6.8 Å channel along the transverse direction perpendicular to graphene walls. (a) Classical MD simulation results. (b) AIMD simulations. Histograms in the upper panel show the distribution of K+ (red) and Cl− ions (black). Distributions of hydrogen (orange) and oxygen (green) atoms in water molecules accompanied by K+ and Cl− ions are provided in the middle and lower panels, respectively. D

DOI: 10.1021/acs.jpcc.8b09742 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

We believe that the subtle difference in hydration structures between K+ and Cl− would play a crucial role in the charge asymmetric effect on ion mobilities. To this end, we investigated their tribological behaviors inside channels. The hydration structures (ion with surrounding water molecules) were extracted from the aqueous environment. It should be noted that the extraction of hydration clusters may weaken the foregoing structural competition with the bilayer water. Nevertheless, this extraction allows us to calculate straightforwardly the friction force with graphene walls. Otherwise, the dynamic exchange of water molecules from hydration shells with others makes the calculation of friction forces impossible. During the simulations, the ion (K+ or Cl−) was connected by a spring whose another end was pulled with a constant velocity (0.005 Å/ps), illustrated in Figure 5a. Unless

that herein OH bonds are directed outward. It is the typical feature of the orientation polarization caused by cations. We also provided simulation results from AIMD in Figure 3b. AIMD simulations are anticipated to yield more accurate results because the interactions are evaluated based on the density functional theory.56 Generally, these two types of simulations are in good agreement and come to an identical judgment; i.e., K+ ions favor the inter-layer location, whereas Cl− ions incline the intra-layer accommodation. The only obvious difference between Figure 3a,b is that there is one peak of location distribution of Cl− ions, rather than two. This is not surprising because of the relatively short simulation duration of AIMD, 15 ps versus 1 ns of the classical MD. Cl− ions need longer time to shuttle between two water layers. The distinct distribution of cations and anions in narrow slits are consistent with previous MD simulations using different ions, water models, and different parameters for ion−wall interactions,17,57,58 which indicate that this finding is robust. In fact, we think it is the balance of structural competition between the hydration shell and the bilayer water that is crucial to interpret the charge asymmetry effect on ion mobilities. On the basis of data at hand, we sketched the energetically favorable hydration structures of K+ and Cl− ions in a channel with an effective height of 6.8 Å, as shown in Figure 4a. In

Figure 4. (a) Schematics of hydration structures of K+ and Cl− ions inside the h = 6.8 Å channel, based on the simulation results of the location distribution of ions and water molecules. (b) Schematics of energetically unstable hydration structures of K+ and Cl− ions under strong confinement. The first and second hydration shells are represented with pink and green shadow, respectively.

Figure 5. Dependence of tribological behaviors of the ion−water cluster on the transverse position and the charge. (a) Snapshot of simulation model. (b) Friction forces on various hydration structures. Upper panel: K+ and Cl− ions located at their most preferable positions. Middle panel: position effect on Cl− ions. Lower panel: opposite charge effect.

order to exclude the antithetical cases, we would discuss the intra-layer hydrated K+ and inter-layer hydrated Cl− ions, respectively. In order to exclude other possibilities, let us first assume that K+ ions stay in one of the water layers. Then the polarization requires that some water molecules settle between two water layers, to satisfy the K+−O spacing of 2.75 Å, given by the first peak of K+−O RDF curve shown in Figure 2b. Such inter-layer water is definitely unstable, because it is in conflict with the bilayer structure (see Figure 4b). On the other hand, oxygen distribution in Figure 3 tells us that the spacing between the two water layers is predominantly less than 4.0 Å. Analogously, Cl− ions staying at the middle of the channel lead to the distance between two opposite water molecules in the first hydration shell larger than 6.5 Å, which are also energetically unstable (see Figure 4b).

otherwise specified, the vertical distance between the ion and graphene walls was fixed in accordance with the peak value of the distribution probability given in Figure 3a. Even for the atomically smooth surfaces of graphene, the corrugation of its potential energy profile would lead to the atomic-scale friction. The friction force on the entire ion−water cluster was calculated and exhibits a stick-slip behavior, as shown in Figure 5b. Because the applied spring force increases from zero, the ion−water cluster would not move until the spring force is larger enough to overcome the energy barrier. After this, the cluster “slips” to the next equilibrium position. Thus, the friction force increases linearly in the short distance. Then the friction force drops because of the shrink of the spring. E

DOI: 10.1021/acs.jpcc.8b09742 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

inclined to stay inside either of the water layers. The discrepancy is because of a balance of competition between the bilayer water and the different polarization of hydration structures caused by the opposite charge. It further leads to the contrasting friction forces of the distorted hydration structures against walls. For this reason, the ion mobility of Cl− under strong confinement is smaller than that of K+. Our simulation results, which are qualitatively consistent with those of the latest experiments, emphasize the decisive influence of subtle difference residing in structures on dynamic properties of hydrated ions, because the channel size is comparable with their hydration diameter. We hope these findings would enhance our understanding of ion transport under strong confinement and provide theoretical guide for nanofiltration and desalination applications.

This process repeats; thus, the friction force as a function of the displacement of the spring’s stretching end exhibits the undulation over larger length scale. We have the following findings from the results shown in Figure 5b. First, our results show that the friction force on the hydration cluster of Cl− is significantly larger than that of K+. This means that hydrated Cl− ions are more difficulty to move inside channels, which is consistent with experiments28 and our simulation results shown in Figure 1c. Second, the preferable location of hydrated ions has a significant impact on the tribological dynamics. Because inter-layer Cl− ions with low probability also exist in Figure 3a, the friction forces on the hydration clusters of inter-layer and intra-layer Cl− ions were calculated. We found that the friction force on the hydration cluster of inter-layer Cl− decreases remarkably (the middle panel of Figure 5b). Third, we further confirm that only the polarization caused by the opposite charge is not enough to explain the difference in the friction forces and ion mobilities. We compared the friction forces on the hydration clusters of inter-layer K+ and K− ions. No significant differences were found, as shown in the lower panel of Figure 5b. Although K− seems to be unphysical, it is very convenient to do this in MD simulations, which allow us to directly compare the effect of asymmetric charge while keeping other parameters unchanged. For the results given in Figures 3a and 5b, we also tested other choices of forcefield parameters59 and obtained similar results as reported here (see the Supporting Information). The aforementioned results and discussion establish that the charge asymmetry effect on ion transport through angstromscale channels can be interpreted by the different tribological interactions with channel walls, exerted on the inter-layer K+ and intra-layer Cl− hydrated ions. These findings are qualitatively consistent with experimental results,28 although with some differences. Experiments show that there are tiny but finite negative surface charge in graphene capillary.28 However, we didn’t discuss the surface charge effect here and the channel walls were modeled to be neutral. It should be mentioned that extra attention is required when dealing with nanochannels with inherent surface change like boron nitride nanotubes,4 graphene oxide laminate,20 and calcium silicate hydrate nanopores.60,61 In another aspect, there have been massive investigations on the cation−π interaction in a system involving graphene sheets or carbon nanotubes.10,62 We have checked that our AIMD simulations are capable of capturing the characteristics of ion−π interactions in our systems. Nevertheless, cation−π interactions would lead to a closer distance and larger friction force with graphene for the K+ ion, which is inconsistent with what is reported in Figure 5. This is not surprising, because we focus on the ion transport through an angstrom-scale capillary with two graphene surfaces. In this scenario, cation−π interactions from the opposite walls may offset each other to some extent. Thus, cation−π interactions are not appropriate to interpret our main finding in the present work.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.8b09742. Simulation modeling used for AIMD calculations; ion mobilities calculated using the Einstein relation; more test with other choice of the force field parameters (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +86 551 6360 1238. ORCID

YinBo Zhu: 0000-0001-9204-9300 FengChao Wang: 0000-0002-5954-3881 Author Contributions

F.C.W. designed and directed the project. Y.Z.Y. performed simulations and their analysis. J.C.F. provided the postprocessing code. F.C.W. and Y.Z.Y. wrote the article. All authors contributed to discussions. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB22040402), the National Natural Science Foundation of China (11772319, 11572307, 11525211, 11472263, 11802302), the Fundamental Research Funds for the Central Universities (WK2090050040, WK2090050043), and the Young Elite Scientists Sponsorship Program by CAST (2016QNRC001). Y.B.Z. acknowledges the National Postdoctoral Program for Innovative Talents (BX201700225), the Anhui Provincial Natural Science Foundation (1808085QA07), and the China Postdoctoral Science Foundation (2017M620264). The numerical calculations have been performed on the supercomputing system in the Supercomputing Center of University of Science and Technology of China.



CONCLUSIONS On the basis of the molecular details already accumulated, we uncover the atomic origin of charge asymmetry effect on ion transport through angstrom-scale channels. To adjust themselves to the inherent bilayer water structure in a h = 6.8 Å channel, hydrated ions choose particular preferable distribution across the channel. K+ ions are more likely to locate at the center of the channel, whereas Cl− ions are



REFERENCES

(1) Cohen-Tanugi, D.; Grossman, J. C. Water Desalination across Nanoporous Graphene. Nano Lett. 2012, 12, 3602−3608. (2) Surwade, S. P.; Smirnov, S. N.; Vlassiouk, I. V.; Unocic, R. R.; Veith, G. M.; Dai, S.; Mahurin, S. M. Water Desalination Using

F

DOI: 10.1021/acs.jpcc.8b09742 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Nanoporous Single-Layer Graphene. Nat. Nanotechnol. 2015, 10, 459−464. (3) Maier, J. Nanoionics: Ion Transport and Electrochemical Storage in Confined Systems. Nat. Mater. 2005, 4, 805−815. (4) Siria, A.; Poncharal, P.; Biance, A.-L.; Fulcrand, R.; Blase, X.; Purcell, S. T.; Bocquet, L. Giant Osmotic Energy Conversion Measured in a Single Transmembrane Boron Nitride Nanotube. Nature 2013, 494, 455−458. (5) Ou, X.; Yu, Y.; Wu, R.; Tyagi, A.; Zhuang, M.; Ding, Y.; Abidi, I. H.; Wu, H.; Wang, F.; Luo, Z. Shuttle Suppression by Polymer-Sealed Graphene-Coated Polypropylene Separator. ACS Appl. Mater. Interfaces 2018, 10, 5534−5542. (6) Gouaux, E.; MacKinnon, R. Principles of Selective Ion Transport in Channels and Pumps. Science 2005, 310, 1461−1465. (7) Lehninger, A. L. Mitochondria and Calcium Ion Transport. Biochem. J. 1970, 119, 129−138. (8) Zwolak, M.; Lagerqvist, J.; Di Ventra, M. Quantized Ionic Conductance in Nanopores. Phys. Rev. Lett. 2009, 103, 128102. (9) Abraham, J.; Vasu, K. S.; Williams, C. D.; Gopinadhan, K.; Su, Y.; Cherian, C. T.; Dix, J.; Prestat, E.; Haigh, S. J.; Grigorieva, I. V.; Carbone, P.; Geim, A. K.; Nair, R. R. Tunable Sieving of Ions Using Graphene Oxide Membranes. Nat. Nanotechnol. 2017, 12, 546−550. (10) Chen, L.; Shi, G.; Shen, J.; Peng, B.; Zhang, B.; Wang, Y.; Bian, F.; Wang, J.; Li, D.; Qian, Z.; Xu, G.; Liu, G.; Zeng, J.; Zhang, L.; Yang, Y.; Zhou, G.; Wu, M.; Jin, W.; Li, J.; Fang, H. Ion Sieving in Graphene Oxide Membranes via Cationic Control of Interlayer Spacing. Nature 2017, 550, 380−383. (11) Frank, H. S.; Wen, W.-Y. Ion-Solvent Interaction. Structural Aspects of Ion-Solvent Interaction in Aqueous Solutions: A Suggested Picture of Water Structure. Discuss. Faraday Soc. 1957, 24, 133. (12) Kielland, J. Individual Activity Coefficients of Ions in Aqueous Solutions. J. Am. Chem. Soc. 1937, 59, 1675−1678. (13) Ohtaki, H.; Radnai, T. Structure and Dynamics of Hydrated Ions. Chem. Rev. 1993, 93, 1157−1204. (14) Tansel, B. Significance of Thermodynamic and Physical Characteristics on Permeation of Ions during Membrane Separation: Hydrated Radius, Hydration Free Energy and Viscous Effects. Sep. Purif. Technol. 2012, 86, 119−126. (15) Sahu, S.; Di Ventra, M.; Zwolak, M. Dehydration as a Universal Mechanism for Ion Selectivity in Graphene and Other Atomically Thin Pores. Nano Lett. 2017, 17, 4719−4724. (16) Corry, B. Mechanisms of Selective Ion Transport and Salt Rejection in Carbon Nanostructures. MRS Bull. 2017, 42, 306−310. (17) Zhou, K.; Xu, Z. Renormalization of Ionic Solvation Shells in Nanochannels. ACS Appl. Mater. Interfaces 2018, 10, 27801−27809. (18) Duan, C.; Majumdar, A. Anomalous Ion Transport in 2-nm Hydrophilic Nanochannels. Nat. Nanotechnol. 2010, 5, 848−852. (19) Humplik, T.; Lee, J.; O’Hern, S. C.; Fellman, B. A.; Baig, M. A.; Hassan, S. F.; Atieh, M. A.; Rahman, F.; Laoui, T.; Karnik, R.; Wang, E. N. Nanostructured Materials for Water Desalination. Nanotechnology 2011, 22, 292001. (20) Joshi, R. K.; Carbone, P.; Wang, F. C.; Kravets, V. G.; Su, Y.; Grigorieva, I. V.; Wu, H. A.; Geim, A. K.; Nair, R. R. Precise and Ultrafast Molecular Sieving through Graphene Oxide Membranes. Science 2014, 343, 752−754. (21) Feng, J.; Liu, K.; Graf, M.; Dumcenco, D.; Kis, A.; Di Ventra, M.; Radenovic, A. Observation of Ionic Coulomb Blockade in Nanopores. Nat. Mater. 2016, 15, 850−855. (22) Hong, S.; Constans, C.; Surmani Martins, M. V.; Seow, Y. C.; Guevara Carrió, J. A.; Garaj, S. Scalable Graphene-Based Membranes for Ionic Sieving with Ultrahigh Charge Selectivity. Nano Lett. 2017, 17, 728−732. (23) Shao, Q.; Zhou, J.; Lu, L.; Lu, X.; Zhu, Y.; Jiang, S. Anomalous Hydration Shell Order of Na+ and K+ inside Carbon Nanotubes. Nano Lett. 2009, 9, 989−994. (24) Choi, W.; Ulissi, Z. W.; Shimizu, S. F. E.; Bellisario, D. O.; Ellison, M. D.; Strano, M. S. Diameter-Dependent Ion Transport through the Interior of Isolated Single-Walled Carbon Nanotubes. Nat. Commun. 2013, 4, 2397.

(25) Lynden-Bell, R. M.; Rasaiah, J. C. Mobility and Solvation of Ions in Channels. J. Chem. Phys. 1996, 105, 9266−9280. (26) Zhou, Y.; Morais-Cabral, J. H.; Kaufman, A.; MacKinnon, R. Chemistry of Ion Coordination and Hydration Revealed by a K+ Channel−Fab Complex at 2.0 Å Resolution. Nature 2001, 414, 43− 48. (27) Sint, K.; Wang, B.; Král, P. Selective Ion Passage through Functionalized Graphene Nanopores. J. Am. Chem. Soc. 2008, 130, 16448−16449. (28) Esfandiar, A.; Radha, B.; Wang, F. C.; Yang, Q.; Hu, S.; Garaj, S.; Nair, R. R.; Geim, A. K.; Gopinadhan, K. Size Effect in Ion Transport through Angstrom-Scale Slits. Science 2017, 358, 511−513. (29) CRC Handbook of Chemistry and Physics, 95th ed.; Haynes, W. M., Ed.; CRC Press: Boca Raton, 2014. (30) Israelachvili, J. N. Intermolecular and Surface Forces, 3rd ed.; Academic Press: New York, 2011. (31) Zhu, Y.; Wang, F.; Bai, J.; Zeng, X. C.; Wu, H. Compression Limit of Two-Dimensional Water Constrained in Graphene Nanocapillaries. ACS Nano 2015, 9, 12197−12204. (32) Ohkubo, T.; Konishi, T.; Hattori, Y.; Kanoh, H.; Fujikawa, T.; Kaneko, K. Restricted Hydration Structures of Rb and Br Ions Confined in Slit-Shaped Carbon Nanospace. J. Am. Chem. Soc. 2002, 124, 11860−11861. (33) Luo, Z.-X.; Xing, Y.-Z.; Liu, S.; Ling, Y.-C.; Kleinhammes, A.; Wu, Y. Dehydration of Ions in Voltage-Gated Carbon Nanopores Observed by in situ NMR. J. Phys. Chem. Lett. 2015, 6, 5022−5026. (34) Tielrooij, K. J.; Garcia-Araez, N.; Bonn, M.; Bakker, H. J. Cooperativity in Ion Hydration. Science 2010, 328, 1006−1009. (35) Scheu, R.; Rankin, B. M.; Chen, Y.; Jena, K. C.; Ben-Amotz, D.; Roke, S. Charge Asymmetry at Aqueous Hydrophobic Interfaces and Hydration Shells. Angew. Chem., Int. Ed. 2014, 53, 9560−9563. (36) Rajamani, S.; Ghosh, T.; Garde, S. Size Dependent Ion Hydration, Its Asymmetry, and Convergence to Macroscopic Behavior. J. Chem. Phys. 2004, 120, 4457−4466. (37) Mukhopadhyay, A.; Fenley, A. T.; Tolokh, I. S.; Onufriev, A. V. Charge Hydration Asymmetry: The Basic Principle and How to Use It to Test and Improve Water Models. J. Phys. Chem. B 2012, 116, 9776−9783. (38) Chen, G.; Guo, Y.; Karasawa, N.; Goddard, W. A., III ElectronPhonon Interactions and Superconductivity in K3C60. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 48, 13959−13970. (39) Pan, F.; Wang, G.; Liu, L.; Chen, Y.; Zhang, Z.; Shi, X. Bending Induced Interlayer Shearing, Rippling and Kink Buckling of Multilayered Graphene Sheets. J. Mech. Phys. Solids 2019, 122, 340−363. (40) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269− 6271. (41) Joung, I. S.; Cheatham, T. E., III Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B 2008, 112, 9020−9041. (42) Wu, Y.; Aluru, N. R. Graphitic Carbon-Water Nonbonded Interaction Parameters. J. Phys. Chem. B 2013, 117, 8802−8813. (43) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (44) Hockney, R. W.; Eastwood, J. W. Computer Simulation Using Particles; Taylor & Francis: New York, 1988. (45) Aksimentiev, A.; Schulten, K. Imaging α -Hemolysin with Molecular Dynamics: Ionic Conductance, Osmotic Permeability, and the Electrostatic Potential Map. Biophys. J. 2005, 88, 3745−3761. (46) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. Comput. Phys. Commun. 2005, 167, 103−128. (47) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (48) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate ab Initio Parametrization of Density Functional Dispersion G

DOI: 10.1021/acs.jpcc.8b09742 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (49) Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703−1710. (50) Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic Separable Dual-Space Gaussian Pseudopotentials from H to Rn. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 3641−3662. (51) VandeVondele, J.; Hutter, J. Gaussian Basis Sets for Accurate Calculations on Molecular Systems in Gas and Condensed Phases. J. Chem. Phys. 2007, 127, 114105. (52) Li, Q.; Dong, Y.; Perez, D.; Martini, A.; Carpick, R. W. Speed Dependence of Atomic Stick-Slip Friction in Optimally Matched Experiments and Molecular Dynamics Simulations. Phys. Rev. Lett. 2011, 106, 126101. (53) Dong, Y.; Vadakkepatt, A.; Martini, A. Analytical Models for Atomic Friction. Tribol. Lett. 2011, 44, 367−386. (54) Koneshan, S.; Rasaiah, J. C.; Lynden-Bell, R. M.; Lee, S. H. Solvent Structure, Dynamics, and Ion Mobility in Aqueous Solutions at 25 °C. J. Phys. Chem. B 1998, 102, 4193−4204. (55) Koneshan, S.; Lynden-Bell, R. M.; Rasaiah, J. C. Friction Coefficients of Ions in Aqueous Solution at 25 °C. J. Am. Chem. Soc. 1998, 120, 12041−12050. (56) Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. CP2K: Atomistic Simulations of Condensed Matter Systems. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4, 15−25. (57) Qiu, H.; Xue, M.; Shen, C.; Guo, W. Anomalous Cation Diffusion in Salt-Doped Confined Bilayer Ice. Nanoscale 2018, 10, 8962−8968. (58) Kalluri, R. K.; Ho, T. A.; Biener, J.; Biener, M. M.; Striolo, A. Partition and Structure of Aqueous NaCl and CaCl2 Electrolytes in Carbon-Slit Electrodes. J. Phys. Chem. C 2013, 117, 13609−13619. (59) Li, P.; Song, L. F.; Merz, K. M., Jr. Systematic Parameterization of Monovalent Ions Employing the Nonbonded Model. J. Chem. Theory Comput. 2015, 11, 1645−1657. (60) Hou, D.; Li, D.; Yu, J.; Zhang, P. Insights on Capillary Adsorption of Aqueous Sodium Chloride Solution in the Nanometer Calcium Silicate Channel: A Molecular Dynamics Study. J. Phys. Chem. C 2017, 121, 13786−13797. (61) Hou, D.; Li, T.; Wang, P. Molecular Dynamics Study on the Structure and Dynamics of NaCl Solution Transport in the Nanometer Channel of CASH Gel. ACS Sustainable Chem. Eng. 2018, 6, 9498−9509. (62) Liu, J.; Shi, G.; Guo, P.; Yang, J.; Fang, H. Blockage of Water Flow in Carbon Nanotubes by Ions due to Interactions between Cations and Aromatic Rings. Phys. Rev. Lett. 2015, 115, 164502.

H

DOI: 10.1021/acs.jpcc.8b09742 J. Phys. Chem. C XXXX, XXX, XXX−XXX