Properties of aqueous electrolyte solutions at carbon electrodes: effects of concentration and surface charge on solution structure, ion clustering and thermodynamics in the electric double layer

Abstract

Surfaces are able to control physical–chemical processes in multi-component solution systems and, as such, find application in a wide range of technological devices. Understanding the structure, dynamics and thermodynamics of non-ideal solutions at surfaces, however, is particularly challenging. Here, we use Constant Chemical Potential Molecular Dynamics (CμMD) simulations to gain insight into aqueous NaCl solutions in contact with graphite surfaces at high concentrations and under the effect of applied surface charges: conditions where mean-field theories describing interfaces cannot (typically) be reliably applied. We discover an asymmetric effect of surface charge on the electric double layer structure and resulting thermodynamic properties, which can be explained by considering the affinity of the surface for cations and anions and the cooperative adsorption of ions that occurs at higher concentrations. We characterise how the sign of the surface charge affects ion densities and water structure in the double layer and how the capacitance of the interface—a function of the electric potential drop across the double layer—is largely insensitive to the bulk solution concentration. Notably, we find that negatively charged graphite surfaces induce an increase in the size and concentration of extended liquid-like ion clusters confined to the double layer. Finally, we discuss how concentration and surface charge affect the activity coefficients of ions and water at the interface, demonstrating how electric fields in this region should be explicitly considered when characterising the thermodynamics of both solute and solvent at the solid/liquid interface.

---

Introduction

Carbon–electrolyte interfaces often feature in technologies and devices designed for energy storage^1–3 and water desalination.^4,5 Moreover, carbon allotropes are increasingly employed as nano-reactors,⁶ as well as supports for liquid-phase catalysts.⁷ A molecular-level picture of the structure and dynamics of multi-component liquid phases at the carbon interface is important to understand the physical chemistry involved in such technologies/devices in order to improve their design for functional applications. Molecular simulations, particularly molecular dynamics (MD), provide powerful tools to investigate such systems at the atomic level.⁸ By explicitly capturing the atomistic details of the solid/liquid interface, MD-based methods enable predictions regarding the effect of changes to the bulk solution composition and the applied interfacial potential (that gives rise to a surface charge) on the properties of the so-called electric double layer (EDL). In turn, this allows for an assessment of the suitability of mean-field models that are commonly used to describe and predict the structure and electrochemical properties of solid–solution interfaces.⁹

Gouy–Chapman theory predicts a monotonically decreasing concentration of ions in the immediate vicinity of electrodes with the same sign of charge, while the concentration of ions with opposite charge to the surface increases smoothly according to a Boltzmann distribution; thus, the solution screens the surface charge by establishing a diffuse EDL.¹⁰ This fundamental model for the structure of charge carriers at electrodes inadequately describes the EDL when large potentials are applied and in the presence of high electrolyte concentrations. By neglecting ion finite-sizes and their correlations, it fails to explain the change in the electrical properties of the graphite–electrolyte interface due to specific ion effects, which was demonstrated across the series of alkali chlorides at graphite.^11,12 The simple picture of the EDL was later developed to address some of these shortcomings, by accounting for the specific adsorption of ions at the electrode and the role that ion solvation spheres play in defining the inner- and outer-Helmholtz plane.¹³ In this framework, the solution-side of the EDL is modelled as a series of plate capacitors; nonetheless, it is assumed that the finite size of charge carriers can be ignored in the diffuse layer. At low concentrations of simple salts—such as NaCl—in water, these simple mean-field-based models were suggested to provide a reasonable approximation of the EDL structure,¹⁴ especially as charge transfer between the electrode and charge carriers in solution is low.¹² However, the combination of high salt concentrations and large surface charge densities results in conditions where the solvation, finite size and cooperative adsorption of ions cannot be neglected. More sophisticated mean-field models of the EDL were developed to address some of these effects.^15–21

Our recent simulations demonstrated how asymmetric electrolyte adsorption gives rise to alternating cation and anion-rich aqueous solution layers perpendicular to planar graphene and graphite substrates at moderate-to-high alkali chloride solution concentrations (∼1 M and above).^22,23 This behaviour is due to the partial saturation of ions in solution layers in contact with the surface that emerges in the EDL.^8,23–26 This picture of the EDL is reminiscent of the structures observed in ionic liquids at charged surfaces and requires a treatment of the EDL that accounts for the finite size of charge carriers accumulating at the interface.^14,27,28 The asymmetric ordering of ions results in charge fluctuations in this region—typically four-to-five liquid layers deep—and a departure from descriptions of the EDL expected from the established mean-field models described above.^9,29

Thanks to the adoption of the Constant Chemical Potential Molecular Dynamics (CμMD) method,^30,31 which maintains a constant thermodynamic driving force associated with ion adsorption, we were able to quantify the electric potential drop across the EDL and the excess chemical potential for ions at the solid–solution interface,^22,32 In CμMD, the use of an explicit molecular reservoir coupled to the model interface prevents any ion depletion in the bulk solution, which would otherwise occur in typical finite-sized MD simulations when ions adsorb at an interface.

Here, we extend our analysis to consider concentrated NaCl(aq) solutions in contact with charged graphite and the resulting properties of the solution side of the EDL. In our analysis of the simulation results, we pay particular attention to the thermodynamic and structural properties of the solvent (as well as ion speciation). Understanding how the presence of ions and surface charge control the thermodynamics of solvent is essential to predict the activity of interfaces for applications in catalysis, and recent computational studies have demonstrated how interfaces impact the ability of the aqueous medium to screen coulombic interactions due to a changing dielectric constant.³³

In what follows, we recap the effect of concentration on the structure and properties of ions in solution at neutral graphite before considering the combined effects of concentration and surface charge. We characterise the structural properties of water molecules in the EDL when compared to bulk solutions, pure liquid water and ice. Finally, we evaluate the electrical properties of the EDL and use this information to calculate how the activity constants for ions and water change on moving from the bulk solution towards the graphite surface.

Computational methods

Following the protocol proposed by Finney et al.,²² all simulations were performed using the Joung and Cheatham³⁴ force field to describe the interactions of ions with SPC/E water.³⁵ Graphite was modelled using the OPLS/AA force field,³⁶ while the intermolecular interactions between carbon and water were modelled using pairwise potentials fitted to water adsorption energies obtained via random phase approximation calculations.^37,38 Several force fields are available to model the interactions of carbon with water, and comparisons of some of the different models are available in the literature.^38–41

The carbon–water model adopted here predicts a water contact angle of ∼40°, with small changes to this mean value being dependent upon the number of carbon layers in the substrate and the truncation distance used in the interaction potential.³⁸ The contact angle is more acute than that predicted by earlier force fields; however, it was shown in experiments that graphene becomes less hydrophilic when exposed to air and the surface becomes populated by contaminants (hence, a smaller contact angle should be reproduced by the model than was initially thought).^42–44 The contact angles for pristine graphene and graphite were found to be 42 ± 7° and 42 ± 3°, respectively.^43,44 A recent exhaustive computational study of the interaction energies of water with graphene using quantum mechanical calculations suggest an upper bound to the contact angle of water on graphene of 56°, as informed by dynamical simulations of coarse-grained water molecules at the carbon surface, where interaction potentials were fitted to the results from calculations at a higher level of theory.⁴⁵ Our model, therefore, captures reasonably well the thermodynamics of water at the carbon interface (as determined by the surface tension); furthermore, it predicts the correct radial breathing mode frequency for carbon nanotubes in water.³⁸

Ion–carbon interactions were modelled using potentials fitted to the results from electronic structure calculations that capture the polarisability of the carbon surface in the presence of ions surrounded by a conductor-like polarisable continuum, mimicking the presence of a solvent.⁴⁶ Despite components of the force field being constructed from various sources, it is important to recognise that consistent descriptions of ions and water molecules (i.e., Joung and Cheatham ions and SPC/E water) were used for the fitting of pairwise potentials throughout.

The GROMACS 2018.6 (ref. 47) MD engine was adopted to perform simulations within the NVT ensemble unless otherwise stated. Atom positions were evolved during the simulations using a leapfrog time integrator with a 2 fs timestep; as such, water intramolecular degrees of freedom were constrained using the LINCS algorithm.⁴⁸ Intermolecular interactions were computed for atoms within 0.9 nm, and long-range electrostatics were treated using smooth particle mesh Ewald summation.⁴⁹ The temperature was held constant at 298 K (within fluctuations) using the Bussi–Donadio–Parrinello thermostat.⁵⁰

The simulation set-up for graphite in contact with NaCl(aq) solutions follows our previous work.²² This involved preparing an eight-carbon layer 2.7 × 5.4 × 5.5 nm (x × y × z) graphite slab with basal surfaces perpendicular to the simulation cell x-axis. 1672 Na⁺ and Cl⁻ ions, as well as 13 819 water molecules, were placed in the orthorhombic simulation cell with periodic boundaries in all three dimensions. With carbon atoms fixed at their lattice positions, a 0.2 ns MD simulation was performed to relax the simulation cell volume in the NPT ensemble using the barostat of Berendsen et al.⁵¹ at a pressure of 1 atm. Following this equilibration step, the ions in solution were accumulated in a reservoir region far away from the graphite surface(s) by applying an external harmonic potential to the distance between the surface and ions (using the PLUMED v2.5 plugin⁵² with force constant, k = 3 × 10⁵ kJ mol⁻¹). The minimum distance between carbon atoms and ions in this external bias was 6 nm. Simulations in the NVT ensemble were performed until the ions were at least 5.9 nm from the carbon slab. The final configuration of the system from this preparatory phase was taken as the starting structure for 100 ns CμMD simulations, where the final 50 ns steady-state trajectory window was used in all analyses of the interfacial properties. A similar procedure was used to prepare simulations of graphite in contact with pure water; here, however, no ions were included in the simulation cell, and it was not necessary to prepare the ionic reservoir.

