Capturing charge and size effects of ions at the graphene–electrolyte interface using polarizable force field simulations

Abstract

We present a systematic investigation capturing the charge and size effects of ions interacting with a graphene surface using polarizable simulations. Our results utilizing the Drude polarizable force field (FF) for ions, water and graphene surfaces, show that the graphene parameters previously developed by us are able to accurately capture the dynamics at the electrolyte–graphene interface. For monovalent ions, with increasing size, the solvation shell plays a crucial role in controlling the ion–graphene interactions. Smaller monovalent ions directly interact with the graphene surface, while larger ions interact with the graphene surface via a well-formed solvation shell. For divalent ions, both interaction modes are observed. For the anion Cl⁻, we observe direct interaction between the ions and the graphene surface. The anion–graphene interactions are strongly driven by the polarizability of the graphene surface. These effects are not captured in the absence of polarization by additive FF simulations. The present study underlines the importance of polarizability in capturing the interfacial phenomenon at the solid–solute interface.

---

A Introduction

Graphene, a 2D allotrope of carbon with a hexagonal unit cell and layered architecture has attracted significant attention from the scientific community with applications in varied fields such as desalination,^1–5 energy storage^6–9 and electrochemical sensing.^10–12 For example, by tuning the interlayer separation in graphene-oxide (GO) membranes Abraham et al. achieved near-perfect salt-rejection, establishing the applicability of graphene and graphene-oxide based membranes for desalination.^13,14 For technological advancement in these applications, it is important to understand the adsorption behaviour of the constituent ions at the electrolyte–graphene interface. It is to be noted that simple theories based on the continuum model, which consider only the size and the magnitude of the charge, cannot capture the dynamics of the ion–surface interactions. The shortcomings of the continuum model are even more pronounced for hydrophobic surfaces as observed in graphene and GO sheets, where the interactions are not purely electrostatic.¹⁵ Experimental studies using deep UV second harmonic generation measurements found that direct ion–graphene interactions were responsible for the adsorption of SCN⁻ ions on the graphene surface, with the interactions being enthalpically driven.¹⁶ A simple continuum model fails to capture these effects. It has been noted that the adsorption of ions on hydrophobic surfaces is affected by various properties of the ions, including their size, polarizability and hydration free energy to name a few.^17,18 Gaining molecular insight into the governing mechanisms of ion-adsorption requires accurate modelling of the electrolyte–graphene interface.

QM calculations can accurately predict such properties, but the scaling of O(N^3–4) for DFT calculations limits the applicability of such methods to very small systems.¹⁹ In one such study the interaction of hydrated cation clusters (cation-(H₂0)₇) with a graphene sheet was studied which revealed reorganization of the hydration shell around K⁺, allowing K⁺ to directly interact with the graphene sheet, while such reorganization of the hydration shell was not observed for Li⁺.¹⁹ These static QM calculations were followed by very short (10 ps) AIMD simulations to comment on the hydration dynamics.^19,20 These studies highlight the need for alternative methods to accurately capture the ion–graphene interaction in aqueous environments.

Molecular dynamics (MD) simulations offer a viable solution to such problems. However, MD simulations based on classical additive force fields (FFs) cannot capture the ion–graphene interactions as they do not account for charge transfer. Gas-phase and polarizable continuum model (PCM) calculations have shown that the polarizability of ions and the surface plays an important role in describing the strength and evolution of the ion–graphene interactions.^21,22 Based on these facts efforts have been focused towards the development of FF parameters that can describe the interactions between ions and the graphene surface.^23,24 In these studies, an additional force field term is added to describe the ion–graphene interactions. The parameters are empirically tuned to reproduce the QM calculated binding energies, thus enabling additive simulations to capture the ion–graphene interactions in an otherwise non-polarizable FF simulation.^23,24 One of the major drawbacks of such implementations is that the target DFT calculations used to tune the ion–graphene interactions over-emphasize the binding energetics as they are based on a mean-field description of the solvation rather than explicit solvation.

To this end polarizable simulations offer a framework for studying the evolution of ion–graphene interactions in explicit solvent, negating the need to empirically tune the FF parameters.^25,26 In our earlier work we described the development and testing of Drude parameters to describe multilayer and monolayer graphene surfaces.^26,27 Here, we study the energetics and dynamics of an aqueous salt solution – multilayer graphene system using Drude polarizable FF simulations, to critically evaluate the ability of FF parameters in describing the interfacial phenomenon.

B Computational methodology

Multilayer graphene comprising four graphene sheets with a radius of 15 Å and combined height of 10.04 Å, where height corresponds to the distance between the lower graphene sheet and the top-most graphene sheet, was modelled using the Inorganic Builder plugin available in the Visual Molecular Dynamics (VMD) package.²⁸ The graphene sheets were stacked in ‘ABAB’ arrangement to simulate multilayer graphene. The multilayer graphene was solvated to bring the final system dimensions to 29.42 × 29.42 × 45 ų. We introduced one molecule (one mole) of the salt of interest (NaCl, KCl, MgCl₂, CaCl₂ and CsCl) to the solvated system using the AutoIonize plugin available in VMD.²⁸ We present a representative image of the system setup in Fig. S1.† The Chemistry at Harvard Molecular Mechanics (CHARMM36) all-atom force field (FF) was used to describe bonded and non-bonded interactions in graphene sheets and ions.²⁹ A TIP3P three-point water model³⁰ was used to describe water molecules in additive FF simulations.

