Introduction

Carbon-based surfaces widely exist in both nanoscale and macroscopic materials and most of them contain π electron-rich structures which is a hexagonal carbon ring where there are rich π electrons, such as aromatic rings in biomolecules1, humus in the soil2,3, polycyclic aromatic hydrocarbons in air pollutants4,5,6 and graphene7,8, carbon nanotubes9,10 and fullerenes11. It has been widely reported that the accumulation of ions at carbon-based surfaces greatly impacts the properties and applications of carbon-based materials, including surface- or interface-charged behaviors12,13,14,15, aqueous suspension and aggregation16,17, surface modification18,19, chemical synthesis20,21 and the function, configuration and stability of biological systems22,23,24. We note that, in nanoscale systems, since the total number of ions is small, the accumulation of ions at a carbon-based surface results in a clear decrease of ions in the solution and at the liquid-gas interface. This affects impact the chemical reactions and charged behaviors in the solution or at the gas-liquid interface on aerosol particles25,26,27,28,29,30,31,32, water films on the soils in the dry environment and aggregations of polymers/macromolecules/nanoparticles in humid environment.

The understanding of the ions accumulation on the surfaces with π electron-rich structures becomes more complex due to the discovery of the cation-π interactions in 1980's33. The cation-π interaction is the non-covalent interaction between a cation and a π electron-rich carbon-based structure. Basing on experiments and quantum mechanical methods computations, it has been found that cation-π interactions play important roles in various systems, including control of the structures and functions of microscale and nanoscale materials, macromolecules and proteins34,35,36,37,38,39,40,41,42,43,44.

Quantum mechanical (DFT and ab initio) methods usually can only compute and simulate small systems with a small number of atoms1,34,35,36,37,38. For the complex systems such as the ions in the solution, we usually use classical analysis and numerical simulations with appropriate force fields. Unfortunately, the cation-π interactions have rarely been explicitly and sufficiently considered in usual all-atoms classic force fields. In the biological systems where the π electron-rich structures (for example aromatic rings) are enriched, in present popular force fields (Amber45, Charmm46 and Opls47), cation-π interactions are mainly considered by increasing the charges of the carbon atoms in π electron-rich structures. This simple treatment leads to a stronger rejection to the anions by the carbon atoms and enhanced interactions between the carbon atoms and surrounding water molecules and thus, affects the configuration and aggregation of protein in the aqueous salt solution22,23,24.

To the worse, although many of the most well-known nano-materials (such as graphene, nanotubes and fullerenes) include the π electron-rich carbon-based surfaces, cation-π interactions have not been included explicitly when considering behaviors of interfacial ions on these nano-materials48,49,50,51. In popular classic force fields such as Amber, Charmm and Opls, the interactions between the ions and these π electron-rich carbon-based surfaces only include the van der Waals interaction; it is undistinguishable for the interactions between ions and carbon-based surfaces, whether the surfaces contain π electron-rich structures or not. This greatly underestimates the interactions energies between the ions and π electron-rich the carbon-based surfaces. For example, the interaction energies of a Na+ in the solution with the graphene surface are −0.245, −0.946 and −0.247 kcal/mol from the classic force fields Amber, Charmm and Opls, respectively. These values are significantly lower than the value of −16.4 kcal/mol for hydrated Na+ in calculations with the density functional theory (DFT) methods52. We note that, very recently, the computation of the interactions between the π electron-rich carbon-based surfaces (graphene, nanotubes and fullerenes) and a single ion in vacuum has been greatly improved by modifying a polarizable force field53. Due to the hydration of the ions by water molecules, this force field cannot be directly used for the study of ions at such liquid-solid interface.