For all CμMD simulations, PLUMED 2 (ref. 52) was utilised to compute the external forces required to control the solute density (n) in a 2.2 nm control region, whose innermost edge (closest to a graphite basal plane) was x_F = 3.7 nm from the centre of the simulation cell x-axis, defining the origin. The CμMD force on ion i takes the functional form,

where ω tunes the width of the force region and was taken to be 0.01% of the cell length in x; superscript CR and t indicate the instantaneous and target value of n in the control region; and k_i = 2 × 10⁵ kJ mol⁻¹ is the force constant for the function that acts like a semi-permeable membrane for the ions.

Standard MD simulations were performed to simulate graphite in contact with water, which we subsequently refer to as 0 mol dm⁻³ (M). In the presence of ions, on the other hand, CμMD simulations were performed where the target ion density was 0.6022, 3.0110 and 6.022 nm⁻³, equating to molar concentrations of 1, 5 and 10 M. When the simulations reached a steady state, the concentrations of ions were maintained at 1.2 ± 0.03, 5.01 ± 0.05 and 9.23 ± 0.07 M. The difference between the target and evaluated concentrations is small and is due to the relative occupancy of the ionic reservoir and the parameters used to apply CμMD forces; this is not particularly important for the current study, and the simulations can be prepared to ensure that concentrations precisely match the target value, if necessary. Both 5 and 10 M cases are beyond the solubility for NaCl(s), determined to be approximately 3.5 M for the Joung and Cheatham force field; the solubility of halite is 3.7 mol kg⁻¹,⁵³ which is approximately 3.5 M according to a 2^nd order polynomial fitting of molarities (c) as a function of molalities (b) obtained from steady-state bulk NaCl(aq) solutions ranging from 1 to 16 mol kg⁻¹ (b = −0.0174c² + 0.9822c + 0.0537, where b and c here refer to numeric values ignoring units).

We explored the effects of applied surface charge in systems at all three sampled bulk solution concentrations as well as in simulations of graphite in contact with pure water. To achieve this, we applied uniform charges to the outermost carbon atoms in the graphite slab; equal charges with the opposite sign were applied to 1144 carbon atoms on each face of the graphite slab to generate charge densities, |σ| = 0.19, 0.39, 0.58 and 0.77e nm⁻². (We use the descriptors positively/negatively charged surface and positive/negative electrode interchangeably throughout.) As such, a single simulation provides information on the effect of positive and negative applied potentials by examining the interface on different sides of the graphite slab. The total charge in the simulation cell was, therefore, zero. We believe that this approach to applying charges to the surface is reasonable for the current study; indeed, when we tested a Drude oscillating charge model for the surface charge polarisability at 1 M, we did not observe significant differences in the properties of the interface when compared with the uniform distribution of charge to carbon centres discussed below. More sophisticated models of the surface charge polarisability have been developed;^24,54 although a comparison of the accuracy and applicability of these methods is beyond the scope of the current work.

Results and discussion

Solution structure at charged graphite

In this section, we discuss the salient features of the steady-state structure of NaCl(aq) solutions at graphite surfaces when negative and positive surface charges are applied to carbon atoms, equating to surface charge densities (σ) in the range 0 to ±0.77e nm⁻². The target bulk ion solution concentrations in our CμMD simulations here were 1, 5 and 10 M. In addition, we also perform simulations of neutral and charged graphite in contact with pure water. A snapshot of a typical CμMD simulation is provided in Fig. 1A, where the ion-rich reservoir can be seen spanning the periodic boundaries, far from the carbon–electrolyte interface. As ions accumulate in the EDL, the solution in contact with the graphite surface is replenished with ions from the reservoir to maintain a constant thermodynamic driving force for the process, ensuring that the solution, several nanometres from the interface—which we refer to as the bulk—is electroneutral. In the following subsections, we focus our analysis on the steady state that emerges between species in solution in this bulk region and the EDL.

Fig. 1 (A) Snapshot from a CμMD simulation of NaCl(aq) at charged graphite projected onto cell x and y axes. The dashed lines indicate the edge of the reservoir (Res.); grey, blue, cyan, red and white spheres indicate C, Na⁺, Cl⁻, O and H atoms, respectively. (B) Mean solution atom densities, ρ, determined as a function of distance from the graphite electrode, Δx, for a range of target bulk solution ion concentrations. The solution concentration and atom type are indicated at the top and left of the grid. Colours blue → red indicate increasing solution charge densities from |σ| = 0 → 0.77e nm⁻².

Density profiles—concentration and charge effects. The one-dimensional atom densities of solution species perpendicular to the graphite basal surface are reported in Fig. 1B. In addition, Fig. S1–S4† provide the same densities on a linear scale over a wider range of Δx. In line with our previous simulation studies^22,32 and those from others using different force fields and graphene,^23–26 the densities indicate a preference for cation adsorption in the first solution layer above the substrate; this is due to the favourable interactions between positively charged ions and the electron-rich carbon surface (implicitly captured by the force field).⁴⁶ At the lowest concentrations, we observe a diffuse anion-rich solution layer adjacent to the first cation-rich layer.

As discussed in detail by Finney et al.,²² the asymmetric adsorption gives rise to a surface potential, the magnitude of which is governed by the sharp cation density in x in the first solution layer, even in the absence of an applied surface charge. This effective surface charge is screened by a diffuse anion layer and, at concentrations below 0.6 M, the concentration profile leading to such charge screening is qualitatively consistent with simple mean-field models of the EDL.²² However, at ∼1 M and above, additional cation and anion density peaks are observed in the EDL, and a complex multi-layered solution structure emerges due to the finite size and cooperative adsorption of ions that is explicitly captured by atomistic simulations and is apparent from the blue curves in Fig. 1B. This picture is reminiscent of the structure of ionic liquids at planar surfaces, where the finite size of charge carriers cannot be ignored.^14,27,28 These results highlight the inadequacy of simple mean-field models to predict the structure and, ultimately, the electrochemical properties of the interface for even simple electrolyte solutions at moderate to high ion concentrations. For a more detailed discussion of the solution structure at uncharged graphite, see the discussion by Finney et al.²²

In this work, we focus on the effect of varying graphite surface charges on water and ion atom densities in the EDL. Such variations are reported more clearly in Fig. S5.† In pure water, increasing the surface charge density to +0.77e nm⁻² results in positive and negative changes to the water oxygen and hydrogen densities (ρ_Ow and ρ_Hw), respectively, in the first solution layer adjacent to the graphite basal surface. Essentially, the surface charge induces an increased ordering of water molecules locally in the vicinity of the positively charged carbon surface. At the negative electrode, however, we observe a restructuring of the liquid in the first two water layers as the magnitude of the surface charge, σ, increases. This can be explained by water molecules reorienting to increase the interactions between H-atoms, bearing a positive partial atomic charge, and the excess negative charge uniformly distributed amongst carbon atoms in the outermost graphite layer. The reorientation of water molecules manifests in the density profiles as a splitting of the first ρ_Hw peak into two peaks separated by ∼0.1 nm. The electrostatic repulsion of water oxygen atoms with the surface also displaces molecules in the first liquid layer, as shown by a decrease in ρ_Ow in the first peak and a shallower minimum in density between the first two water layers. Regardless of the sign or magnitude of σ, perturbations to the liquid structure encompass approximately three water layers, up to ∼1 nm from the graphite basal plane, consistent with other studies of carbon–solution interfaces.^24,33

Interestingly, as ions are added to the system, the amplitude for the fluctuations in water density at the interface somewhat diminish, although perturbations to the water structure cf. the bulk extend further into the solution as the magnitude of the applied surface potential increases. This can be explained by the screening of the surface charge due to the ions accumulating in the EDL. For example, at 10 M, the ordering of water is observed four-to-five water layers from the surface, and the spacing between peak centres in ρ_Hw decrease, most notably when σ is negative, due to the complex ion layering that is found under these conditions at the interface. This can be reconciled by considering how the surface displaces ions with associated water molecules in their solvation spheres—cations, in particular, have a relatively strong solvent coordination sphere, and so any change to the steady-state structure of ρ_Na will affect the density of water molecules in the EDL.

Fig. 1B and S4† show that a single peak in ρ_Hw is observed in the first solution layer at all values of σ at 10 M, denoting an inhibition of the water structuring found in the absence of electrolyte. Moreover, a large cation density in the first solution layer gives rise to a large increase in the local anion density, which is not apparent at 1 M; hence, the cooperative accumulation of ions displaces water in the second solution layer. Furthermore, the screening of the surface charge by ions at high concentrations mitigates any reorientation of water dipoles. Fig. S5† indicates that at low concentrations, the density of water in the first solution layer increases at the negatively charged surface due to an increased cation density; however, at 10 M, the density change is negative, demonstrating the complex restructuring of water that occurs in the EDL at high solution concentrations.