We employed in-house scripts to add Drude particles to the final systems to generate the corresponding Drude FF files. The classical Drude oscillator polarizable FF hereto referred to as the Drude polarizable FF was used to describe non-bonded interactions in ions.^25,26 Drude polarizable FF parameters previously tested by us were used to describe bonded and non-bonded interactions in graphene sheets.²⁶ The SWM4-NDP polarizable water model³¹ was used to describe water molecules in Drude polarizable FF simulations. We restrained the bond lengths and bond angles in water molecules using the SETTLE algorithm³² in Drude polarizable FF simulations.

Molecular dynamics simulations were performed using the nanoscale molecular dynamics^33,34 (NAMD) package in an isobaric-isothermal (NPT) ensemble. Particle mesh Ewald³⁵ (PME) summation was used to evaluate the electrostatic interactions with a real-space cut-off of 9 Å. All simulations were performed under NTP conditions (298 K and 1 atm pressure) with Langevin dynamics and the Nose–Hoover Langevin piston to maintain the NTP conditions. We employed an additional dual thermostat at 1 K to maintain the Drude particles at 1 K during Drude polarizable FF simulations. All systems were minimized for 4000 steps and equilibrated for 1 ns in NPT and 1 ns in NVT ensembles respectively. The central atom of the graphene sheet was restrained using a harmonic potential of 1 kcal mol⁻¹ Å⁻² to arrest the sliding motion of the graphene sheets. No additional restraints were employed during the simulations, and the sheets were allowed to relax and breathe during the simulations. We had previously shown that constraining the graphene surface affects the energetics of the system.²⁶ The equations of motion were integrated using a time step of 1 fs in additive FF simulations and 0.5 fs in Drude polarizable FF simulations. Production simulations were performed for a total of 100 ns for each of the six salt–graphene systems considered.

We performed simulations using adaptive biasing force (ABF) to estimate the binding free energies for ion–graphene interactions. ABF scans were performed in a window spanning 3 Å to 15 Å above the topmost graphene sheet. A biasing force of 0.20 kcal mol⁻¹ Å⁻² was employed at both ends of the scan length. The scan window was divided into bins of 0.05 angstrom length. The ABF runs were performed in a hexagonal cell of dimensions 29.42 × 29.42 × 45 ų. ABF runs were performed for a total of 100 ns, and the runs were deemed convergent when the number of observations in each bin was >1.5 × 10⁶. The methodology for the calculation of ABF was adopted from the work by Comer et al.^36,37

The densities of ions and water at the graphene–ion interfaces were calculated using in-house scripts. Densities were computed in a window spanning 15 Å from the surface of the multilayer graphene. The calculation window was divided into bins of 0.1 Å length, and the number of observations in each bin was computed.

C Results and discussion

