Diffusion coefficient of ions through graphene nanopores

Ion transport through nanoscale confinements has attracted growing interest in recent years due to tremendous applications in many areas, including chemistry,¹ energy conversion and storage,^2–4 biology,^5–7 and desalination.⁸ Different from macroscale systems, ion transport at the nanoscale is nontrivial because of molecular interactions, which affect the dynamics of ions^9,10 and fluid properties.¹¹ As a special case of ion transport in nanoporous media, the diffusion of ions through graphene nanopores has diverse applications,¹² such as porous graphene-based electrodes^13,14 and batteries,¹⁵ owing to the excellent electrical and physical properties of graphene.¹⁶ To promote the applications of graphene in various areas, the transport of ions through graphene nanopores needs to be intensively studied.

In the past years, efforts have been made to investigate transport phenomena through graphene nanopores.^17–20 For example, the translocation of DNA molecules through graphene nanopores has been studied, which offers new opportunities for DNA sequencing^21–23 based on the current changes caused by the passage of DNA molecules through atomic-scale nanopores.²⁴ Water flows through graphene nanopores have also been examined.^20,25 It is found that the size of graphene nanopores and pore functionalization significantly affect the flow rate, whose dependence on the pressure gradient is different from the predictions of classic theories. Furthermore, ion transport through graphene nanopores has been explored.^4,25–27 It is revealed that the size and geometry, edge functionalization, and charge distribution of graphene nanopores play important roles in ion transport.^26–30 When the pore size becomes sufficiently small, approaching the diameter of ions, it requires substantial energy to overcome the electrostatic interaction between ions and graphene pores. Moreover, the dehydration of ions may be necessary for ion migration through nanopores, which hinders the transport of ions.^25,31,32 Nevertheless, most of the previous studies on ion transport are driven by external forces, such as pressure gradients and electric fields, where only the streaming velocities of ions are obtained.^26,27,33 In many applications, such as electrodes and batteries, however, the diffusion coefficient of ions, which determines ion diffusion under concentration gradients, plays an essential role in the design and fabrication of electrodes and may affect the performance of certain energy systems. For graphene nanopores, unfortunately, the diffusion coefficient of ions has not been extensively studied. The effects of pore size and porosity on the ion diffusion coefficient are unclear. Although porosity is an important parameter for ion transport in porous media, it is affected by pore size and pore–pore spacing. Therefore, the roles of both pore size and porosity require extensive investigation.

In this work, we investigate ion transport through graphene nanopores driven by concentration gradients through molecular dynamics (MD) simulations. KCl solutions are considered. The diffusion coefficients, D, of K⁺ and Cl⁻ are obtained by varying the pore size and porosity. It is revealed that the diffusion coefficients of both K⁺ and Cl⁻ strongly depend on the pore size when the pore diameter is less than 3.0 nm and increase with increasing porosity in a near-linear fashion. The transport mechanisms are examined, and a scaling law for the diffusion coefficient is developed.

The MD simulations are carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package.³⁴ The simulation system is confined by two rigid walls and partitioned by a graphene sheet with a circular pore at the center into two equal parts, which are used as reservoirs, as illustrated in Fig. 1. Initially, a KCl solution of 1.0M is filled in the left reservoir and the right reservoir is filled with water. The length of the system in the y-direction, L_y, is 20.0 nm, and the lengths in the x and z directions, L_x and L_z, range from 3.2 to 10.7 nm, depending on the pore diameter, d, and porosity, ϕ, of the graphene. Details about the system configuration are summarized in Table I. The pore in the graphene sheet is constructed by removing carbon atoms in a circular region of diameter d at the center of the graphene sheet. The pore diameter ranges from 1.0 to 4.9 nm, which corresponds to porosities from 0.2 to 0.6 (the porosity is calculated as the ratio of pore area to the area of the whole graphene sheet).

Molecular interactions are generally described by the Lennard-Jones (LJ) 12-6 and Coulomb potentials,³⁵ 

where r_ij is the distance between atoms/ions i and j, q_i and q_j are the charges carried by the interacting atoms/ions, and ε_ij and σ_ij are the binding energy and collision diameter of the LJ potential. Water molecules are modeled by the extended simple point charge (SPC/E) water model.^36,37 The SHAKE algorithm is applied to maintain water molecules to be rigid.³⁸ The graphene sheets considered are neutral and rigid (carbon atoms are fixed). The LJ parameters and charges for various atoms/ions are summarized in Table II. The LJ parameters for different interacting species are obtained using the Lorentz–Berthelot mixing rules.^39,40 The particle–particle particle–mesh (PPPM) method is used for long-range Coulombic interactions. The cutoff distance for the LJ potential and electrostatic interactions is set as 1.0 nm. The left wall is fixed such that the volume of the left reservoir is constant. The pressure of the right reservoir is maintained at 1 atm by applying a force on the right wall. Periodic boundary conditions are employed in the x and z directions. The time step is set as 2.0 fs. Simulations are performed in the canonical (N, V, T) ensemble at T = 300 K, which is controlled by the Nosé–Hoover thermostat. Initially, the graphene sheet is complete and the system is relaxed for 3 ns. Then, a pore is constructed at the center of the graphene, which allows ions to diffuse to the right reservoir. Data collection is performed for 30–100 ns, depending on the pore size and porosity. To obtain the average properties, at least three independent simulations with different initial conditions are conducted for each case.