The perturbation of the solution structure under the effect of increasing surface charge is analogous to an ‘accordion-like’ deformation, where solution layers are compressed under the action of the additional coulombic forces. This compression of the solution layers gives rise to further perturbations to the solution structure, as increased ion ordering propagates the effect of the surface charge into solution: a radically different picture of the EDL than those predicted by simple Poisson–Boltzmann-based models. This feature is also observed as a function of concentration, where partial saturation of the solution with ions in the first layers above the surface result in further deviations from the bulk structure moving away from the interface (compare the blue curves for ion density in Fig. 1B).

In general, we conclude that increasing the solution concentration and surface charge have analogous effects on the solution structure in the EDL. While it is possible to capture the effects of asymmetric adsorption and ion correlations in mean-field models of the EDL^15–18 and the role that water structuring can play in screening the surface potential,¹⁹ the complex solution structure observed at the high concentrations and surface charges here suggests that an explicit model for atoms in the EDL is necessary to capture these cooperative, emergent effects.

EDL relaxation in the presence of surface charge. In order to study the collective dynamics of ionic species at the interface, it is useful to consider how the application of a surface charge changes the composition in the EDL as a function of time. As such, we performed five additional simulations at 1, 5 and 10 M, where the initial configurations were taken from 50, 60, 70, 80 and 90 ns time points in simulations where the graphite had no applied charge. In these simulations, however, we set σ to ± 0.77e nm⁻² on opposite surfaces of the graphite slab. By starting from an equilibrated steady-state structure in the absence of surface charge, these simulations allow us to investigate how the EDL evolves in time when a surface charge is instantaneously applied.

Fig. 2A shows the average change in the number of cations and anions within 2.5 nm from the outermost carbon layer of the graphite surface. At all simulated concentrations, there is a rapid change in the number of ions in the EDL (ΔN_ions). Indeed, on the log scale provided, the change in ΔN_ions from zero is not apparent, as this occurs during the first 0.05 ns. The change is, therefore, extremely rapid as ions are displaced to minimise electrostatic forces. This behaviour might be expected for this system which is often adopted as a model system to study electric double-layer capacitors (EDLCs).^1,2 EDLCs offer a high power density, able to deliver and absorb electrical energy at a much higher rate than typical batteries through rapid charge/discharge cycles. The excellent cycling capability of EDLCs—which typically undergo millions of charge–discharge cycles (instigated by changing the applied potential)—is possible without significant degradation of the interface, making them attractive options for energy storage devices.⁵⁵

Fig. 2 (A) Changes to the mean number of ions (ΔN_ions) in the EDL (defined as a 2.5 nm region adjacent to the carbon basal surface) when a surface charge is switched on with σ = ±0.77e nm⁻². Solid and dashed lines pertain to positively and negatively charged surfaces, and the blue and red colours indicate cations and anions, respectively. The data are averages from five independent simulations with the uncertainties calculated as the standard error of the mean in the data, as shown by the shaded regions. (B) A summary of the data in panel A, where the green and cyan points provide the differences between ΔN^Na and ΔN^Cl at the counter and co-electrode in the EDL. Dashed lines provide a linear interpolation between the data as a guide, while the blue dashed line provides the magnitude of the surface charge in units of e.

Despite the rapid change in ΔN_ions, it is clear from Fig. 2A that a longer-timescale relaxation of the EDL composition develops over tens of nanoseconds. At 1 M, the preference for cations to adsorb in the first solution layer means that there is an asymmetry in the displacement of Cl⁻ at the negative electrode when compared with Na⁺ at the positive electrode. This means that after 30 ns, ΔN^Na increases on both electrodes and so does ΔN^Cl, such that ΔN^Na_σ>0 ≈ 0. At all stages throughout the relaxation, the surface charge is screened by ions in the EDL (ΔN^Na_σ<0 − ΔN^Cl_σ<0 ≈ ΔN^Cl_σ>0 − ΔN^Na_σ>0 ≈ |σ|).

A different behaviour is observed as the concentration of ions is increased. At 5 M, Fig. 2A shows that after 30 ns, the surface charge is almost completely screened by incorporation of cations or their removal from the EDL, whereas ΔN^Cl ≈ 0, regardless of the sign of the applied potential. When the concentration is increased to 10 M, ΔN^Na_σ<0 increases and ΔN^Cl_σ<0 is positive; essentially, as more cations are accumulated in the EDL beyond the number necessary to screen the surface charge, ion–ion correlations also induce an increase to the number of anions in the EDL. Irrespective of this change in behaviour, the surface charge remains screened by changes to ion concentrations in the EDL throughout the period of relaxation.

Fig. 2B summarises the changes to the EDL ion concentrations. Here, we present the differences in absolute changes to the ion concentrations at the counter- (where the sign of the ion charge is opposite to the surface charge) and co-electrodes (where the sign of the charge on ions matches that of the surface charge) for Na⁺ and Cl⁻ as |ΔN^Na_σ<0| − |ΔN^Na_σ>0| and |ΔN^Cl_σ>0| − |ΔN^Cl_σ<0|. The data indicate the affinity of the charged surfaces for cations or anions. A value of zero in Fig. 2B indicates equivalent displacement of the ions at the positive and negative electrode, which is the case only when c(NaCl) is 5 M. At 1 M, on the other hand, the surface can accumulate a net excess of ions in equal amounts at both positive and negative σ. Finally, at 10 M, the net accumulation of cations exceeds that of anions when comparing surfaces with opposite charges. Together, these data highlight the multifaceted relaxation of the EDL structure due to asymmetric ion effects and cooperative changes that result, which are typically neglected in analytical models of the EDL.

Ion association in solution

In our previous work, we demonstrated how correlations in bulk solutions give rise to liquid-like NaCl assemblies, also referred to as clusters (across the concentration range sampled in this work), that can reach substantial sizes at high concentrations.⁵⁶ Essentially, NaCl(aq) solutions become increasingly non-ideal as the solution concentration is increased. This results in small ion associates containing up to three or four ions at ∼1 M, but these clusters can encompass hundreds of ions at the high end of solution concentration, i.e., ∼10 M.

Experimental studies of levitated droplets recently found large liquid-like NaCl clusters at high supersaturations in NaCl(aq), with MD simulation results—using the same force field as the one adopted here—supporting the experimental observations.⁵⁷ The authors speculated that these clusters could play a role in NaCl crystallisation, which was later shown to be the case in simulations of high concentration metastable solutions from our own studies of homogeneous NaCl(aq),⁵⁶ and in solutions at and beyond the limit of solution stability.^58–60 Cutting-edge experiments have also demonstrated that NaCl crystals can emerge from disordered ion associates confined to aminated conical carbon nanotubes.⁶¹ In addition, amorphous NaCl solids have been isolated using supersonic spray-drying techniques, where the rapid removal of water from dense ion assemblies occurs before crystal nucleation can occur.⁶²

It was recently shown, using a sophisticated machine learning force field trained on ab initio MD simulation trajectories, that the final stage during NaCl dissolution involves the dissipation of amorphous ion clusters,⁶³ and MD simulations using the Joung–Cheatham force field also indicate this mechanism for cluster dissolution.^58,64 These studies combined, therefore, suggest that disordered ion assemblies are potentially involved in both the formation and dissolution of solid NaCl and that the Joung–Cheatham force field can capture these mechanisms reasonably well. In the presence of graphite, we demonstrated how disordered ion assemblies are stabilised in the EDL due to the increased ion densities in this region and postulated that by catalysing the formation of these clusters, surfaces might control the pathway for NaCl crystallisation.²² In this section, we explore how applied surface charges affect liquid-like ion assemblies in the EDL.

We determined the connectivity between ions in their first coordination sphere using a truncation distance based on the first minimum in the radial distribution functions between pairs of atoms (∼0.35 nm) and a continuous rational switching function to smoothly decorrelate ions according to this definition (details are provided in the files that can be obtained by the link in the Data Availability statement below). Fig. 3A shows a typical cluster residing in the first solution layer at an uncharged graphite surface when c(NaCl) is 5 M. These clusters evolve their topology over ∼ps timescales due to density fluctuations in solution.⁵⁶

Fig. 3 (A) A [Na₃Cl₄]⁻ cluster at the carbon surface surrounded by its solvation sphere. Na⁺, Cl⁻, O and H atoms are shown by blue, cyan, red and white spheres, respectively. Green lines highlight ion connections within 0.4 nm, and H-bonds are indicated by the red dashed lines. (B) Average Na–Cl coordination number (〈CN_Na–Cl〉) in the EDL as a function of surface charge density, σ. (C) Average maximum (max) CN_Na–Cl in the EDL. In B and C, ○, △ and ▽ symbols indicate data for 1, 5 and 10 M with dashed, dotted-dashed and solid lines highlighting the bulk values at 1, 5 and 10 M. (D) Number of clusters and largest cluster size in the EDL at the positive (×) and negative electrode (▽) with the bulk (○) values also provided.

Fig. 3B provides the average cation–anion coordination number, 〈CN_Na–Cl〉, in the EDL when compared to the bulk values. At 1 M, there is no clear surface effect; however, as the bulk solution concentration increases, the coordination of ions in the EDL exceeds the bulk, in line with our previous observations.²² When a surface charge is applied, Fig. 3B shows a clear bias for higher levels of ion coordination on the negative electrode. At 5 and 10 M, there is a monotonic increase in the average coordination number when σ < 0, while this remains roughly constant when σ > 0, matching bulk values at 5 M, but increasing compared to the bulk at 10 M.