In this article, based on the DFT computations, we incorporate cation-π interactions to the classic force field to study the behavior of interfacial ions on the surfaces with π electron-rich carbon-based structures in the aqueous solution. On the basis of the example of the uncharged hydrophobic graphite surface and the NaCl solutions, we show that large amounts of ions in the aqueous solution are adsorbed on the carbon-based surface with π electron-rich structures, when the cation-π interaction is explicitly account for. Interestingly, although the interaction of the hydrated Cl− with the carbon-based surface is small, Cl− is also enriched on the surface due to the electrostatic interaction between Na+ and Cl−. This difference of Na+ and Cl− interaction strength with carbon-based surface due to the cation-π interactions leads to more Cl− in the solution, resulting in a significant negatively charged behavior of solution (especially for nanoscale systems). We also show that the enrichment of cations on the graphite surface is universal and that other cations (Li+, K+, Mg2+, Ca2+, Fe2+, Co2+, Cu2+, Cd2+, Cr2+ and Pb2+) have similar behavior as Na+. This finding provides new insight for understanding and controlling surface-charged behaviors and points to additional applications such as ion storage and detection, carbon nanotube-based water purification, graphene-based ion filtration and the development of low-cost carbon-based materials for water and soil pollution from heavy metal cations.

Results and discussion

Due to strong Na+-π interactions, the hydrated Na+ still significantly interact with the graphite surface, although the hydration interaction reduces the Na+-π interactions52. For example, the interaction energy of a Na+ around nine water molecules at its equilibrium position, including the Na+-π interaction, is 16.4 kcal/mol (the detailed computation can be found in the Supplementary Information). This is much larger than the value of 0.9 kcal/mol obtained if the Na+-π interaction is neglected.

In order to quantitatively study the effect of the cation-π interaction on the ion distributions on the graphite surface, using density functional theory (DFT) (see details in the Supplementary Information), we modified the NAMD program by introducing the Na+-π interactions. We modeled the interaction as

where ε and zm are the adsorption energy and balance position distance (the vertical dimension between the Na+ and the surface) of Na+ on the graphite surface, respectively and z is the distance of vertical dimension between the Na+ and the surface. We note that the hydrated Cl−-π interaction is much smaller, which is only of ~1/10 of the hydrated Na+-π interaction, so it is neglected in the computation.

We first consider how the hydrated Na+-π interaction affects the distribution of Na+ and Cl− in the NaCl solution. Figure 1A and 1B shows a typical snapshot in which many Na+ and Cl− in the NaCl solution are adsorbed on the graphite surface and others are free in the bulk water. In order to obtain a quantitative description, we compute the distribution probabilities of Na+ and Cl− in the water along the z-direction on the graphite surface, which are shown in Figure 1C. We can see that there is a peak for Na+ at ~0.39 nm. We also compute the average number of Na+ at the graphite surface, which is defined as the number of Na+ with z positions less than zC; zC = 0.50 nm is the location of the first valley in the Na+ distribution, with a value of 48.0. Considering that the total number of Na+ in the solution is only 70, it is clear that most of Na+ is adsorbed on the graphite surface. For comparison, we also performed MD simulations with the hydrated Na+-π interactions turned off. As shown in Figure 1C, without the hydrated Na+-π interactions, Na+ has a flat distribution along the z-direction.

Figure 1
figure 1

Na+ and Cl− in the NaCl solution on the graphite surface.

(A) Snapshot. (B) Close-up of the lower left corner. The gray structures depict the graphite sheets; water molecules and ions are shown with oxygen in red, hydrogen in white, Na+ in blue and Cl− in green, respectively. (C and D) Distribution probabilities of Na+ (C) and Cl− (D) along the z-direction in the NaCl solution on the graphite surface, with the normal hydrated Na+-π interaction (ε = ε0, red solid line), reduced interaction of ε = ε0/2 (green solid line) and ε = ε0/4 (blue solid line) and without the hydrated Na+-π interaction (ε = 0.0, black dashed line).

Full size image