To understand the mechanism of ion transport through pores, the potential of mean force (PMF) of ions is computed for which the simulation system is evenly split into many boxes in the y direction with an ion constrained at the center of each box using a harmonic potential. The PMF is obtained based on the positions of ions through umbrella sampling with the combination of the weighted histogram analysis method (WHAM).^41,42 This method has been used and refined for studying ion permeation through porous media.^27,43,44 The PMF is calculated from y = −10 to 10 Å with 40 sampling boxes. The time step for the PMF calculation is set as 1 fs. Data are collected for 2.0 ns after 100 ps equilibration for each box, which corresponds to a total simulation time of 84 ns.

One way to determine the diffusion coefficient of ions, D, is to use Fick’s first law of diffusion, which is given by

where J is the flux of ions and ∇c is the ion concentration gradient. Usually, a constant ∇c needs to be maintained to obtain D directly. Numerically, however, this is nontrivial because the ion concentrations in the two reservoirs keep changing as ions diffuse from high to low concentration and the flux J also varies with time. The time-dependent format of Eq. (2) can be written as

where t is time, $N′t$ is the number of ions diffusing through the pore per unit time at t = t, and N(t) is the number of ions transporting through the pore during the time interval, $t∈0,t$⁠. N₀ is the initial number of ions in the left reservoir, and ∆y is the thickness of the graphene. Figure 2 shows N(t) as a function of time for pores of various diameters. It is seen that N(t) increases with time at the beginning due to large concentration gradients. As the concentration difference between the two reservoirs decreases, N(t) gradually levels off and asymptotically approaches a constant value, for which the ion concentrations in the two reservoirs become the same. The time variation of N(t) and the boundary conditions, N(0) = 0 and N(∞) = N₀/2, suggest that N(t) can be expressed as

where a is a positive parameter. As $N′t$ = dN(t)/dt, $N′t$ = aN₀e^−at/2, which, together with Eq. (4), give

Using Eq. (5), with a being determined by fitting the data of N(t), the diffusion coefficient is calculated, which is shown in Fig. 3 as a function of pore diameter d and porosity ϕ. In Figs. 3(a) and 3(b), it is seen that the diffusion coefficients for both K⁺ and Cl⁻ increase with increasing pore diameter when d < 3 nm and roughly remain constant for d ≥ 3 nm, regardless of the value of porosity. For a given pore size, the dependence of the diffusion coefficient on the porosity shows an approximately linear increase with increasing porosity, as illustrated in Figs. 3(c) and 3(d).

The transport of ions through nanopores is governed by several factors. An important one is the energy barrier, ∆G, which regulates the probability of ions, P, migrating through nanopores via P ∼ exp(−∆G/k_BT), where k_B is the Boltzmann constant,⁴⁵ i.e., a large ∆G corresponds to a low diffusion coefficient. ∆G can be obtained by calculating the PMF of ions and is determined as the difference between the minimum and maximum of PMF next to and at the pore, respectively. As an example, the PMF and ∆G for the pore with d = 1.0 nm for K⁺ are shown in the inset of Fig. 4(a). Figure 4 shows the energy barrier ∆G for K⁺ and Cl⁻ as a function of pore diameter [Figs. 4(a) and 4(b)] and porosity [Figs. 4(c) and 4(d)]. In Figs. 4(a) and 4(b), it is seen that ∆G decreases with increasing pore size, indicating that the transport of ions through nanopores is relatively easy for large pores, which is expected. It is also seen that ∆G decreases more rapidly for pores with d < 3 nm compared with relatively large pores (d ≥ 3 nm). For a given pore size, however, ∆G is not very sensitive to the porosity, as shown in Figs. 4(c) and 4(d). As the porosity increases, ∆G decreases mildly, especially for large pores. ∆G in Fig. 4 is consistent with and explains the variation of the diffusion coefficient in Fig. 3.

In addition to the ion–pore interaction, which is the major factor determining the PMF, the dehydration of ions also affects the PMF, especially for small pores. It is known that an ion in a solution is usually wrapped by layers of water molecules due to strong ion–water interactions, which are termed the first and the second hydration shells, as illustrated in the insets of Figs. 5(a) and 5(b) (certain ions may only have one hydration shell, depending on the hydration energy). The numbers of water molecules in these hydration shells are called the coordination number, C₁ and C₂, respectively. Herein, we examine the total number of water molecules in the two hydration shells and denote it as the total coordination number, C, i.e., C = C₁ + C₂. As an ion diffuses through a small nanopore, some water molecules in the hydration shells need to be removed, which is called dehydration. A dehydration process associates with the decrease of the total coordination number and requires certain dehydration energy. Therefore, it hinders the transport of ions. Previous studies show that dehydration is the key contribution to the energy barrier when ions diffuse through narrow pores.^25,46,47