The asymmetric ion coordination can be reasoned by considering the changes in the EDL atom densities as a function of σ, as provided in Fig. S5.† The affinity for cations to adsorb in the first solution layer increases as negative charges are applied to carbon atoms, and this results in increasing ρ_Cl in the vicinity of the surface, particularly as the bulk solution concentration is increased. This effect is most apparent at 10 M, where the first peaks in both ρ_Na and ρ_Cl increase substantially, in contrast to the positively charged surface, where the increase in the first ρ_Cl peak is a small fraction of that on the negative electrode and Δρ_Na for the first peak is negative. Because the highest atom densities occur close to the surface in the EDL, increasing the density in these regions means that the distribution of ions at the negative electrode is more disproportionate than at the positive electrode (see also Fig. S4†), facilitating greater ion coordination close to the surface.

Fig. 3C provides the maximum cation–anion coordination number: max(CN_Na–Cl). The plot indicates that there is a greater propensity to form ion pairs in the EDL than in the bulk solution, but there is no clear surface charge effect at the lowest concentration. Two-coordinate cations are most likely to be observed at 5 M, and this increases to three-fold cation–anion coordination in the EDL when σ ≪ 0, indicating a change from linear ion coordination to branched coordination, consistent with changes to the structure of clusters that are found on increasing ion concentrations at uncharged graphite.²² The complex, multi-layered solution structure emerging at 10 M means that very high levels of ion coordination are observed in the EDL, irrespective of the sign of the applied charge.

A value of CN_Na–Cl = 6 is consistent with the levels of coordination in the rock salt crystal. The fact that the max(CN_Na–Cl) values approach this limit in the EDL at 10 M supports the hypothesis that crystallisation is promoted in this region, with order emerging from the liquid-like clusters. Indeed, when σ = +0.58e nm⁻², a high ion density, anhydrous region of the extended cluster could potentially progress to a close-packed crystal structure. The solutions are highly metastable at this bulk concentration (using the adopted force field, the bulk solutions are metastable at 3.7–15 mol kg⁻¹ (ref. 53 and 59) which is approximately 3.5–10.9 M). Nonetheless, crystal nucleation is a rare event that is unlikely to occur over the simulation times sampled in this work.

To evaluate the size of the ion clusters, we performed a graph analysis to identify the subsets of connected components, considering ions as nodes in the graph.⁶⁵Fig. 3D provides the average largest cluster size as a function of the total number of clusters at each concentration. This confirms that ion pairs are likely to form in the EDL at 1 M. When the concentration of ions is increased to 5 M, we observed clusters in the EDL containing around four ions, and the number of these clusters substantially increases when σ < 0. At 10 M, many clusters are observed in all regions of the solution; however, the largest clusters are typically found in the EDL. Moreover, the EDL at the negatively charged graphite surface contains clusters which are significantly larger than the largest clusters observed at the positively charged surface. A snapshot of one of these extended ion networks is provided inset in Fig. 3B; this highlights the chemical heterogeneity and liquid-like ion connectivity that is typically observed in the assemblies.

Water structure at graphite

Following the evaluation of how surface charge and concentration change the structural properties of ions at the interface, in this section, we discuss how the interface affects the microscopic water structure in the EDL when compared to the bulk solution. To this aim, we evaluated variables which are functions of the positions of water O atoms that quantify the relative order of the molecules, as well as the H-bond network in the solvent. It is useful to compare these analyses of the CμMD simulations to liquid and solid forms of water. As such, additional simulations of bulk liquid water, cubic ice (ice I_c) and hexagonal ice (ice I_h), containing 4000, 2744 and 2880 molecules, respectively, were performed for 5–10 ns.

Water ordering at the interface. To quantify the local ordering of water molecules approaching the solid–liquid interface, we computed the approximate two-body excess entropy (S₂), adopting the position of oxygen atoms in water as a proxy for the centre of mass of the molecules:⁶⁶Here, k_B is Boltzmann’s constant, and ρ_Ow is the atom density in the simulation cell. g(r) is a radial pair distribution function of the distances, r, between i and j pairs of water O atoms:Here, N_Ow is the total number of water oxygen atoms. We chose r_lim to be 0.5 nm and the broadening parameter, ξ = 0.015. We obtained local averages⁶⁷ of S₂ according to,Here, f(r^ij) is a sharp but continuous switching function that identifies water molecules in the first coordination sphere according to,where r⁰ = 0.35 nm, p = 50 and q = 100.

Fig. 4A provides the approximate excess entropy probability distributions, f(S̄₂), when c(NaCl) is 5 M at the most extreme values of σ (±0.77e nm⁻²). All liquid water states are clearly separated from the S̄₂ range of values calculated for solid water phases, demonstrating that the variable adopted differentiates water molecules with different levels of local order. The presence of ions leads to a shifting of the median S̄₂ to larger values when compared to the case of pure water. In the context of the variable adopted, this suggests that the water network in solution is less ordered than in liquid water. In the EDL, the distributions are shifted to even higher values of S̄₂, particularly so in the case of the positive electrode, indicating further loss of order when compared to molecules in the bulk solution.

Fig. 4 Water structure in the bulk and EDL regions of a CμMD simulation compared to pure water (water) and hexagonal ice (I_h). (A) Water local average approximate pair entropy (S̄₂) calculated for 5 M NaCl(aq) at graphite with charge density, |σ| = 0.77e nm⁻². (B) Number of donated hydrogen bonds per water molecule (n_H). Colours blue → red indicate increasing solution charge density as shown by the scale inset. Solid lines show the mean n_H values at each concentration for water in the bulk (circles), negatively (inverted triangles) and positively (triangles) charged interface regions. Uncertainties in the data are on the scale of the size of data points. (C) H-bond lifetimes, as indicated by the time-dependent H-bond autocorrelation function (ACF). Data are provided for the case where |σ| = 0.77e nm⁻² and ion concentrations are 0 M (dashed lines; the curves for bulk and EDL σ ≤ 0 are nearly perfectly overlaid on the graph) and 10 M (solid lines). (D) Water molecule orientation as indicated by the average angle between the water molecule dipole moment (shown by the arrow on the water molecule inset) and the normal to the carbon basal surface (θ_D). The colours indicate surface charge density (see the scale in B). The data in D were smoothed using a zeroth-order Savitsky–Golay filter with a 0.05 nm window size.

In these analyses, the EDL was taken to be the region above the carbon surface encompassing only the first two solution layers; hence, this is the region of the double layer where the water structure is most perturbed when compared to the bulk. As shown in Fig. S6B,† which provides f(S̄₂) in slices throughout the entire simulation cell x-axis, the local structure of water is only significantly perturbed in the immediate vicinity of the carbon substrate. Given this observation, it is important to assess how significant the presence of ions and surface charge change the water structure, as opposed to the excluded volume effects associated with the water void space occupied by the graphite slab.

Clearly, Fig. 4A identifies a surface charge effect, but to consider the role that ions play in changing the local water order, we computed additional variables, namely the third-order Steinhardt bond orientational order parameter (q₃),⁶⁸ as well as its local (lq₃) and local average (q̄₃) values, for all water molecules comprising the bulk solution. For the functional form of q₃ and lq₃, please refer to the PLUMED documentation.⁵² These variables were previously combined with S̄₂ to identify water order in different physical states.⁶⁹ All of the distributions for these variables, provided in Fig. S6A,† indicate a small deviation from the pure water case, although this is minimal when compared to ice and amorphous water (liquid water crash cooled to 100 K during a 10 ns simulation). From this analysis, it would appear that the surfaces have a greater effect on the local ordering of water molecule centres than the presence of structure-breaking ions.

Water H-bonds. The intermolecular structure of water is usually described in terms of the H-bond network. We, therefore, calculated—using MDAnalysis^70,71—the number of donated H-bonds per water molecule, n_H, in the bulk and EDL regions of all simulations using a simple geometric criterion for these bonds of the type, O_D—H⋯O_A, where D and A subscripts refer to the H-bond donor and acceptor, respectively. H-Bonds were assigned when O_D and O_A were within 0.3 nm and the angle, ∠O_DHO_A > 150°. Fig. S7† provides the distributions for H-bond distances and angles in the EDL and bulk regions of CμMD simulations, as well as for ice I_h and pure liquid water.

Fig. 4B shows that irrespective of the applied surface charge, the mean n_H for water in the bulk and EDL regions of CμMD simulations in the absence of ions are close to the values in pure liquid water. The small deviation from the homogeneous liquid case is likely to result from the excluded volume effects associated with the selection of water molecules in regions of the simulation cell x-axis that creates (artificial) excluded volumes, even in the case of what we describe as the bulk. Slightly more than one H-bond per water molecule is observed, which is understandably lower than the expected value of two for ice I_h.

In contrast to the structural variables discussed above, n_H is far more sensitive to the solution concentration and less sensitive to the magnitude and sign of σ. At 1 M, we find a difference between n_H in the EDL when compared to the bulk by around n_H = 0.1. This difference was approximately the same at all levels of concentration; however, n_H in all regions of the simulations decreased as the concentration of ions increased. At 10 M, n_H is approximately 0.35–0.4, suggesting a near complete breaking of the H-bond structure at the highest concentrations. The presence of ions and their assemblies, as well as associated local electric fields, greatly perturbs the water structure from the pure solvent, which ultimately has implications for these systems and their performance as conductors of electrical charge.

The lifetime for H-bonds was determined according to the autocorrelation function (ACF):

where τ is a time lag in the data, t₀ indicates a time origin and H_ij signifies the presence of an assigned H-bond between molecules i and j, taking a binary value of zero or one according to the H-bond distance and angle cut-offs described above. The H-bond lifetimes are provided in Fig. 4C for pure liquid water, ice and the bulk and EDL regions of simulations at 0 and 10 M with σ = ±0.77e nm⁻².