Interestingly, there is also a peak for Cl− even when the hydrated Cl−-π interaction is neglected. Moreover, the peak for Cl− is located at z = 0.37 nm, which is lower than the location of the peak for Na+, as shown in Figure 1D. Careful examination shows that the Na+-Cl− electrostatic interaction together with the van der Waals interaction between Cl− and the graphite surface results in a large amount of Cl− preferring to remain closer to the graphite surface. As with the two typical cases shown in Figure 2A and 2B, a Cl− close to the graphite surface always has a very close neighboring Na+ and the nearby water molecules also provide an indirect interaction between Cl− and Na+. Our analysis shows that the average distance between a Cl− close to the graphite surface and the nearest Na+ is only 0.26 nm, which is much smaller than the average distance (0.34 nm) between a Cl− in the solution and its nearest Na+. However, the z value of the peak for Cl− should be larger than the z value of the peak for Na+ if the Cl− close to the solid surface only results from the electrostatic attraction by Na+. We note that the Cl− have a much stronger van der Waals interaction with the graphite surface compared to Na+, since the Lennard-Jones (LJ) parameters are εCCl = 0.103 kcal/mol and εCNa = 0.057 kcal/mol, respectively. Thus, we postulate that it is the combination of the van de Waals interaction of Cl− with the graphite surface and the Na+-Cl− electrostatic interaction that makes the Cl− accumulate closer to the graphite surface. The van der Waals distance σGCl of 0.37 nm, consistent with the z value of the peak for Cl−, further verifies this assumption. Since the “apparent” adsorption of Cl− by the graphite surface is mediated by the Na+-Cl− electrostatic interaction, the peak for Na+ is higher than the peak for Cl−. Moreover, the average number of Cl− adsorbed at the graphite surface, defined as the number of Cl− with z positions less than zC, is only 40.3 and is smaller than the corresponding value of 48.0 for Na+. This leads to the solution, especially near the liquid-gas surface, having negatively charged behavior since there is more Cl− floated in the solution. As shown in Figure 1D, without the hydrated Na+-π interactions, Cl− also has only a flat distribution along the z-direction.

Figure 2
figure 2

Ions at the liquid-graphite interface.

(A and B) Snapshots of two typical cases of a Cl− with its neighboring Na+ and water molecules on the graphite surface. In (A, top view), a water molecule (marked by the red arrow) directly interacts with both Na+ and Cl− simultaneously (the interactions are shown by blue and green dotted lines). In (B, top view), one water molecule (marked by the blue arrow and blue dotted line) interacts with Na+ and another water molecule (marked by the green arrow and green dotted line) interacts with Cl−. These two water molecules form a hydrogen bond (marked by the red dotted line) to bridge the indirect interaction between Na+ and Cl−. In the side view of (A,B), the horizontal line of the center Na+ (blue) and Cl− (green) are shown. The color settings are the same as in Figure 1. (C) The numbers of ions adsorbed on the graphite surface (black) and floated in the solution (green with right coordinates) with respect to the total numbers (N) of ions in the aqueous salt solution.

Full size image

We also performed simulations for systems with different amounts of NaCl solution at the same concentration. Figure 2C shows the average numbers of Na+ and Cl− adsorbed on the graphite surface and floated in the solution. We can see that the numbers of Na+ and Cl− adsorbed on the surface are saturated with a total number of NaCl, N = 70. It should be noted that all Na+ are adsorbed on the surface, but some Cl− still float in the solution, N = 30; this results in a significant negatively charged behavior in the solution and gas-solution interface. The numbers of Na+ and Cl− floated in the solution increase as N increases. From N = 110, the amounts of Na+ and Cl− floated in the solution are larger than those adsorbed on the surface. This indicates that the accumulation of ions on the surface will largely reduce the ion distribution in the solution and gas-liquid interface, especially for small nanoscale systems.