To understand the dehydration effect, the total coordination number of ions is studied. Figure 5 plots the total coordination number of K⁺ and Cl⁻, C, at the pores for various pore diameters and porosities. In Figs. 5(a) and 5(b), it is seen that C is smaller than the bulk value (solid lines), which is 26.7 and 30 for K⁺ and Cl⁻, respectively, for small pores. For the pore with d = 1.0 nm, the values of C indicate that around six to seven water molecules are removed from the hydration layers when K⁺ and Cl⁻ migrate through the pore. This makes ion transport difficult and, consequently, reduces the diffusion coefficient. As the pore size is increased, C approaches the bulk value and remains almost constant, suggesting that the dehydration effect is unimportant. For a given pore size, the values of C in Figs. 5(c) and 5(d) show that the porosity plays a minor role in the dehydration of ions compared with the pore size. The variations of C in Fig. 5 show that the dehydration only takes place in small pores, where it contributes to the energy barrier. For large pores, the energy barrier is mainly caused by ion–pore intermolecular interactions.

The results for the diffusion coefficient in Fig. 3 can be further analyzed. For the dependence of the diffusion coefficient on the pore diameter, it is easy to find that d = 3.0 nm appears to be a critical pore diameter, below and above which the diffusion coefficient behaves differently. For d < 3.0, the variation of D in Figs. 3(a) and 3(b) suggest an exponential relationship, $D∼1−exp−d2$⁠, where the term d² is employed to reflect that ion diffusion is regulated by the pore area. For d ≥ 3.0 nm, the pore size seems to be unimportant and the diffusion coefficient mainly depends on the porosity. The fashion of D in Figs. 3(c) and 3(d) indicates a near-linear relationship between D and porosity ϕ. If dimensionless quantities are used, a general scaling law for the diffusion coefficient can be proposed as

where D₀ is the bulk diffusion coefficient (D₀ = 0.1048 and 0.1043 m² s⁻¹ for K⁺ and Cl⁻, respectively, calculated by using a system without graphene), λ, α, and β are fitting parameters, and $d*=d−d0/dion$ with d₀ and d_ion being the minimum pore diameter allowing ion transport and the diameter of ions (d_ion = 0.284 and 0.483 nm for K⁺ and Cl⁻, respectively). Herein, d₀ is determined as the diameter of ions plus the collision diameter of the LJ potential for ion–carbon interactions (d₀ = 0.5955 and 0.894 nm for K⁺ and Cl⁻, respectively).

The best fit between Eq. (6) and the results in Fig. 3 gives λ = 0.058, α = 0.812, β = 0.101 for K⁺ and λ = 0.321, α = 0.804, β = 0.105 for Cl⁻. Figure 6 shows D/(D₀ϕ) as a function of d* for d < 3.0 nm and D/D₀ as a function of ϕ for d ≥ 3.0 nm. It is seen that all the data collapse onto master curves. Therefore, Eq. (6) can be used as a scaling law for predicting the diffusion coefficient. It is noted that more accurate modeling would be achieved using a power term of ϕ, i.e., D/D₀ ∼ ϕ^γ, for which γ ≈ 1. To make the scaling law simple, a linear dependence on ϕ is employed in Eq. (6) and the error is less than 3.5%.

Finally, it is worth mentioning that the results in Fig. 3 are only valid for the transport of K⁺ and Cl⁻ through graphene pores. The general format of the scaling law in Eq. (6), however, may apply to other monovalent ions, such as Na⁺. For other ions (e.g., bivalent ions), the diffusion coefficient requires further investigation.

In summary, we have studied the diffusion of K⁺ and Cl⁻ through graphene nanopores for various pore sizes and porosities. The diffusion coefficient is sensitive to the pore size for small pores. For relatively large pores (d ≥ 3 nm), the diffusion coefficient becomes roughly independent of the pore size. The dependence of the diffusion coefficient on the porosity is nearly linear. The diffusion process is determined by the free energy barrier at the pore, which is caused by the ion–graphene molecular interactions. In addition, the dehydration of ions at the pore also hinders the transport of ions, which tends to reduce the diffusion coefficient and is important for small pores. A general scaling law has also been developed to model the dependence of the diffusion coefficient on pore size and porosity.

This work was partially supported by the Research Grants Council of the Hong Kong Special Administrative Region under Grant No. 16209119 and the project of Hetao Shenzhen-Hong Kong Science and Technology Innovation Cooperation Zone (Grant No. HZQB-KCZYB-2020083).

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article.