In the absence of ions, the H-bond lifetime for water in the bulk and at the negative electrode is identical to the lifetime of H-bonds in pure water. At the positive electrode, however, the H-bond lifetime is extended upon the application of a large surface charge. This is most likely due to the increased binding strength of water O atoms to the charged graphite surface. At 10 M, we find that the H-bond lifetimes in the bulk are extended cf. 0 M, and there is a divergence in the lifetimes in the bulk and at the negative electrode. These lifetimes, however, are still lower than the average H-bond lifetime at the positive electrode, which can extend to ∼10² ps. Simulations indicate that the orientation of water molecules at charged surfaces determines the propensity for heterogeneous ice crystallisation, with positively charged silver iodide surfaces suggested to promote ice nucleation.⁷² A stabilisation of the H-bond lifetime at positively charged surfaces potentially has additional ramifications on the ability of these charged substrates to promote ice nucleation. In light of our findings, it would be useful to test how these phenomena and the presence of ions in solution control ice crystallisation rates.

Water dipole moments. We characterised the orientation of water molecules with respect to the graphite surface plane by computing the angle between water molecule dipole moments and the normal to the surface, θ_D. The calculation was implemented such that if the dipole moment was perfectly perpendicular to the surface normal and pointing towards the surface (see the water dipole moment arrow inset of Fig. 4D) the value of θ_D is zero and 180° if the dipole moment vector points away from the surface. θ_D = 90°, therefore, indicates water molecules with dipole moments that are, on average, aligned parallel with the basal surface, and/or indicates no preference for the orientation of water molecules at the surface.

Fig. 4D provides the mean θ_D as a function of Δx. Before discussing the effect of surface charge and solution concentration on these results, it is useful to discuss an important feature of the curves at 0 M. In this subset of simulations, no ions are present, and we did not adopt CμMD to study the effect of charged surfaces on the water structure. As such, there is no ionic reservoir in this system. When equal but opposite signs of surface charge are applied to opposite faces of the graphite slab, electric fields are induced that span the periodic boundaries in x. This effect is evident in the θ_D curves at 0 M, where the neutral graphite case (see the symmetric blue curves at both surfaces) indicates that θ_D = 90° in the bulk, but where positive and negative deviations in the mean θ_D are found at the positive and negative electrode, respectively. It is possible to avoid these electrical artefacts by removing the periodic boundaries in x and/or by including artificial electrical insulating layers parallel to the graphite surface. For the purposes of this study, this was not necessary, and, importantly, these effects are removed by the presence of the ionic reservoir. In terms of the discussion of θ_D at 0 M, we compare the features of the distributions in the EDL with respect to the values in the bulk in order to understand how the surface controls the orientation of water; furthermore, we do not believe that the dipole associated with the graphite slab has a significant effect on the analysis of EDL solution thermodynamics, discussed in the following section.

The θ_D distributions indicate that the surface, in the absence of ions, leads to no significant ordering of water molecule dipole moments with respect to the surface plane, as expected with the adopted force field³⁸ and shown from simulations elsewhere.²⁴ Furthermore, negative surface charges give rise to increased θ_D in the first solution layer(s), confirming that water molecules tend to point away from the surface (compared to the bulk) in order to maximise the interactions between H atoms and the surface. At the positive electrode, a maximum occurs at Δx ≈ −0.5 nm; here, there is a depletion in the water density (see Fig. S1†), which perhaps allows water to restructure more freely to screen the surface potential.

In the presence of ions, θ_D angles are obtuse in the first two solution layers, regardless of the sign of the surface charge; however, maximum values are observed on the negative electrode. In all curves, a minimum occurs at Δx ≈ −0.65 nm, which is where the second peak in water density profiles is observed when moving away from the surface. An acute θ_D was also observed for 1 M solutions in contact with graphene using quantum mechanical MD simulations.²⁴ This result indicates that the mean orientation of water dipoles in the first and second solution layers differs by 40–50°. The complex EDL structure is evident in the θ_D curves, notably at 10 M, where large fluctuations in the angle distributions are evident within 0.5 nm from graphite. Finally, very little ordering of water is observed beyond ∼2 nm from the substrate at the highest concentrations.

Solution thermodynamics

As well as the capability to undergo rapid charge/discharge cycling, carbon–electrolyte interfaces can be exploited for applications to promote chemical reactions.^6,7 Understanding the thermodynamic properties of the interface is essential in this regard. In the following section, we characterise the electrochemical properties of the interface, with a focus on the solution side of the EDL.

Electric potential at the interface. The capacity for the interface to store charge is C = σ/Δψ⁰, where Δψ⁰ is the electric potential change (usually termed the ‘potential drop’) across the interface with an applied surface charge minus the potential drop in the absence of a surface charge. Poisson’s equation relates the electric potential to the charge density (ρ_q) according to,where E is the electric field and ε is the permittivity of the medium. Here we take ε = ε₀, the permittivity of free space, because the full solution charge density is used in the calculation. It is important to recognise that this analysis provides the capacitance associated with the solution side of the EDL in contact with a uniformly charged surface. An additional contribution to the total capacitance comes from the density of electron states in the substrate, which represents a minor contribution to the interfacial capacitance under the conditions studied.^23,32

Fig. 5 provides ψ(x) for all systems and all applied charges. In the case of pure water, we find that the potential drop across the interface, defined as Δψ = ψ(Δx = 0) − ψ^b (where ψ^b is the electrostatic potential in the bulk), in the absence of surface charge is 0.23 V. As the surface charge is increased to ± 0.77e nm⁻², Δψ was calculated as −1.2 V and 2.6 V on the negative and positive electrode, respectively. Two minima are observed in Fig. 5A, in accordance with the maxima in the water density profiles (see Fig. 1). The application of surface charge perturbs these densities, as discussed above, which changes the position of the first minima in ψ(x) (from the surface) and makes the second minimum shallow when σ < 0.

Fig. 5 Interface electric potential as a function of distance from the electrode, Δx. Colours blue → red indicate increasing solution charge densities from |σ| = 0 → 0.77e nm⁻². (A)–(D) are taken from simulations where the target bulk ion concentration was 0, 1, 5 and 10 M, respectively.

The presence of ions induces additional fluctuations in ψ(x) when compared to the pure water case. These extend to around 1.5 nm from the carbon surface at the highest concentrations, although the amplitude of the fluctuations is not particularly correlated with concentration, due to the fact that local electric fields are determined by the total solution charge density. Some notable features of the ψ(x) curves are the fact that cation adsorption in the absence of surface charge increases Δψ to ∼0.42 V, in good agreement with studies of aqueous solutions at carbon surfaces using different models for the interface.⁸ In addition, a minimum emerges at Δx ≈ −0.3 nm at the positive electrode, associated with a depletion of cations and accumulation of anions. Furthermore, a deep minimum is observed at 10 M when Δx ≈ −0.4 nm on the negative electrode, which can be attributed to the increasing anion concentration and water restructuring that occurs at this surface.

When σ = +0.77e nm⁻², Δψ at the positive electrode ranges from 2.85–2.9 V across the range of concentrations investigated. Using the expression for capacitance reported above, this equates to a solution-side capacitance of ∼5 μF cm⁻². Δψ at the negatively charged electrode, instead, ranges from −0.92 to −1.34 V, corresponding to a capacitance of C = 7–9 μF cm⁻², indicating that graphite has a greater capacity to store charge when negative charges are applied. C evaluated in these simulations is consistent with estimates from experiments.²² Moreover, the increased capacity to store ionic charge at the negative electrode is consistent with results elsewhere.^23,25 A result worthy of note is that the concentration of ions has very little effect on the ability of graphite to store charge over the range of molarities considered, which was also found for graphene.²³ Indeed, increasing the charge capacity of the carbon–solution interface typically requires tuning the properties of the solute and the substrate to increase the overall capacitance of the system, rather than simply changing the concentration of the solution.^14,73

Ion activity coefficients. The chemical potential of species i, μ_i, is defined as the change in free energy associated with a variation in the number of i molecules and represents the ability of that species to undergo a physical–chemical transformation. In the presence of an electric field, when i is a charged species, the same information is captured by the electrochemical potential, which accounts for the additional energetic contributions to insert/remove a charged particle to/from the system:In the above equation, μ⁰ is a reference chemical potential. The second term on the right of the equation provides the energy associated with particle exchange in non-ideal solutions, where R, T and a indicate the gas constant, temperature and solute activity, respectively. This term can be expanded to account for the ideal and excess chemical potential, which are functions of the total concentration, in this case, the solution molality, m (strictly, this is a unitless quantity defining the mole fraction of solute in solution compared to the standard state of 1 mol kg⁻¹), and the activity coefficient, γ. The final term defines the work to transfer a particle with charge z into the system with electrostatic potential, ψ. Faraday’s constant, F, ensures that the term has the correct energy units.

In order to determine [small mu, Greek, tilde]_i for ions in our simulations, an activity model is required. Zimmermann et al. provided an analytical formula to calculate μ for ions in NaCl(aq) as a function of ion molality, m_ion, by fitting to simulation data:^74–76

μ_ion = μ⁰_ion + 2 RT ln m_ion + 2RT ln γ_ion (9)