The adsorption of Na+ and Cl− on the graphite surface in the NaCl solution is robust. We performed new simulations with smaller values of ε, i.e., ε = ε0/2 and ε = ε0/4. As shown in Figure 1C and 1D, the peaks are still clear for both values of ε; the peak heights only decrease about one third value for ε = ε0/2. For Na+, the height of the peak only decreases a little when the values of ε decrease from ε0 to ε0/2. Correspondingly, the height of the distribution probability of Na+ increases slightly over the graphite surface distance z > 0.5 nm. We can see a significant decrease in the peak height for ε = ε0/4 and a clear increase in the distribution probabilities for z > 0.42 nm. For Cl−, as shown in Figure 1D, we can see a little change in the peak when the values of ε decrease from ε0 to ε0/2. The peak height decreases to its quarter value for ε = ε0/4. In other words, Na+ and Cl− are still abundant on the graphite surface when the adsorption energy ε is reduced to half ε0, even one-quarter ε0.

The adsorption of cations on the graphite surface is universal, in that other cations have similar behavior as Na+. We performed DFT computations of the other cations (Li+, K+, Mg2+, Ca2+, Fe2+, Co2+, Cu2+, Cd2+, Cr2+ and Pb2+) on the graphite surface, as shown in Figure 3. We find that the adsorption energy of K+ is the smallest, at 28.7 kcal/mol and still reaches 48.1 kBT at T = 300 K, indicating that all of these metal ions can be adsorbed on the surface at room temperature. We note that the adsorption stability of divalent cations is significantly stronger than that of monovalent cations, so we expect that the graphite surface has a special adsorbing capability for heavy cations, such as Fe2+, Co2+, Cu2+, Cd2+ and Cr2+, which are common heavy metal pollutants in soil and water.

Figure 3
figure 3

Adsorption energies of various cations on the graphite surface.

Full size image

In summary, with a careful consideration of the cation-π interaction, we conclude that Na+ is enrichment on the typical hydrophobic carbon-based surface with the π electron-rich structures in aqueous salt solutions. Although the interaction of hydrated Cl− with the carbon-based surface is weak, the electrostatic Na+-Cl− interaction also enriches Cl− at the interface. The average number of Na+ adsorbed at the graphite surface is significantly higher than that of Cl−, which results in more Cl− floated in the solution. We emphasize that, in nanoscale systems, a small total number of ions leads to a significant negatively charged behavior in the solution and enhances the negative charge at the liquid-gas interface.

These observations are robust in that similar results are obtained even when the strength of the hydrated Na+-π interaction is reduced to half of the original values in the numerical simulations. We note that other cations (Li+, K+, Mg2+, Ca2+, Fe2+, Co2+, Cu2+, Cd2+, Cr2+ and Pb2+) have similar behavior to Na+, in that the cations adsorb on the carbon-based surface with π electron-rich structures.

These findings provide new insight into the understanding and control of surface-charged behaviors and they may lead to potential applications such as ion storage and detection, carbon nanotube-based water purification, graphene-based ion filtration and bio-related processes. Finally, we note that graphite surfaces have a wonderful adsorbing capability for heavy cation pollution, such as Fe2+, Co2+, Cu2+, Cd2+ and Cr2+, indicating that a low-cost carbon-based material with π electron-rich structures may have widespread potential for dealing with water or soil pollution from heavy metals.

Methods

The B3LYP54,55 method in the framework of DFT, which have been widely used in the study of water molecules on solid surfaces or inside carbon nanotubes52,56,57, is used to study the intermolecular interactions. For geometry optimizations, the double-ζ basis is employed and a d-polarization function is added (marked with 6–31 G(d)). At the same time, the pseudopotential function with Lanl2dz is introduced into the basis set for some metal ions (Fe2+, Co2+, Cu2+, Cd2+, Cr2+ and Pb2+). A two-dimensional graphite surface of 12.275 × 15.658 Å2 is used for the DFT study as shown in Figure 4A, which is large enough to obtain results with a tolerable error, whereas the detailed discussions are shown in our former article52. All calculations are carried out using the Gaussian-03 program (Revision D.01, details in Supplementary Information).

Figure 4
figure 4