Before discussing the characteristics of the ion–graphene interactions at the graphene–solvent interface we first establish the energetics of ion–graphene interactions. To this end, we calculate the binding free energy for Li⁺, Na⁺, K⁺, Cs⁺, Mg²⁺, Ca²⁺ and Cl⁻ using ABF (adaptive biasing force) simulations for both polarizable and additive force fields. We chose four monovalent cations (Li⁺, Na⁺, K⁺, and Cs⁺) and two divalent cations (Mg²⁺ and Ca²⁺) to evaluate the influence of size and charge on the binding free energetics. Cl⁻ is used as the common counterion. For the ABF simulations we use one molecule of the salt of interest in the solution. Representative images depicting the system setup for ABF simulations are presented in Fig. S1(a) of the ESI file.† We employed the z-axis projection of the center-of-mass distance between the ions of interest and the top layer of the multilayer graphene as the scan coordinate. The scan coordinate used for the ABF simulations is presented in Fig. S1(b) of the ESI.† Additional computational details are presented in the ESI file.† We present the potential of mean force (PMF) scans obtained from Drude polarizable and additive simulations for the monovalent cations in Fig. 1 and for divalent cations and Cl⁻ ions in Fig. 2. We also present the binding free energies obtained from Drude polarizable and additive FF simulations in Table 1. The binding free energies were estimated with respect to the minimum in the PMF scan and the well separated structure. For the monovalent ions, Li⁺, Na⁺, K⁺ and Cs⁺, the binding free energies from Drude polarizable (additive) FF simulations are estimated to be −6.91 (−7.84), −5.83 (−1.10), −7.13 (−0.56) and −13.26 (−5.38). All energies are reported in kJ mol⁻¹. In Table 1 we also report the difference (ΔE₁) between the Drude polarizable and the additive binding energies. ΔE₁ for Li⁺, Na⁺, K⁺ and Cs⁺ is observed to be 0.93 kJ mol⁻¹, −4.73 kJ mol⁻¹, −6.57 kJ mol⁻¹ and −7.88 kJ mol⁻¹ respectively. We observe that except for Li⁺ the binding free energies estimated from Drude polarizable simulations are higher than the corresponding values obtained from the additive simulations. This indicates a stronger binding between the graphene surface and the monovalent ions in Drude polarizable simulations when compared to the additive simulations. Earlier efforts to study ion–graphene interactions using MD simulations relied on capturing DFT derived binding profiles.^23,24 These have been found to overestimate the binding energetics.³⁸ We present a comparison of the binding free energies obtained by us with those obtained in these earlier studies to put the values obtained by us into perspective. In Table 1 we present the adsorption energies estimated by Carbone et al. using CPCM DFT calculations (E_DFT) and from PMF scans obtained after tuning the additive FF parameters (E_additive^#) to reproduce the CPCM DFT calculations.²³ The E_additive^# values for Li⁺, Na⁺ and K⁺ were reported to be −10.7 kJ mol⁻¹, −14.5 kJ mol⁻¹ and −12.3 kJ mol⁻¹. The values obtained by us using the Drude polarizable FF are lower than those estimated by Carbone et al. In a follow up study Carbone et al. highlighted two potential shortcomings of the Lennard–Jones parameters optimized by them: (i) the model was parametrized using a single ion thereby the effect of ionic screening due to multiple ions was only included via the standard Lorentz–Berthelot combination rules and (ii) the model did not account for the polarization of the graphene surface due to the specific arrangement of the water molecules at the interface.³⁸ We note that these effects might have contributed to the very high adsorption energies observed by Carbone et al. DFT typically overestimates the adsorption of the ions onto the graphene sheet as it does not account for the screening of the charges by the water molecules. Fang et al. studied the interaction of hydrated cations, Li⁺-(H₂O)_n, Na⁺-(H₂O)_n and K⁺-(H₂O)_n with the graphene sheet.¹⁹ The values for Li⁺-(H₂O), Na⁺-(H₂O) and K⁺-(H₂O) were observed to be around −42.0, −35.0 and −29.0 kcal mol⁻¹ respectively, while the same for Li⁺-(H₂O)₉, Na⁺-(H₂O)₉ and K⁺-(H₂O)₉ were observed to be −21.12, −22.27 and −22.98 kcal mol⁻¹. We clearly observe a screening effect upon the inclusion of water which is inversely correlated to the size of the cation; hence parameterizing the FF using a simple screened model would result in overstabilization.¹⁹ In contrast in the Drude polarizable FF the ions were parametrized to be consistent with aqueous bulk thermodynamic properties, such as hydration free energies, self-diffusion coefficients and the energetics of small ion–water clusters, thereby capturing the screening effects introduced by solvation.³⁹ This is also reflected in the PMF curves obtained from the Drude polarizable simulations wherein with increasing size of the monovalent cation, we observe signatures of the stabilization of the hydrated species.

Fig. 1 Potential of mean force (PMF) obtained from additive and Drude polarizable FF simulations of monovalent cations, Li⁺, Na⁺, K⁺ and Cs⁺, considered in this study. All energies are reported in units of kJ mol⁻¹. All distances are presented in units of Å. Representative structures corresponding to the minima present in the PMF surface obtained from Drude polarizable simulations are presented to highlight the interaction modes. Ions, the underlying graphene sheet and water molecules in the 1st hydration shell are presented in VdW representation. The remaining solvent molecules are presented using a solvent environment to illustrate the solvent environment. Only a small section of the sheet and the solvent are presented for clarity. All figures are prepared using VMD.²⁸

Fig. 2 Potential of mean force (PMF) obtained from additive and Drude polarizable FF simulations of divalent cations, Mg²⁺ and Ca²⁺ and counter-ion Cl⁻, considered in this study. All energies are reported in units of kJ mol⁻¹. All distances are presented in units of Å. Representative structures corresponding to the minima present in the PMF surface obtained from Drude polarizable simulations are presented to highlight the interaction modes. Ions, the underlying graphene sheet and water molecules in the 1st hydration shell are presented in VdW representation. The remaining solvent molecules are presented using a solvent environment to illustrate the solvent environment. Only a small section of the sheet and the solvent are presented for clarity. All figures are prepared using VMD.²⁸ The representative structure corresponding to the minima present in the PMF surface obtained from additive simulations of Cl⁻ ions is also presented.

Table 1 Binding free energies of ions obtained from Drude polarizable and additive FF ABF simulations. Binding energies are calculated as the difference between the energies of the equilibrium structure (global minimum) and the well-separated structure. ΔE₁ is calculated using the formula ΔE₁ = E_Drude − E_additive. E_DFT corresponds to the binding free energies obtained by Carbone et al. using CPCM DFT calculations.²³E_additive^# corresponds to the binding free energies obtained by Carbone et al. using modified additive FF parameters.²³ Ionic radii in solution.⁴⁷ All energies are presented in kJ mol⁻¹. Radii are presented in Å

Ion	E Drude (kJ mol^−1)	E additive (kJ mol^−1)	ΔE1	E DFT ^23	E additive ^# ^23 (kJ mol^−1)	Ionic radii (Å)^47
Li^+	−6.91	−7.84	0.93	−10.4	−10.7	0.71
Na^+	−5.83	−1.10	−4.73	−13.8	−14.5	0.97
K^+	−7.13	−0.56	−6.57	−12.6	−12.3	1.41
Cs^+	−13.26	−5.38	−7.88	—	—	1.73
Mg^2+	−1.0	−6.59	5.59	−16.5	−16.3	0.70
Ca^2+	−2.93	−1.03	−1.9	−15.7	−16.3	1.03
Cl^−	−5.26	−5.76	0.5	−6.90	−7.0	1.80

---

From Fig. 1, we note that the binding energy curves for all the monovalent ions are predominantly characterized by two minima, the first one at a distance close to ∼4 Å and the other minima at a distance ≥6 Å. These minima correspond to two distinct interaction modes, (i) wherein the ions interact directly with the graphene surface or as a partially solvated species and (ii) wherein the ions interact with the graphene surface through a solvation shell. The distances corresponding to the global minima and the next minima from both the Drude polarizable and additive simulations are tabulated in Table 2 for all the systems. For Li⁺, we observe that the ions favor interacting directly with the graphene sheet at a distance of 4.20 (4.05) Å in Drude polarizable (additive) FF simulations, with binding free energies of −6.91 (−7.84) kJ mol⁻¹. The solvent separated interaction mode is observed at a distance of 7.40 (7.25) Å with binding free energies of −3.09 (−2.10) kJ mol⁻¹. For the remaining monovalent cations Na⁺, K⁺ and Cs⁺ the solvent separated minima is observed to be the global minima. For Na⁺ ions the global minimum is observed as a solvent separated interaction at a distance of 7.75 Å in both the Drude polarizable and additive simulations. The binding free energy corresponding to this minimum is observed to be −5.83 kJ mol⁻¹ and −1.1 kJ mol⁻¹ form Drude polarizable and additive simulations respectively. We observe that the Drude polarizable simulations also favour the direct interaction of the Na⁺ ions and the graphene sheet with a minimum in the binding free energy profile appearing at a distance of 5.10 Å, with the binding free energy being −3.15 kJ mol⁻¹. This interaction is not favoured in the additive simulations with the binding free energy being 7.13 kJ mol⁻¹. With increasing size of the ionic radii of the ions we observe the stabilization of the solvent separated minima for both the K⁺ and Cs⁺ species. From the Drude polarizable simulations the minima for K⁺ are observed at 5.95 Å, with the binding free energy being −7.13 kJ mol⁻¹. On the other hand, from the additive simulations we observe two shallow minima at 5.85 Å and 8.10 Å, with the corresponding binding energy being −0.53 kJ mol⁻¹. For Cs⁺, a single global minimum is observed at 6.35 Å and 6.25 Å for the Drude polarizable and additive simulations with the binding free energy being −13.26 kJ mol⁻¹ and −5.38 kJ mol⁻¹ respectively. Overall, for the monovalent ions the binding free energy follows the pattern Cs⁺ > K⁺ > Li⁺ > Na⁺ for the Drude polarizable simulations, while for the additive simulations the pattern is Li⁺ > Cs⁺ > Na⁺ > K⁺. It can be inferred that the solvent separated interactions are stabilized in the Drude polarizable FF when compared to the additive FF, which is indicative of differential solvation dynamics observed in the Drude polarizable and the additive simulations.

Table 2 Location of the global minima and next minima observed in the potential energy surface obtained from Drude polarizable and additive PMF simulations. All distances are presented in units of Å

Ion	Drude polarizable FF (Å)		Additive FF (Å)
	Global minima	Next minima	Global minima	Next minima
Li^+	4.20	7.40	4.05	7.25
Na^+	7.75	5.10	7.75	5.10
K^+	5.95	—	8.10	5.85
Cs^+	6.35	—	6.25	—
Mg^2+	7.55	4.80	4.75	7.46
Ca^2+	5.95	5.05	7.75	4.70
Cl^−	3.90	—	6.10	—

---

The PMF scans obtained from Drude polarizable and additive simulations for divalent cations and Cl⁻ ions are presented in Fig. 2. For the divalent ions Mg²⁺ and Ca²⁺, we observe interaction patterns that are consistent with the dependence of the ionic radii of these species. The ionic radii of Mg²⁺ (0.70 Å) are similar to the ionic radii of Li⁺ (0.71 Å). From the additive simulations we observe that similar to Li⁺ the Mg²⁺ ions favour interacting directly with the graphene sheet with the global minima appearing at 4.75 Å. The binding free energy corresponding to this interaction is found to be −6.59 kJ mol⁻¹. The solvent-separated minima are observed at 7.46 Å, with the binding free energy being −4.45 kJ mol⁻¹. From the Drude polarizable simulations we observe two shallow minima at 4.80 Å and 7.55 Å, which correspond to the direct interaction and solvent-separated interaction of Mg²⁺ with the graphene surface. The binding free energy corresponding to these minima is found to be −0.85 kJ mol⁻¹ and −1.00 kJ mol⁻¹ respectively. We observe that for both Li⁺ and Mg²⁺ the additive FF stabilizes the direct interaction of the ions with the graphene surface when compared to the Drude polarizable FF. We notice that this is directly related to the ability of the ions to polarize their surrounding solvation shell. In the Drude polarizable FF simulations the water model captures the influence of polarization while in the additive simulations the rigid water model does not account for the change in the polarizability. This is discussed subsequently in the manuscript. For Ca²⁺ (1.03 Å) the behaviour is similar to the Na⁺ (0.97 Å), which share comparable ionic radii. Similar to Na⁺ the additive FF favours only the solvent-separated interaction between the Ca²⁺ ions and the graphene sheet with the global minima being observed at 7.75 Å and the corresponding binding free energy being −1.03 kJ mol⁻¹. The direct interaction of Ca²⁺ ions with the graphene sheet is not favoured with the free energy cost of the interaction being 2.21 kJ mol⁻¹ for the minima appearing at 4.70 Å. From the Drude polarizable simulations we observe a free energy distribution that is different from that of the additive simulations. Two closely spaced minima are observed at 5.05 Å and 5.95 Å, with the binding free energies being −2.21 kJ mol⁻¹ and −2.93 kJ mol⁻¹. This corresponds to the stabilization of the partially solvated cation interacting with the graphene surface.

The major differences between the additive and Drude polarizable FF are observed for the chloride anions (Cl⁻). From the additive simulations we observe a minimum at a distance of 6.1 Å from the graphene surface, with the binding free energy being −5.76 kJ mol⁻¹. On the other hand, from the Drude polarizable simulations we observe the minimum at a distance of 3.9 Å from the graphene surface, with the binding free energy being −5.26 kJ mol⁻¹. In the Drude polarizable simulations we observe the preference for a direct interaction between Cl⁻ and the graphene surface, while in additive simulations the interaction between the Cl⁻ ions and the graphene surface is mediated via a solvation shell. Both, experimental studies⁴⁰ and ab initio calculations^41–44 have shown that anions interact directly with the graphene surface. Experimental studies using deep UV second harmonic generation measurements found direct ion–graphene interactions to be responsible for the adsorption of SCN⁻ ions on the graphene surface, with the interactions being enthalpically driven.¹⁶ The free energy of adsorption of the thiocyanate to the graphene was estimated to be −8.5 ± 1.1 kJ mol⁻¹. The experimental studies also point towards a direct interaction between the anions and the graphene surface which is captured only by the Drude polarizable FF simulations. It was noted that the anions strongly polarize the graphene surface and the water molecules, which results in the direct interaction between the ions and the graphene surface. This behaviour was not captured by additive simulations and required the inclusion of explicit polarization in molecular dynamics simulations to capture the effect. The anion–graphene interactions have been described earlier by accounting for the polarization using a QM/MD coupling method,³⁸ force field based on a neural network model⁴⁵ and explicitly polarizable force fields.⁴⁶ Our results agree with those of earlier studies on the explicit inclusion of polarization.

In order to capture the influence of the FF on interfacial dynamics in a realistic system we study the dynamics of 1 M salt solutions of LiCl, NaCl, KCl, CsCl, MgCl₂ and CaCl₂. A representative image depicting the system setup is presented in Fig. S1(c) of the ESI.† We compute the density profiles of the ions as a function of the distance from the topmost graphene surface to establish the effect of polarization on the ion–graphene interactions. The density profiles obtained from additive and Drude polarizable FF simulations are presented in Fig. 3. Before commenting on the distribution of the cations we analyse the distribution of the water and the counter-ion Cl⁻ in these simulations. The distribution of the water molecules as a function of distance from the topmost graphene sheet for all the systems is presented in Fig. S2 of the ESI file.† For all the systems, we observe a bimodal distribution with peaks at 3.30 Å and 6.10 Å for the TIP3P additive water model and at 3.40 Å and 6.20 Å for the SWM4 Drude polarizable water model. The first peak corresponds to the direct interaction of the water molecules with the graphene sheet. In Fig. 4 we present the distribution of the counterion Cl⁻ as a function of the distance from the topmost graphene surface for all the systems. In the additive simulations we observe a significant peak at around 6.40 Å for all the systems. This peak is indicative of fully solvated Cl⁻ ions interacting with the graphene surface in the additive simulations. For the Drude polarizable simulations we observe a significant peak at around 3.80 Å which implies a direct interaction between the Cl⁻ ions and the graphene surface. These observations are directly related to the underlying free energy distribution, wherein the free energy minimum is observed at around 6.10 Å for the additive FF and around 3.90 Å for the Drude polarizable FF (Fig. 2). For the Drude polarizable simulations we also observe a dependence on the cations from the peak height of the distribution at 3.80 Å. The height of the peak follows the trend Ca²⁺ > Mg²⁺ > Li⁺ ∼ Na⁺ > K⁺ > Cs⁺. For Ca²⁺ we also observe a strong bimodal distribution of the Cl⁻ ions with a second significant peak at 6.80 Å.

Fig. 3 Density profiles for all ions calculated from additive (black) and Drude polarizable (red) FF simulations of 1 M salt solutions. In each panel we also present the location of the peak corresponding to the Cl⁻ ion density using dashed lines. All densities are presented in units of Å⁻³ and all distances are presented in units of Å.

Fig. 4 Number density (η) for Cl⁻ as a function of distance from the topmost graphene sheet, from 1 M additive and Drude polarizable FF salt simulations. All distances are presented in units of Å.

The ion density distribution for the cations depends on both the size of the cation as well as the underlying FF. For all the systems we observe that the additive simulations favour the interactions between the solvated cation and the graphene sheet, with the major peak in the distribution being observed in the range between 5.7–7.7 Å. We observe that this is being driven by the distribution of the counter ion Cl⁻ in the system which is observed at around ∼6.40 Å for all the systems. For Li⁺, Na⁺, Mg²⁺ and Ca²⁺ we also observe minor peaks in the distribution at 4.9 Å, 4.8 Å, 4.8 Å and 4.7 Å, which correspond to the direct interaction of the ions with the graphene surface or the interaction of a partially solvated ion with the graphene surface. For Li⁺, Na⁺, Mg²⁺ and Ca²⁺ this interaction is favoured due to the favourable free energy associated with such interactions (Fig. 1) in the additive FF. For the Drude polarizable FF we observe two distinct favourable interactions depending on the cations. Li⁺, Mg²⁺ and Ca²⁺ favour a direct interaction between the cations or partially solvated cations and the graphene surface with the major peak in the distribution being observed at 4.0 Å, 4.8 Å and 5.1 Å, respectively. On the other hand, Na⁺, K⁺ and Cs⁺ favour interaction between a fully solvated cation and the graphene surface, with the major peaks being observed at 7.7 Å, 6.1 Å and 6.4 Å respectively. Interestingly for Ca²⁺ we observe the next peak in the distribution at 6.1 Å. This close pacing of the peaks is correlated to the broad shallow minimum observed in the free energy profile for Ca²⁺ obtained using the Drude polarizable FF parameters (Fig. 2). The increased propensity of Li⁺, Mg²⁺ and Ca²⁺ to reside close to the graphene surface in Drude polarizable FF simulations is also driven by the prominent density of Cl⁻ near the graphene surface with the Cl⁻ distribution being observed at 3.8 Å. A decrease in the Cl⁻ ion density at 3.8 Å, which is correlated to the peak height at 3.8 Å (Ca²⁺ > Mg²⁺ > Li⁺ ∼ Na⁺ > K⁺ > Cs⁺) results in a shift in the Na⁺, K⁺ and Cs⁺ densities away from the graphene surface.

Before concluding we analyse the hydration dynamics of the water molecules around the ions from the Drude polarizable FF simulations. In Fig. 5a we present the average hydration number computed for ions within blocks of 1 Å as a function of distance from the graphene sheet. For the monovalent ions we observe that the hydration number increases as a function of the size of the ion. The average hydration number for Li⁺ and Na⁺ was found to be 3.96 and 5.62 irrespective of the distance from the graphene sheet indicative of a near tetra-coordinated and near hexa-coordinated structure for the first hydration shell in these systems. This is also observed in the representative structures presented in Fig. 1. The larger monovalent ions K⁺ and Cs⁺ exhibit variable hydration numbers depending upon the distance from the graphene sheet. The K⁺ ions exhibit a distorted pentagonal bipyramidal geometry in the bulk with the hydration number being 6.92 for the first hydration shell. Closer to the graphene sheet the hydration umber is found to be 5.25 indicating a loss of one/two water molecules from the bulk hepta-coordinated structure. For the largest ion Cs⁺ we observe variable hydration numbers in the bulk, with the hydration number varying from 9.21 to 11.05. This indicates an unstructured dynamic first hydration shell around Cs⁺. Close to the graphene surface the hydration shell easily loses water molecules and the hydration number drops to 8.21. For the divalent ion Mg²⁺ the average hydration number was found to be 6.00 irrespective of the distance from the graphene sheet. However, a closer inspection of the structure of the hydration shell reveals that near the graphene sheet the hydration shell rearranges into a pentagonal pyramidal structure while in the bulk the hydration shell resembles an octahedral structure (Fig. 2). The larger divalent ion Ca²⁺ exhibits a variable hydration number like K⁺ and Cs⁺. Close to the graphene sheet the hydration number is found to be 3.18, while in the bulk solution the hydration number was found to vary between 3.00 and 4.65. This points to a significant deviation from the octa-coordinated crystal environment for the Ca²⁺ ions in CaCl₂ salt. The hydration shell around Cl⁻ also undergoes partial reorganization upon interacting with the graphene surface. The hydration number drops to 6.02 from a bulk value of 6.65. Visualizing the hydration shell structure reveals a pentagonal pyramidal arrangement of water molecules around Cl⁻ similar to the arrangement observed for Mg²⁺ (Fig. 2).

Fig. 5 (a) First hydration number (η₁) for Li⁺, Na⁺, K⁺, Cs⁺, Mg²⁺, Ca²⁺ and Cl⁻ ions from 1 M ion–graphene simulations as a function of distance from the graphene surface. (b) Average dipole moment (μ) for water molecules in various 1 M salt solutions discussed in the study. All distances are presented in units of Å. Dipole moments are presented in units of Debye.

Finally, we analyse the dipole moment of the water molecules as a function of the distance from the graphene surface. In Fig. 5b we present the average dipole moment computed for water molecules within blocks of 0.5 Å as a function of distance from the graphene sheet. We note that, from additive FF simulations, the dipole moment of water molecules is constrained to be 2.374 D due to the rigid water model used in the simulations. However, the SWM4 water model used in the Drude polarizable FF simulations is able to respond to the changes in the local environment. Very close to the graphene sheet the average dipole moment of water molecules stabilizes at around 2.41 D, wherein the water molecules are in close contact with the hydrophobic graphene and the influence of ions is not felt. Moving away from the graphene sheet we begin to observe an increase in the ion densities of both the anions and the cations (Fig. 3). For MgCl₂ and LiCl simulations, we observe a significant change in the average dipole moment of water with the dipole moment increasing to 2.58 D and 2.51 D at a distance of 5.7 Å from the graphene sheet. This is due to the small size of Mg²⁺ (0.70 Å) and Li⁺ (0.71 Å) ions, which favours a large charge/surface area ratio resulting in significant polarization of the local environment. On the other hand, for the larger cations we observe a reduced impact of the ions on the average dipole moment of water.

To comment on ion–graphene interactions and how the charges on the ions and the graphene surface respond to such interactions we analyse the dipole moment fluctuations for both the ions and the interacting graphene surface. In Fig. S3 of the ESI† we present the instantaneous dipole moment of the ions as a function of distance from the graphene surface. In Fig. S4 of the ESI† file we present the instantaneous dipole moment of the interacting graphene surface as a function of distance from the ions. The instantaneous dipole moments of the ions (Fig. S3†) and interacting graphene surface (Fig. S4†) are found to be invariant with respect to the distance between the ions and the graphene surface. This indicates an absence of direct interaction between the ions and the graphene surface. This is highlighted in the observed distribution of the instantaneous dipole moment of the graphene sheet (Fig. S4†), which appears independent of the chemical nature of the ions. However for both the ions and the graphene sheet we observe a spread in the dipole distribution. For the ions (Fig. S3†) we observe a strong correlation between the distribution of the instantaneous dipole moments and the charge to surface ratio of the ions, with small ions such as Li⁺ and Mg²⁺ exhibiting a very small dipole distribution and larger ions such as Cs⁺ and Ca²⁺ exhibiting a broader distribution. This distribution can be traced back to the ion–water interactions in the first solvation shell of the ions. For the graphene sheet we observe a distribution from 0.0 D to 2.4 D (Fig. S4†). This dipole response of the graphene surface can be traced back to the interactions with the water molecules present in the 1st monolayer of water observed at 3.40 Å (Fig. S2†). Thus, we only observe solvent driven interactions between the ions and the graphene surface.

D Conclusions

In conclusion, we show that the Drude polarizable FF parameters for graphene²⁶ along with the parameters for ions³⁹ reliably capture the dynamics at the graphene electrolyte interface. The Drude parameters are able to capture both the size and charge dependent ion–graphene interactions, which cannot be captured by additive simulations. In particular the Drude parameters are able to capture the anion–graphene interactions, which are severely underestimated in the additive simulations. Additionally, the Drude polarizable simulations also capture the influence of ions on solvation dynamics. These results establish the applicability of Drude parameters for studying ion–graphene interface interactions, like those observed in graphene-based membranes for desalination,^1–5 or graphene-based electrodes in energy storage.^6–9 The Drude parameters also present an effective alternative to computationally expensive first-principles MD simulations.

Author contributions

Hemanth H. (HH) and Rohan Mewada (RM) have contributed equally to the work. The work was designed by Sairam S. Mallajosyula (SSM). HH and SSM analyzed the results and wrote the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

H. H. thanks the Indian Institute of Technology Gandhinagar for the SRF fellowship. H. H. and R. M. thank the Indian Institute of Technology Gandhinagar for infrastructure and facilities. S. S. M. thanks the Department of Science and Technology, India (DST) and National Supercomputing Mission (NSM) for the research grant (DST/NSM/R&D_HPC_Applications/Sanction/2021/20). S. S. M. thanks the Indian Institute of Technology Gandhinagar for High Performance Computing Facility. S. S. M. thanks the Japan Advanced Institute of Science and Technology for High Performance Computing Facility. S. S. M. thanks the Indian Institute of Technology Banaras Hindu University for providing accesses to PARAM SHIVAY computing facility.

Notes and references

1. P. Sun, F. Zheng, M. Zhu, Z. Song, K. Wang, M. Zhong, D. Wu, R. B. Little, Z. Xu and H. Zhu, ACS Nano, 2014, 8, 850–859.
2. Q. Yang, Y. Su, C. Chi, C. T. Cherian, K. Huang, V. G. Kravets, F. C. Wang, J. C. Zhang, A. Pratt, A. N. Grigorenko, F. Guinea, A. K. Geim and R. R. Nair, Nat. Mater., 2017, 16, 1198–1202.
3. R. K. Joshi, P. Carbone, F. C. Wang, V. G. Kravets, Y. Su, I. v. Grigorieva, H. A. Wu, A. K. Geim and R. R. Nair, Science, 2014, 343, 752–754.
4. S. Homaeigohar and M. Elbahri, NPG Asia Mater., 2017, 9, e427.
5. A. Boretti, S. Al-Zubaidy, M. Vaclavikova, M. Al-Abri, S. Castelletto and S. Mikhalovsky, NPJ Clean Water, 2018, 1, 5.
6. M. Liang and L. Zhi, J. Mater. Chem., 2009, 19, 5871–5878.
7. G. Wang, B. Wang, X. Wang, J. Park, S. Dou, H. Ahn and K. Kim, J. Mater. Chem., 2009, 19, 8378–8384.
8. Y. Wang, Z. Shi, Y. Huang, Y. Ma, C. Wang, M. Chen and Y. Chen, J. Phys. Chem. C, 2009, 113, 13103–13107.
9. M. D. Stoller, S. Park, Z. Yanwu, J. An and R. S. Ruoff, Nano Lett., 2008, 8, 3498–3502.
10. Z. H. Sheng, X. Q. Zheng, J. Y. Xu, W. J. Bao, F. bin Wang and X. H. Xia, Biosens. Bioelectron., 2012, 34, 125–131.
11. P. Martin, Chem. Rec., 2009, 9, 211–223.
12. M. Pumera, A. Ambrosi, A. Bonanni, E. L. K. Chng and H. L. Poh, Trac. Trends Anal. Chem., 2010, 29, 954–965.
13. R. R. Nair, H. A. Wu, P. N. Jayaram, I. v. Grigorieva and A. K. Geim, Science, 2012, 335, 442–444.
14. J. Abraham, K. S. Vasu, C. D. Williams, K. Gopinadhan, Y. Su, C. T. Cherian, J. Dix, E. Prestat, S. J. Haigh, I. v. Grigorieva, P. Carbone, A. K. Geim and R. R. Nair, Nat. Nanotechnol., 2017, 12, 546–550.
15. B. Mennucci, J. Tomasi, R. Cammi, J. R. Cheeseman, M. J. Frisch, F. J. Devlin, S. Gabriel and P. J. Stephens, J. Phys. Chem. A, 2002, 106, 6102–6113.
16. D. L. McCaffrey, S. C. Nguyen, S. J. Cox, H. Weller, A. P. Alivisatos, P. L. Geissler and R. J. Saykally, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 13369–13373.
17. P. Iamprasertkun, W. Hirunpinyopas, A. Keerthi, B. Wang, B. Radha, M. A. Bissett and R. A. W. Dryfe, J. Phys. Chem. Lett., 2019, 10, 617–623.
18. E. D. Glendening and D. Feller, J. Phys. Chem., 1996, 100, 4790–4797.
19. L. Mu, Y. Yang, J. Liu, W. Du, J. Chen, G. Shi and H. Fang, Phys. Chem. Chem. Phys., 2021, 23, 14662–14670.
20. H. Li, J. S. Francisco and X. C. Zeng, Proc. Natl. Acad. Sci. U.S.A., 2015, 112, 10851–10856.
21. J. Sunner, K. Nishizawa and P. Kebarle, J. Phys. Chem., 1981, 85, 1814–1820.
22. K. Zhou and Z. Xu, Phys. Rev. Res., 2020, 2, 042034.
23. C. D. Williams, J. Dix, A. Troisi and P. Carbone, J. Phys. Chem. Lett., 2017, 8, 703–708.
24. L. Chen, Y. Guo, Z. Xu and X. Yang, ChemPhysChem, 2018, 19, 2954–2960.
25. F. Y. Lin, P. E. M. Lopes, E. Harder, B. Roux and A. D. Mackerell, J. Chem. Inf. Model., 2018, 58, 993–1004.
26. H. Hemanth and S. S. Mallajosyula, Nanoscale, 2021, 13, 4060–4072.
27. H. Hemanth, P. K. Yadav and S. S. Mallajosyula, J. Phys. Chem. C, 2022, 126, 13122–13131.
28. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graph., 1996, 14, 33–38.
29. K. Vanommeslaeghe, E. Hatcher, C. Acharya, S. Kundu, S. Zhong, J. Shim, E. Darian, O. Guvench, P. Lopes, I. Vorobyov and A. D. Mackerell, J. Comput. Chem., 2010, 31, 671–690.
30. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1998, 79, 926.
31. G. Lamoureux, E. Harder, I. v. Vorobyov, B. Roux and A. D. MacKerell, Chem. Phys. Lett., 2006, 418, 245–249.
32. S. Miyamoto and P. A. Kollman, J. Comput. Chem., 1992, 13, 952–962.
33. J. C. Phillips, D. J. Hardy, J. D. C. Maia, J. E. Stone, J. v. Ribeiro, R. C. Bernardi, R. Buch, G. Fiorin, J. Hénin, W. Jiang, R. McGreevy, M. C. R. Melo, B. K. Radak, R. D. Skeel, A. Singharoy, Y. Wang, B. Roux, A. Aksimentiev, Z. Luthey-Schulten, L. v. Kalé, K. Schulten, C. Chipot and E. Tajkhorshid, J. Chem. Phys., 2020, 153, 044130.
34. J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kalé and K. Schulten, J. Comput. Chem., 2005, 26, 1781–1802.
35. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1998, 98, 10089.
36. J. Comer, R. Chen, H. Poblete, A. Vergara-Jaque and J. E. Riviere, ACS Nano, 2015, 9, 11761–11774.
37. H. Poblete, I. Miranda-Carvajal and J. Comer, J. Phys. Chem. B, 2017, 121, 3895–3907.
38. J. D. Elliott, A. Troisi and P. Carbone, J. Chem. Theory Comput., 2020, 16, 5253–5263.
39. H. Yu, T. W. Whitfield, E. Harder, G. Lamoureux, I. Vorobyov, V. M. Anisimov, A. D. MacKerell and B. Roux, J. Chem. Theory Comput., 2010, 6, 774–786.
40. H. Lee, J. il Choi, J. Park, S. S. Jang and S. W. Lee, Carbon, 2020, 167, 816–825.
41. G. Shi, Y. Ding and H. Fang, J. Comput. Chem., 2012, 33, 1328–1337.
42. D. Kim, P. Tarakeshwar and K. S. Kim, J. Phys. Chem. A, 2004, 108, 1250.
43. X. Liu and G. Shi, Chin. Phys. B, 2021, 30, 046801.
44. F. Xiaozhen, L. Xing, H. Zhenglin, Z. Kaiyuan and S. Guosheng, J. Mol. Model., 2022, 28, 1–9.
45. N. di Pasquale, J. D. Elliott, P. Hadjidoukas and P. Carbone, J. Chem. Theory Comput., 2021, 17, 4477–4485.
46. R. P. Misra and D. Blankschtein, J. Phys. Chem. C, 2021, 125, 2666–2679.
47. Y. Marcus, Chem. Rev., 1988, 88, 1475–1498.

---

Footnote

† Electronic supplementary information (ESI) available: Representative structures depicting various systems investigated. Number density (η) plots from additive and Drude polarizable FF 1 M simulations for water molecules as a function of distance from the graphene surface. Analysis of instantaneous dipole moments for ions and graphene surfaces. See DOI: https://doi.org/10.1039/d2na00733a

---