where,In these equations, μ⁰_ion = −391.6 kJ mol⁻¹,⁷⁵a = 0.568 mol^−1/2 kg^1/2, b = 1.17769 mol^1/2 kg^−1/2 and c = 0.177157 mol⁻¹ kg. It is important to recognise that this activity model assumes changes to the solution density and dielectric constant are only a function of the solution composition. This model can be extended²² to account for the effect of electric fields and associated varying ion molalities that occur on approach to the graphite surface according to,

[small mu, Greek, tilde]_ion(x) = μ⁰_ion + RT ln m_Na(x) + RT ln γ_ion(m_Na(x)) + ωFψ(x) + RT ln m_Cl(x) + RT ln γ_ion(m_Cl(x)) − (1 − ω)Fψ(x) (11)

where subscript labels indicate Na⁺ or Cl⁻ molalities and ω(x) = m_Na(x)/(m_Na(x) + m_Cl(x)). In the limiting case where the molalities of cations and anions are equal locally, eqn (11) reduces to eqn (9).

For the case of 1 M NaCl(aq), we determined [small mu, Greek, tilde]_ion(x) ≈ −392 kJ mol⁻¹ in the bulk where ψ(x) = 0, in good agreement with the expected chemical potential from the model by Zimmermann et al.⁷⁶ for homogeneous solutions with m_ion ≈ 1.3 mol kg⁻¹. γ_ion = 0.9 under these conditions, which we approximate to a value of one for the subsequent analyses, such that [small mu, Greek, tilde]^b_ion ≈ μ⁰_ion, with [small mu, Greek, tilde]^b_ion representing the electrochemical potential of ions in the bulk. The energy change associated with 2RT ln γ_ion when γ_ion = 0.9 is ∼0.5 kJ mol⁻¹.

In our simulations, the chemical potential of ions and water as a function of x in the steady state is constant.³¹ Therefore, with knowledge of the ion molalities and electric potential in the EDL, eqn (11) can be rearranged to determine how the presence of the surface—where the density of ions and dielectric constant of the solution are changing compared with the bulk—affects the activity of ions as captured by γ_ion:

−RT ln[γ_ion(m_Na(x))γ_ion(m_Cl(x))] = RT[ln m_Na(x)m_Cl(x)] + (2ω − 1)Fψ(x) (12)

We label this quantity Δ[small mu, Greek, tilde]^E_ion. Fig. 6A provides −Δ[small mu, Greek, tilde]^E_ion/RT = 2 ln(γ_ion) at graphite with varying charge density when c(NaCl) = 1 M. On approach to the surface, there is a small minimum around Δx = −0.8 nm when σ = 0, which is consistent with the position of a second cation-rich solution layer above the surface. This is followed by a gradual decrease in 2 ln(γ_ion) towards a second minimum around Δx = −0.4 nm. This minimum resides between the maxima for the densities of the first cation- and anion-rich solution layers in the EDL (see Fig. S2†).

Fig. 6 (A) The excess ion electrochemical potential −Δμ^E_ion/RT = 2 ln(γ_ion) in solution (1 M) as a function of distance from the graphite basal surface (Δx). (B–D) The excess water chemical potential −Δμ^E_wat/RT = ln(γ(x)/γ^b) for water molecules as a function of Δx when the target bulk solution concentration was 1, 5 and 10 M, respectively. In all plots, the top and bottom panels provide data for solutions at negatively and positively charged surfaces, respectively; blue → red colours indicate |σ| = 0–0.77e nm⁻². Green dashed lines highlight the position of the first two minima in −Δμ^E_ion/RT when σ = 0.

When a positive charge is applied to graphite, as shown in the bottom panel of Fig. 6A, 2 ln(γ_ion) becomes more negative to around −25RT; this is due to the positive surface charge pushing and pulling Na⁺ and Cl⁻ away from and towards the surface (see Fig. S5†). In addition, there is a small increase in 2 ln(γ_ion) at Δx = −0.65 nm, due to a relatively high mole fraction of Na⁺ in this region. On the negatively charged surface, the applied potential displaces Cl⁻, which leads to a decrease in the anion density at Δx = −0.4 nm and an increase at Δx = −0.55 nm compared to the case when σ = 0 (see Fig. S5†); this results in positive and negative increases to 2 ln(γ_ion), respectively.

Fig. S8† provides the contributions to Δ[small mu, Greek, tilde]_ion(x) for a single case where c = 1 M and σ = +0.77e nm⁻². This shows that it is essential to account for the effect of the surface excluded volume and charge on the structure of the solution when determining the activities of ions. In particular, in the presence of an applied surface potential, the contribution to Δ[small mu, Greek, tilde]_ion(x) from the electric potential drop can be as significant as the changing mole fraction of solute at the surface.

Water activity coefficients. We now consider how charged surfaces affect water activity coefficients in the EDL when compared to the bulk. We note that, for an electrically neutral species such as water, the electrochemical potential reduces to the chemical potential even in the presence of a charged surface. To determine the chemical potential of SPC/E water in NaCl(aq) without an accurate activity model, we make use of the fact that in our CμMD simulations, the chemical potential for water molecules in the EDL and bulk are equal. We can estimate the potential of mean force () to transfer water molecules from the bulk to the EDL according to,where p_wat represents the probability density of observing water molecules at position x, and the superscript b indicates the probability density at a point in x representative of the bulk solution. As such, provides a proxy for ΔA: the Helmholtz free energy change for the transformation under question.

Given that ΔA = nΔμ_wat = −RT ln K, where n is the number of moles of water and K = a_wat(x)/a^b_wat (where a indicates activity) which is the equilibrium constant for the transfer of one water molecule from position x to the solution bulk, we can also write,

where χ_wat and γ_wat are the mole fraction and activity coefficient for water molecules in solution, respectively. Hence, by combining eqn (13) and (14) we can evaluate Δμ^E_wat, which indicates how γ_wat changes in comparison to γ^b_wat.

Fig. 6B–D provides −Δμ^E_wat(x)/RT for systems where the ion concentration varies from 1 to 10 M. At 1 M, two minima are observed in ln(γ_wat/γ^b_wat) at Δx = −0.5 and −0.75 nm that are within RT of the bulk value. The contributions to in the case where c = 1 M and σ = +0.77e nm⁻² are provided in Fig. S9,† which indicate that the first minimum from the surface arises due to a relatively high value of Δμ^I_wat(x) in this region, associated with a decreased water mole fraction, as well as effects associated with the electric fields close to the carbon basal plane. The minimum at Δx = −0.75 nm, however, occurs in a region where the water mole fraction is not significantly different from the bulk value and is due to the structuring of ions and electric fields locally. As the magnitude of the applied charge increases, small shifts to the position of the minima occur, associated with changes to the water density. The maximum in ln(γ_wat/γ^b_wat) occurs around Δx = −0.3 nm; this conforms to the minimum in , and represents the first water layer adsorbed at the graphite surface (see Fig. S2 and S9†). This indicates that γ_wat is 3–4 greater in the innermost EDL solution layer than in the bulk.

As the concentration of ions is increased, the positions of the maxima and minima in −Δμ^E_wat(x)/RT are unchanged. Additional fluctuations beyond Δx ≈ −1 nm are observed at the highest concentrations, which are only partly associated with changes to Δμ^I_wat(x) (see Fig. S10†). At 10 M, the most negative minimum in Fig. 6D suggests that γ_wat/γ^b_wat = 0.2. It is also apparent, at the highest concentration, that an additional minimum in −Δμ^E_wat(x)/RT occurs around Δx = −0.35 nm. This can be attributed to an increase in Δμ^I_wat when compared with lower concentrations, concomitant with a decrease in χ_wat compared with χ^b_wat. Changes to the local structure of the solution on increasing surface charge density tend to be greatest at the highest bulk concentrations (see Fig. S5†), so it is perhaps not surprising that a richer behaviour in the Δμ^E_wat curves emerges at 10 M as a function of σ.

Conclusions

Understanding the properties of the carbon–electrolyte interface is important for a range of applications of these systems to facilitate, e.g., charge storage and chemical reactions. The CμMD simulations we have performed in this work provide an atomic scale resolution of the interface of graphite with NaCl(aq) where the concentration of ions and surface charge was varied. The simulations allow us to investigate how the asymmetric but cooperative adsorption of ions in the EDL, under the effect of an applied potential, affects the structure, dynamics and thermodynamic properties of the interface at a constant thermodynamic driving force for adsorption (defined by the chemical potential of ions in the bulk solution).

Our simulations indicate that increasing the magnitude of the surface charge is analogous to increasing the concentration of ions in solution; both changes give rise to a complex, multi-layered solution structure comprising cation- and anion-rich solution layers due to the finite size of ions and the partial saturation of solution layers with ions, that is not readily captured by mean-field models of the EDL. Perturbations to the solution structure typically extend 1–2 nm from the surface. Interestingly, the presence of a relatively low concentration of ions decreases the intensity of fluctuations in water densities in the EDL compared to the case where no ions are present.

Liquid-like NaCl clusters have been observed in bulk NaCl(aq) solutions at relatively high concentrations,^56,57 and our previous work on graphite identified that these clusters are stabilised in the EDL.²² In the present study, we demonstrated that negative surface charges increase the number and size of these networks, which can include tens of ions at moderate supersaturations. Furthermore, at negative electrodes, the local ion density in these networks increases, as indicated by changes to the average cation–anion coordination number. This result raises important questions regarding the ability of charged surfaces to induce NaCl crystallisation.