(A) Schematic description of a two-dimensional graphite surface (C84H24) of 12.275 × 15.658 Å2 (84 carbon atoms and 24 hydrogen atoms). (B) The interaction potential between Na+ and graphite surface (red line) in the present well-used force field (CHARMM) without including Na+-π interactions. The adsorption energies of Na+ with nine water molecules on the graphite surface (black square) on the different distance (z, the vertical dimension between the Na+ and the surface) at the B3LYP/6–31 G(d) level and the fitting (blue line) potential. The interaction potential between Na+ and graphite surface (green line) in our used force field incorporated Na+-π interactions to CHARMM force field. (C) Hydrated Na+-π interactions energies ε (black squares) and fitting curve (pink line) and the H2O-π interactions energies ε (red filled circles).

Full size image

The cation-π interactions between Na+ and the graphite surface are represented by a model potential in the classical MD simulation. Here, we fitted the overall profile a Na+ with nine water molecules separating the nonpolar hydrophobic graphite surface, as shown in Figure 4B. The fitted potential function presents in Equation 1, where the parameters ε and zm are the adsorption energy and balance position distance (the vertical dimension between the Na+ and the surface) of Na+ with the graphite surface and z is the distance of the vertical dimension between the Na+ and the surface, which are the main potential parameters describing the cation-π interaction. The zm = 3.8 Å and ε = ε0 = −16.4 kcal/mol are used, corresponding to the case of a Na+ with nine water molecules adsorbed on the graphite surface (see Figure 4C, detailed computation can be found in the Supplementary Information). These cation-π interactions of Na+ with the nonpolar hydrophobic graphite surfaces were added in MD dynamics simulations programs. Figure 4B shows the present well-used force field (CHARMM) without including Na+-π interactions, our used force field incorporated Na+-π interactions to CHARMM force field and the DFT computational results as well as its fitting curve. The DFT computational results are much larger than the present well-used force field obtained if the Na+-π interaction is neglected and our used force field incorporated Na+-π interactions to CHARMM force field is consistent with the DFT computational results.

Molecular dynamics simulations are carried out using the program NAMD2/VMD1.9 packages58, with the CHARMM force field46, at time steps of 2 fs with the O-H bonds and C atoms held fixed. The TIP3P water model59 is used with parameters of εOO = 0.152 kcal/mol, εOH = 0.084 kcal/mol, εHH = 0.046 kcal/mol, σOO = 3.15 Å, σOH = 1.77 Å and σHH = 0.40 Å, respectively. The carbon atoms are modeled as uncharged Lennard-Jones particles with parameters of σCC = 3.55 Å and εCC = 0.070 kcal mol−1. Na+ and Cl− ions are assigned charge to 1.0 e and −1.0 e, respectively, with Lennard-Jones parameters of σNaNa = 2.43 Å, σClCl = 4.04 Å, εNaNa = 0.047 kcal mol−1 and εClCl = 0.150 kcal mol−1. The Lennard-Jones parameters of water-ions are εONa = 0.084 kcal/mol, εOCl = 0.151 kcal/mol, εHNa = 0.046 kcal/mol, εHCl = 0.083 kcal/mol, σONa = 2.79 Å, σOCl = 3.60 Å, σHNa = 1.41 Å and σHCl = 2.22 Å, respectively. In order to speed up computation, the particle-mesh Ewald method with a grid spacing on 1 Å or less is used to full electrostatic interactions computation, while the van der Waals (vdW) interactions are shifted smoothly at 12.0 Å. and constant-volume molecular dynamics simulation is performed. The simulations are performed at NVT with a constant-temperature (300 K) and a Langevin damping coefficient 5 ps−1. The selection of a vapor-liquid coexistence system is used to maintain the ambient condition. Periodic boundary conditions are applied in all directions.

For the NaCl solution, the initial simulation box size is Lx = 4.2 nm, Ly = 4.2 nm and Lz = 20.000 nm, where Lz is set to be sufficiently large to eliminate the image effect in the z direction. During the simulations, we fix the graphite surface, which contains two layers graphite and there are 598 carbon atoms of every layer and every carbon-carbon bond is 1.420 Å. All of Na+, Cl− and water molecules are free to move in this simulation. All systems are equilibrated for 4 ns.