Our analyses of water structure in the EDL and the lifetime of H-bonds indicate that positive electrodes can induce the reorientation of water molecules at the surface and increase the lifetime of H-bonded networks in this region. These effects, in turn, have potential implications for the crystallisation of ice in the presence of charged carbon substrates, as well as for the role that these interfaces play in catalysis. It is important to note that the water model we adopt is constrained to its equilibrium, bulk liquid water geometry and the partial charges on O and H atoms are fixed. Future studies should consider how constrained geometry, non-polarisable water models affect the trends found in this work, and how different models for carbon–water interactions affect the thermodynamic properties evaluated here.

Although our analysis of the electrical properties of the interface indicates a small increase in the capacity of the negative electrode to store charge, the difference in capacitance was 2–4 μF cm⁻², with solution concentration playing only a small role in increasing the ability for the negative electrode to accumulate ions. This is perhaps not surprising because, at molar concentrations, the affinity of graphite for cations and the cooperative adsorption of anions lends the first solution layers already partially saturated with ions in the absence of surface charge. Analysis of how water and ion activity coefficients deviate from the bulk values, when the electric potential in the EDL is accounted for, indicates how the excess chemical potential for water decreases in the first solution layers adjacent to the surface, concomitant with an increase in the water density in this region, as well as how the changing ion densities induce fluctuations in water activity ratios in the EDL as the concentration of ions increases.

In summary, the complex interplay of solution concentration and surface charge effects provides a picture of the EDL that is difficult to obtain in experiments and from mean-field models. We hope that the questions raised in this work provide inspiration for further simulation and experimental studies of this system.

Data availability

GROMACS input and example output files, including the force field parameters necessary to reproduce the simulation results reported in this paper, are available on github (see https://github.com/aaronrfinney/CmuMD-NaCl_at_graphite). The PLUMED input files are also accessible via PLUMED-NEST (https://www.plumed-nest.org (ref. 77)), the public repository for the PLUMED consortium, using the project ID, plumID: 23.027. Details on how to use and implement the CμMD method within PLUMED are available on github (see https://github.com/mme-ucl/CmuMD).

Author contributions

A. R. F. and M. S. designed the research. A. R. F. performed the research and analyses. A. R. F. and M. S. wrote and edited the paper.

Conflicts of interest

The authors declare no conflict of interest.

Acknowledgements

The authors acknowledge funding from the Crystallisation in the Real World EPSRC Programme grant (grant EP/R018820/1) and the ht-MATTER UKRI Frontier Research Guarantee grant (EP/X033139/1). The authors acknowledge the use of the UCL Myriad High Throughput Computing Facility (Myriad@UCL), and associated support services, in the completion of this work.

References

1. E. Frackowiak and F. Béguin, Carbon materials for the electrochemical storage of energy in capacitors, Carbon, 2001, 39, 937–950.
2. Y. Wang, Y. Song and Y. Xia, Electrochemical capacitors: mechanism, materials, systems, characterization and applications, Chem. Soc. Rev., 2016, 45, 5925–5950.
3. P. Simon and Y. Gogotsi, Materials for electrochemical capacitors, Nat. Mater., 2008, 7, 845–854.
4. S. Porada, R. Zhao, A. Van Der Wal, V. Presser and P. Biesheuvel, Review on the science and technology of water desalination by capacitive deionization, Prog. Mater. Sci., 2013, 58, 1388–1442.
5. D. Cohen-Tanugi and J. C. Grossman, Water Desalination across Nanoporous Graphene, Nano Lett., 2012, 12, 3602–3608.
6. H. Tian, J. Liang and J. Liu, Nanoengineering Carbon Spheres as Nanoreactors for Sustainable Energy Applications, Adv. Mater., 2019, 31, 1903886.
7. N. M. Julkapli and S. Bagheri, Graphene supported heterogeneous catalysts: an overview, Int. J. Hydrogen Energy, 2015, 40, 948–979.
8. J. D. Elliott, A. A. Papaderakis, R. A. W. Dryfe and P. Carbone, The electrochemical double layer at the graphene/aqueous electrolyte interface: what we can learn from simulations, experiments, and theory, J. Mater. Chem. C, 2022, 10, 15225–15262.
9. D. N. Petsev, F. van Swol and L. J. Frink, Molecular Theory of Electric Double Layers, IOP Publishing Ltd, 2021, OCLC: 1280155254.
10. A. J. Bard and L. R. Faulkner, Electrochemical Methods: Fundamentals and Applications, Wiley, New York, 2nd edn, 2001.
11. P. Iamprasertkun, W. Hirunpinyopas, A. Keerthi, B. Wang, B. Radha, M. A. Bissett and R. A. W. Dryfe, Capacitance of Basal Plane and Edge-Oriented Highly Ordered Pyrolytic Graphite: Specific Ion Effects, J. Phys. Chem. Lett., 2019, 10, 617–623.
12. C. Zhan, M. R. Cerón, S. A. Hawks, M. Otani, B. C. Wood, T. A. Pham, M. Stadermann and P. G. Campbell, Specific ion effects at graphitic interfaces, Nat. Commun., 2019, 10, 4858.
13. D. C. Grahame, The Electrical Double Layer and the Theory of Electrocapillarity, Chem. Rev., 1947, 41, 441–501.
14. M. V. Fedorov and A. A. Kornyshev, Ionic Liquids at Electrified Interfaces, Chem. Rev., 2014, 114, 2978–3036.
15. I. Borukhov, D. Andelman and H. Orland, Steric Effects in Electrolytes: A Modified Poisson–Boltzmann Equation, Phys. Rev. Lett., 1997, 79, 435–438.
16. Z. A. Goodwin, G. Feng and A. A. Kornyshev, Mean-Field Theory of Electrical Double Layer In Ionic Liquids with Account of Short-Range Correlations, Electrochim. Acta, 2017, 225, 190–197.
17. L. Yin, Y. Huang, H. Chen and T. Yan, A mean-field theory on the differential capacitance of asymmetric ionic liquid electrolytes. II. Accounts of ionic interactions, Phys. Chem. Chem. Phys., 2018, 20, 17606–17614.
18. Y. Uematsu, R. R. Netz and D. J. Bonthuis, The effects of ion adsorption on the potential of zero charge and the differential capacitance of charged aqueous interfaces, J. Phys.: Condens. Matter, 2018, 30, 064002.
19. J. G. Hedley, H. Berthoumieux and A. A. Kornyshev, The Dramatic Effect of Water Structure on Hydration Forces and the Electrical Double Layer, J. Phys. Chem. C, 2023, 127, 8429–8447.
20. M. McEldrew, Z. A. H. Goodwin, S. Bi, A. A. Kornyshev and M. Z. Bazant, Ion Clusters and Networks in Water-in-Salt Electrolytes, J. Electrochem. Soc., 2021, 168, 050514.
21. Z. A. H. Goodwin, M. McEldrew, J. Pedro De Souza, M. Z. Bazant and A. A. Kornyshev, Gelation, clustering, and crowding in the electrical double layer of ionic liquids, J. Chem. Phys., 2022, 157, 094106.
22. A. R. Finney, I. J. McPherson, P. R. Unwin and M. Salvalaglio, Electrochemistry, ion adsorption and dynamics in the double layer: a study of NaCl(aq) on graphite, Chem. Sci., 2021, 12, 11166–11180.
23. N. Di Pasquale, A. R. Finney, J. D. Elliott, P. Carbone and M. Salvalaglio, Constant Chemical Potential-Quantum Mechanical-Molecular Dynamics simulations of the Graphene-electrolyte double layer, J. Chem. Phys., 2023, 158, 134714.
24. J. D. Elliott, A. Troisi and P. Carbone, A QM/MD Coupling Method to Model the Ion-Induced Polarization of Graphene, J. Chem. Theory Comput., 2020, 16, 5253–5263.
25. J. Dočkal, F. Moučka and M. Lísal, Molecular Dynamics of Graphene–Electrolyte Interface: Interfacial Solution Structure and Molecular Diffusion, J. Phys. Chem. C, 2019, 123, 26379–26396.
26. J. Dočkal, M. Lísal and F. Moučka, Molecular dynamics of the interfacial solution structure of alkali-halide electrolytes at graphene electrodes, J. Mol. Liq., 2022, 353, 118776.
27. A. A. Kornyshev, Double-Layer in Ionic Liquids: Paradigm Change?, J. Phys. Chem. B, 2007, 111, 5545–5557.
28. M. V. Fedorov and A. A. Kornyshev, Ionic Liquid Near a Charged Wall: Structure and Capacitance of Electrical Double Layer, J. Phys. Chem. B, 2008, 112, 11868–11872.
29. W. Schmickler, Interfacial Electrochemistry, Oxford University Press, New York, 1996.
30. C. Perego, M. Salvalaglio and M. Parrinello, Molecular dynamics simulations of solutions at constant chemical potential, J. Chem. Phys., 2015, 142, 144113.
31. T. Karmakar, A. R. Finney, M. Salvalaglio, A. O. Yazaydin and C. Perego, Non-Equilibrium Modeling of Concentration-Driven processes with Constant Chemical Potential Molecular Dynamics Simulations, Acc. Chem. Res., 2023, 56, 1156–1167.
32. A. R. Finney and M. Salvalaglio, Bridging the gap between mesoscopic and molecular models of solid/liquid interfaces out-of-equilibrium, Chem. Eng. Res. Des., 2022, 180, 285–295.
33. J.-F. Olivieri, J. T. Hynes and D. Laage, Confined Water’s Dielectric Constant Reduction Is Due to the Surrounding Low Dielectric Media and Not to Interfacial Molecular Ordering, J. Phys. Chem. Lett., 2021, 12, 4319–4326.
34. I. S. Joung and T. E. Cheatham, Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations, J. Phys. Chem. B, 2008, 112, 9020–9041.
35. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, The missing term in effective pair potentials, J. Phys. Chem., 1987, 91, 6269–6271.
36. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids, J. Am. Chem. Soc., 1996, 118, 11225–11236.
37. J. Ma, A. Michaelides, D. Alfè, L. Schimka, G. Kresse and E. Wang, Adsorption and diffusion of water on graphene from first principles, Phys. Rev. B, 2011, 84, 033402.
38. Y. Wu and N. R. Aluru, Graphitic Carbon–Water Nonbonded Interaction Parameters, J. Phys. Chem. B, 2013, 117, 8802–8813.
39. T. Werder, J. H. Walther, R. L. Jaffe, T. Halicioglu and P. Koumoutsakos, On the Water–Carbon Interaction for Use in Molecular Dynamics Simulations of Graphite and Carbon Nanotubes, J. Phys. Chem. B, 2003, 107, 1345–1352.
40. H.-Y. Kim, M. C. Dos Santos and M. W. Cole, Wetting Transitions of Water on Graphite and Graphene, J. Phys. Chem. A, 2014, 118, 8237–8241.
41. J. Li and F. Wang, Water graphene contact surface investigated by pairwise potentials from force-matching PAW-PBE with dispersion correction, J. Chem. Phys., 2017, 146, 054702.
42. M.E. Schrader, Ultrahigh vacuum techniques in the measurement of contact angles. IV. Water on graphite (0001), J. Phys. Chem., 1975, 79, 2508–2515.
43. M. E. Schrader, Ultrahigh-vacuum techniques in the measurement of contact angles. 5. LEED study of the effect of structure on the wettability of graphite, J. Phys. Chem., 1980, 84, 2774–2779.
44. A. V. Prydatko, L. A. Belyaeva, L. Jiang, L. M. C. Lima and G. F. Schneider, Contact angle measurement of free-standing square-millimeter single-layer graphene, Nat. Commun., 2018, 9, 4185.
45. J. G. Brandenburg, A. Zen, M. Fitzner, B. Ramberger, G. Kresse, T. Tsatsoulis, A. Grüneis, A. Michaelides and D. Alfè, Physisorption of Water on Graphene: Subchemical Accuracy from Many-Body Electronic Structure Methods, J. Phys. Chem. Lett., 2019, 10, 358–368.
46. C. D. Williams, J. Dix, A. Troisi and P. Carbone, Effective Polarization in Pairwise Potentials at the Graphene–Electrolyte Interface, J. Phys. Chem. Lett., 2017, 8, 703–708.
47. B. Hess, C. Kutzner, D. Van Der Spoel and E. Lindahl, GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation, J. Chem. Theory Comput., 2008, 4, 435–447.
48. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, LINCS: a linear constraint solver for molecular simulations, J. Comput. Chem., 1997, 18, 1463–1472.
49. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, A smooth particle mesh Ewald method, J. Chem. Phys., 1995, 103, 8577–8593.
50. G. Bussi, D. Donadio and M. Parrinello, Canonical sampling through velocity rescaling, J. Chem. Phys., 2007, 126, 014101.
51. H. J. C. Berendsen, J. P. M. Postma, W. F. Van Gunsteren, A. DiNola and J. R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys., 1984, 81, 3684–3690.
52. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, PLUMED 2: new feathers for an old bird, Comput. Phys. Commun., 2014, 185, 604–613.
53. A. L. Benavides, J. L. Aragones and C. Vega, Consensus on the solubility of NaCl in water from computer simulations using the chemical potential route, J. Chem. Phys., 2016, 144, 124504.
54. A. Coretti, C. Bacon, R. Berthin, A. Serva, L. Scalfi, I. Chubak, K. Goloviznina, M. Haefele, A. Marin-Laflèche, B. Rotenberg, S. Bonella and M. Salanne, MetalWalls: simulating electrochemical interfaces between polarizable electrolytes and metallic electrodes, J. Chem. Phys., 2022, 157, 184801.
55. T. Funabashi, Integration of Distributed Energy Resources in Power Systems, Elsevier, 2016, pp. 1–14.
56. A. R. Finney and M. Salvalaglio, Multiple pathways in NaCl homogeneous crystal nucleation, Faraday Discuss., 2022, 235, 56–80.
57. H. Hwang, Y. C. Cho, S. Lee, Y.-H. Lee, S. Kim, Y. Kim, W. Jo, P. Duchstein, D. Zahn and G. W. Lee, Hydration breaking and chemical ordering in a levitated NaCl solution droplet beyond the metastable zone width limit: evidence for the early stage of two-step nucleation, Chem. Sci., 2021, 12, 179–187.
58. G. Lanaro and G. N. Patey, Birth of NaCl Crystals: Insights from Molecular Simulations, J. Phys. Chem. B, 2016, 120, 9076–9087.
59. H. Jiang, P. G. Debenedetti and A. Z. Panagiotopoulos, Nucleation in aqueous NaCl solutions shifts from 1-step to 2-step mechanism on crossing the spinodal, J. Chem. Phys., 2019, 150, 124502.
60. P. S. Bulutoglu, S. Wang, M. Boukerche, N. K. Nere, D. S. Corti and D. Ramkrishna, An investigation of the kinetics and thermodynamics of NaCl nucleation through composite clusters, PNAS Nexus, 2022, 1, pgac033.
61. T. Nakamuro, M. Sakakibara, H. Nada, K. Harano and E. Nakamura, Capturing the Moment of Emergence of Crystal Nucleus from Disorder, J. Am. Chem. Soc., 2021, 143, 1763–1767.
62. E. Amstad, M. Gopinadhan, C. Holtze, C. O. Osuji, M. P. Brenner, F. Spaepen and D. A. Weitz, Production of amorphous nanoparticles by supersonic spray-drying with a microfluidic nebulator, Science, 2015, 349, 956–960.
63. N. O’Neill, C. Schran, S. J. Cox and A. Michaelides, Crumbling Crystals: On the Dissolution Mechanism of NaCl in Water, arXiv, 2022, preprint, arXiv:2211.04345 [cond-mat, physics:physics],  DOI:10.48550/arXiv.2211.04345.
64. G. Lanaro and G. N. Patey, Molecular Dynamics Simulation of NaCl Dissolution, J. Phys. Chem. B, 2015, 119, 4275–4283.
65. G. A. Tribello, F. Giberti, G. C. Sosso, M. Salvalaglio and M. Parrinello, Analyzing and Driving Cluster Formation in Atomistic Simulations, J. Chem. Theory Comput., 2017, 13, 1317–1327.
66. P. M. Piaggi, O. Valsson and M. Parrinello, Enhancing Entropy and Enthalpy Fluctuations to Drive Crystallization in Atomistic Simulations, Phys. Rev. Lett., 2017, 119, 015701.
67. W. Lechner and C. Dellago, Accurate determination of crystal structures based on averaged local bond order parameters, J. Chem. Phys., 2008, 129, 114707.
68. P. J. Steinhardt, D. R. Nelson and M. Ronchetti, Bond-orientational order in liquids and glasses, Phys. Rev. B, 1983, 28, 784–805.
69. M. Fulford, M. Salvalaglio and C. Molteni, DeepIce: A Deep Neural Network Approach To Identify Ice and Water Molecules, J. Chem. Inf. Model., 2019, 59, 2141–2149.
70. R. Gowers, M. Linke, J. Barnoud, T. Reddy, M. Melo, S. Seyler, J. Domański, D. Dotson, S. Buchoux, I. Kenney and O. Beckstein, MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations, Austin, Texas, 2016, pp. 98–105.
71. R. J. Gowers and P. Carbone, A multiscale approach to model hydrogen bonding: the case of polyamide, J. Chem. Phys., 2015, 142, 224907.
72. B. Glatz and S. Sarupria, The surface charge distribution affects the ice nucleating efficiency of silver iodide, J. Chem. Phys., 2016, 145, 211924.
73. H. Ji, X. Zhao, Z. Qiao, J. Jung, Y. Zhu, Y. Lu, L. L. Zhang, A. H. MacDonald and R. S. Ruoff, Capacitance of carbon-based electrical double-layer capacitors, Nat. Commun., 2014, 5, 3317.
74. F. Moučka, I. Nezbeda and W. R. Smith, Molecular simulation of aqueous electrolytes: water chemical potential results and Gibbs–Duhem equation consistency tests, J. Chem. Phys., 2013, 139, 124505.
75. Z. Mester and A. Z. Panagiotopoulos, Mean ionic activity coefficients in aqueous NaCl solutions from molecular dynamics simulations, J. Chem. Phys., 2015, 142, 044507.
76. N. E. R. Zimmermann, B. Vorselaars, J. R. Espinosa, D. Quigley, W. R. Smith, E. Sanz, C. Vega and B. Peters, NaCl nucleation from brine in seeded simulations: sources of uncertainty in rate estimates, J. Chem. Phys., 2018, 148, 222838.
77. The PLUMED Consortium, Promoting transparency and reproducibility in enhanced molecular simulations, Nat. Methods, 2019, 16, 670–673.

---

Footnotes

† Electronic supplementary information (ESI) available: Additional figures. See DOI: https://doi.org/10.1039/d3fd00133d

‡ Current address: Department of Materials, Design and Manufacturing Engineering, The University of Liverpool, Liverpool, L69 3GH. Email: E-mail: aaron.finney@liverpool.ac.uk

---