Introduction

Ionic permeation through narrow water-filled channels and pores is of profound importance in both biophysics and technology1,2, finding applications in fuel cells3, biological ion channels4, water desalination5,6,7,8, gas9 and isotope10 separation, DNA sequencing11,12, and "blue energy” harvesting12,13,14. Angström-sized pores in artificially fabricated membranes are promising candidates for reproducible on-demand engineering and control of selective ionic and molecular transport.

It is well-known that the structure of ionic hydration shells is one of the key factors controlling ionic permeation in nanopores8,15,16,17. In confined geometries, some of the water molecules are lost from the hydration shells, and the resultant dehydration barrier is the subject of extensive research8,17,18,19,20,21,22,23,24,25,26,27,28. One outcome is an appreciation of the extraordinary complexity of the hydration patterns around an ion in the pore18,22,29, as has been confirmed by numerous molecular dynamics (MD) simulations5,8,22,30,31. The latter reveals that the patterns are highly inhomogeneous, anisotropic and spatially clustered near membranes32, inside nanotubes33, and in nanoslits16,34. This clustering determines the distribution of charge around the ion, the local dielectric permittivity35,36, and the strength of the ions’ electrostatic interactions with their surroundings37. However, to understand, characterise, and predict the hydration patterns theoretically is a challenging problem2 that hitherto has remained unresolved.

We now develop an analytic model that can predict the hydration patterns of ions in nano-confinement through the use of radial distribution functions (RDFs), evaluated in the bulk electrolyte. We show that the patterns result from interference between the ionic hydration shells, water layering near the membranes, and the hydration cloud around the nanopore. We compare the analytically predicted patterns with the results of MD simulations for different ionic locations near graphene nanopores with a diversity of shapes, geometries and strains. The results are applicable to narrow biological channels24,38,39, functionalised sub-nanopores24,25,40,41, nanotubes42 and nanoslits16,34. Thus, this work takes an important step towards controlling the function of nanoscale devices by tailoring the hydration patterns.

Results

Model

We consider graphene lattice43 with nanopore and an ion (e.g. K+) in a water box (see Fig. 1). The pore is created by removing one hexagon of carbon atoms from near the geometrical centre of the lattice.

Fig. 1: Model setup.
figure 1

A hexagonal pore in pristine graphene (grey balls and sticks) and a K+ ion (purple sphere) centred at Z = −0.4 nm. The colour plane represents the density pattern of water oxygen atoms in the plane at Z = −0.28 nm. (There will be comparable colour planes at other Z values.) The colours indicate density values relative to the bulk, brown being the highest (4) and dark blue being the lowest (0) value. The density has been evaluated by MD simulation, and it is compared with theory, Eq. (4), in Supplementary Fig. 1.

Full size image

For simplicity, we consider the carbons and the ion as being clamped. The spatial distribution of water oxygens is given by

$$\rho \left({{\bf{r}}}_{1}^{{\rm{w}}}| \left\{{{\bf{r}}}_{m}^{{\rm{c}}}\right\},{{\bf{r}}}^{{\rm{i}}}\right) \\ =\frac{1}{Z}\int \ldots \int {{\rm{e}}}^{-\beta H\left({{\bf{r}}}_{1}^{{\rm{w}}},\ldots ,{{\bf{r}}}_{{N}^{{\rm{w}}}}^{{\rm{w}}},{{\bf{p}}}_{1}^{{\rm{w}}},\ldots ,{{\bf{p}}}_{{N}^{{\rm{w}}}}^{{\rm{w}}}| \{{{\bf{r}}}_{m}^{{\rm{c}}}\},{{\bf{r}}}^{{\rm{i}}}\right)}\\ \times {\rm{d}}{{\bf{r}}}_{2}^{{\rm{w}}}\ldots {\rm{d}}{{\bf{r}}}_{{N}^{{\rm{w}}}}^{{\rm{w}}}\,{\rm{d}}{{\bf{p}}}_{1}^{{\rm{w}}}\ldots {\rm{d}}{{\bf{p}}}_{{N}^{{\rm{w}}}}^{{\rm{w}}}$$
(1)

where H(⋅) is the all-atom Hamiltonian, comprising the water–ion Uiw, water–lattice Uic and water–water Uww interactions; Z is the partition function, β = 1/kBT, kB is Boltzmann’s constant and T is the absolute temperature. Integration over momenta \({{\bf{p}}}_{n}^{{\rm{w}}}\) leads to a factor that subsequently cancels out due to normalisation. Here \(\{{{\bf{r}}}_{n}^{{\rm{w}}}\}\), \(\{{{\bf{r}}}_{m}^{{\rm{c}}}\}\), ri represent the 3D coordinates of the water molecules, atoms of the lattice, and the ion, respectively. Water molecules can be chosen arbitrarily, and so we omit the sub-index in r1 for clarity, so that r is used hereafter. Supplementary Discussions 1 and 2 in Supplementary Methods provide a more detailed derivation of the quantities in this section.

Direct application of the analytic formula (1) is not feasible due to a large number of degrees of freedom in the solution. However, statistical averaging over all water molecules generates rigorously the potential of the mean force (PMF)44,45,46 which can be approximated as47,p. 77

$$W({\bf{r}}| \left\{{{\bf{r}}}_{m}^{{\rm{c}}}\right\},{{\bf{r}}}^{{\rm{i}}})={{\mathcal{W}}}_{{\rm{iw}}}({\bf{r}}-{{\bf{r}}}^{{\rm{i}}})+\mathop{\sum }\limits_{m}^{{N}^{{\rm{c}}}}{{\mathcal{W}}}_{{\rm{cw}}}({\bf{r}}-{{\bf{r}}}_{{m}}^{{\rm{c}}}).$$
(2)

Note that the interactions of the ion or carbon atoms with other water molecules are already incorporated implicitly within the PMF. Note also that a similar strategy—construction of a multi-ion PMF from spherically symmetrical pairwise components measured in separate MD simulations— has been used successfully in an atomic-resolution Brownian dynamics study of the ionic selectivity of α-Hemolysin48 and to describe DNA translocation through nanopores49.

Importantly, the water PMF is related rigorously to the correlation function \({\mathcal{W}}\) (CF) via44,45,46

$${\mathcal{W}}\;\;({\bf{r}})\;=-{\beta }^{-1}{\mathrm{ln}}\,g({\bf{r}}),$$
(3)

assuming that \({\mathcal{W}}\to 0\) and g(r) → 1 as r → ∞. We evaluate the RDFs in bulk through separate MD simulations (see Discussions 1 and 2 of Supplementary Methods for details). Note that this needs to be done only once for a given density ρ0 and temperature T50.

The relations (2) and (3) allow one to rewrite the water spatial density ρ(r) in the form

$$\frac{\rho ({\bf{r}}| \left\{{{\bf{r}}}_{m}^{{\rm{c}}}\right\},{{\bf{r}}}^{{\rm{i}}})}{{\rho }_{0}}={g}_{{\rm{iw}}}\left(| {\bf{r}}-{{\bf{r}}}^{{\rm{i}}}| \right)\cdot \mathop{\prod }\limits_{m}^{{N}^{{\rm{c}}}}{g}_{{\rm{cw}}}\left(| {\bf{r}}-{{\bf{r}}}_{{m}}^{{\rm{c}}}| \right).$$
(4)

where ρ0 is the bulk water density (see Supplementary Discussion 1 for details). Note that the modulus sign reduces the 3D ion–water and carbon–water interactions to ion–oxygen and carbon–oxygen interactions, respectively. Thus, one is relying on the radial density functions (RDFs) g(r) which are one-dimensional functions of distance, in contrast to three-dimensional CFs.

In the derivation, we take no explicit account of water–water interactions, in particular between their oxygen atoms. This simplification leads to the water being represented as an ideal gas, and it overestimates the water density near the graphene walls. This can clearly be seen near the graphene lattice in the lower panels of Fig. 2.

Fig. 2: Comparison of analytic theory with MD simulations.
figure 2

The graphene sheet is shown in section by the vertical line of carbon atoms, in which the pore appears as a centrally placed gap with the horizontal axis of symmetry. The O− density in the vicinity of an ion (purple sphere) near the pore, derived from the MD simulations (a)–(c) is compared with the predictions of Eq. (4) (d)–(f). Panels a and d show the oxygen density around an empty pore; c and f show the density for a K+ ion in the centre of the pore; b and e show the density when the ion is centred 0.4 nm from the entrance to the pore. The fine structure of the oxygen density pattern around a K+ ion at the pore centre, obtained from MD simulations (g), (i) and from theory Eq. (4) (h), (j). The concentric black circles represent the first two maxima of the K+−O− RDF, and are given as guides to the eye. The graphene lattice is indicated in a grey ball-and-stick representation.

Full size image

Equation (4) allows one to compute the density of the water wetting a given set of atoms as a product of component RDFs, and to evaluate the water PMF via the Boltzmann inversion (3). Equation (4) also shows that the water pattern results from interference between the ionic hydration shells (giw term) and the water layers surrounding the graphene membrane (∏gcw term). Supplementary Discussion 3 illustrates the emergence of the hydration pattern when groups of carbon atoms are added sequentially to compose a nanopore. Equation (4) represents the well-known Kirkwood superposition approximation50–52,57,51,,p.186. It has also been generalised to describe the ice–water interface50 and solvated DNA53. Here, we apply this relation for the first time to describe the water density patterns around an ion near an artificial nanopore.

The usefulness of an approach based on Eq. (4) is twofold. First, the use of the bulk RDF values to describe water density near the pore provides for a 102−104 speed-up in comparison with all-atom MD simulations50. In turn, the single-component bulk RDFs g(r) can be obtained in a number of ways. The most straightforward is to use experimental values derived from the structure factor measured by X-ray54 or neutron55 diffraction. Thus one can avoid uncertainties related to the properties of the force field (parametrisation, polarisability) and thus connect the experimentally measurable structural properties of the bulk electrolyte solution with those near a nanopore. Other approaches rely on MD, Monte-Carlo or DFT simulations, or on the integral theory of liquids46,56,57,58,59. In the present work, MD-derived RDFs have been chosen in order to simplify the comparison between the MD simulations and analytics.

Secondly, the method can be readily extended to any material or pore geometry. The latter paves the way towards changing the pore parameters—size, shape, number of layers, strain tensor, rim charge, functional groups and choice of lattice type (e.g. MoS2, hBN, WS2, etc. as well as graphene)—with the purpose of changing the dehydration pattern. The pattern defines the dehydration energy profile and thus determines the pore’s permeability and selectivity, which are vital for inverse-designing the function of Angström-scale ionic devices. In Supplementary Discussion 4, we consider a tentative analytical way of connecting the hydration patterns and the single-ion PMF. Below we compare the ionic hydration patterns, obtained from MD simulations, with those from theory (Eq. (4)).

MD simulation results

Figure 2 compares the MD results (top row) with the corresponding theoretical predictions, Eq. (4) (middle row). They agree well in that the qualitative features of the water distributions are well reproduced: as described by the ∏giw term in Eq. (4), superposition of the contributions from all the carbon atoms results in a decaying plane density wave in the z-direction, matching the troughs and peaks of the C–OH RDFs (see Supplementary Fig. 2). The presence of the pore in an otherwise intact lattice distorts the water layers near the membrane and allows water molecules to exist in the pore as revealed by the strength of their distribution.

An ion in the pore (Fig. 2c), suggests spherically symmetric spatial density waves, as implied by the term giw in Eq. (4). The resultant superposition of the plane wave from the graphene with the spherically symmetrical hydration shells from the ion leads to a reshaping of the density such that the first hydration shell, while partially remaining, becomes significantly altered. In particular, it acquires maxima where the O–C and O–K+ RDF maxima coincide minima when both RDFs reach minima and mid-values where the maximum of one RDF overlaps with the minimum of another. Circular waves in Fig. 2b, c, e, f around the ion corresponding to the first two maxima of the intact bulk O–K+ RDF, indicated by the concentric circles at their respective radii. The density around the centrally located ion preserves the first hydration shell, but with a modulation due to the graphene atoms; the structure of the second shell is significantly affected by interactions with carbon atoms. Thus, the observed water patterns arise from interference between the ionic hydration shell and the hydration cloud around the nanopore. Analysis of the hydration patterns around Na+ and Cl− ions, shown in Supplementary Note 1, reveals similar features.

The central parts of Fig. 2b, e are shown at higher resolution in panels (g) and (h). One can see the fine structure in the water density pattern resulting from the superposition of the ion–oxygen and carbon–oxygen interactions. The density crescent at Z ≈ 0.6 nm in Fig. 2g, h is induced by the first maximum of the RDF, as indicated by the concentric circles, while the two peaks appear due to the oxygen interaction with the carbon atoms. The low-density area around [0.45, ± 0.25] nm emerges due to the minimum of the K+–O− RDF located at around 0.35 nm. The “island” around [0.175, 0] nm indicates the position of the trapped water molecule, and highlights the predictive power of the proposed analytical method (see Supplementary Note 2 for more details). The cross-sections of the crescents in Fig. 2g, h, made at 0.25 nm, are shown in panels (i) and (j). These arise from an interplay between the oxygen–ion and oxygen–pore interactions: the minimum at the origin arises due to the convex shape of the first hydration shell near the ion, while the hexagonal structure is inherited from the hexagonal geometry of the graphene pore (grey balls and sticks). Similar water density patterns in the plane of the nanopores have been observed in numerous MD simulation works8,17,18,22,30 recently.

Figure 3 combines the previous results by illustrating the 3D isosurfaces of the hydration shells around the ion and the pore, taken at a level of 1.15. Panels a and b show that the results from MD and theoretical prediction Eq. (4), respectively, agree well. These figures reproduce the layering of water near the graphene sheet, with the first three layers visible (see Supplementary Note 3 for details), as well as the spherical hydration shells around the K+ ion. Parts of the preserved 1st and 2nd hydration shells are clearly evident. The figure emphasises the inhomogeneity and anisotropy of the water patterns emerging from the combined effects of the ion and the graphene lattice. The origin of multiple water layers is illustrated in Supplementary Fig. 12.

Fig. 3: Three-dimensional ionic hydration patterns.
figure 3

Water distribution around a K+ ion (purple sphere) at Z = 0.4 nm near the pore in a graphene lattice (black ball-and-stick representation), obtained from a MD simulations and b theory (Eq. (4)). An isovalue of 1.15 was chosen to reveal the layered structure of the hydration around the ion and the graphene lattice. Both panels have been smoothed with a 5-point window.

Full size image

We have also studied how the hydration structure evolves as the ion’s position changes. The results are summarised in Supplementary Note 4 and the two Supplementary frame series (Supplementary Movies 1 and  2). These sets clearly illustrate that the water structure is significantly modified, even when the ion is located more than 1 nm away from the nanopore. In other words, when the lattice-induced (hydration layers) and ion-induced (hydration shells) inhomogeneities begin to overlap, the ion–pore interaction starts to change. This result supports the notion that solvent-mediated interactions extend the range of the interactions with the pore by a large factor compared with the bare ion radius60.

Finally, in Fig. 4 we illustrate the predicted effects of the set of tools that we propose for controlling the dehydration pattern. They can be categorised into two groups. The first set is extrinsic and includes (a) stretching, (b) skewing or asymmetrical stretching, (c) bending, and (d) twisting of the lattice (see Supplementary Notes 8–11). Secondly, we propose adjustments of some intrinsic features: (e) the nanopore geometry, (f) rim charge, (g) atom type and (i) layer number (see Supplementary Notes 12–17). Our theory yields good agreement with the published crown-like water density patterns (panel (h)) from ref. 61. More examples of these tools in action are given in Section 3 of the Supplementary Information.

Fig. 4: Tools to control the hydration pattern.
figure 4

These include the extrinsic processes of: a stretching, b skewing or asymmetric stretching, c bending, d twisting, and the intrinsic features of—e choice of pore isomer, f charging of the rim atoms (red balls), g functionalising the rim with oxygen atoms (black balls), h use of a multi-layer pore, as exemplified in a molecular dynamics simulation (h adapted from ref. 61, copyright 2019, with permission from Elsevier), and i theoretical prediction, Eq. (4). Density slices have been taken at Z = −0.23 nm (a)–(d) and Z = −0.28 nm (e)–(g). The graphene lattice is indicated in a grey ball-and-stick representation.

Full size image

It is evident that Eq. (4) can reproduce well the qualitative features of the water distribution near the pore given in Fig. 2. We point out, however, that there is a numerical discrepancy—a factor of ~3—between the peak values of the maps, clearly seen between the MD-generated and theoretical density of the “ridges” near the graphene lattice. We attribute this discrepancy to our omission of the water–water interactions in Eq. (4), which would imply an absence of short-range repulsion between oxygen or hydrogen atoms of separate water molecules, so that arbitrarily small separations can be realised. The consequence is denser water distribution, localised near the carbon atoms. Also, the reduction from a full 3D model to the radial CF picture overlooks the water–water relative orientations and thus limits the applicability of the method in carbon nanotubes. We are currently working on eliminating these shortcomings of the technique by using higher-dimensional CFs50 and applying machine learning (see Supplementary Discussion 5 for a toy example).

Discussion

Our approach can readily be extended to account for the density of hydrogen atoms, and thus to evaluate the average charge density in the domain. Doing so promises deep insight into water-mediated interactions, including the pore rim’s hydration62 and reorientation of the water dipoles in confinement which reduces63,64 the dielectric permittivity. The latter is needed for quantifying the ion and polymer interactions with the pore and with the external electric field and for understanding the ionic diffusivity65 in biological channels, as well as the properties of the electric double layer66,67. The method can also be used to define the demarcation interface between the pore and the protein in continuous theories of ion transport such as the Poisson–Nernst–Planck (PNP) theory68. The method can thus provide fundamental insights into the electrostatics and dynamics of ions in sub-nanopores, where the dielectric properties on the nanoscale remain a challenging unsolved problem that is subject to intensive research35,36,69,70.

Our results demonstrate the applicability of the method to describing water density patterns under quite general conditions, e.g. for arbitrary ionic positions, pore geometry (size, shape, number of layers, offset eclipse71), stretching, skewing, bending, and twisting of the lattice, the type of atom around the rim of the pore, and the nature of the 2D material72. This opens the way to multi-parameter optimisation73 of the energy landscape by adjustment of the corresponding hydration patterns. For example, by iteratively choosing the most appropriate combination of pore geometric parameters, and subsequently manipulating the pore shape externally (e.g. by stretching), the nanopore’s permeability and selectivity can be optimised and regulated, providing for permeation along specified pathways74,75,76. This approach could prove useful in addressing two highly topical technological challenges: first, water desalination, to maximise the quality and energy efficiency while minimising the environmental impact of the process7); and secondly, energy harvesting, to maximise energy conversion by minimising Joule heating, thereby raising the output power density above the level of 5 W/m2 needed to make harvesting economically viable12,14.

The same design and optimisation strategy should serve to slow down the translocation and guide the orientation of DNA (see Fig. 5a for a possible example), which constitutes one of the current challenges to fast genome sequencing11,77. Likewise, it is believed that the accuracy of nucleobase detection using the Coulter principle can benefit from an optimisation of the features of the nanopore, Fig. 5b (see also Supplementary Movies 3, 4, and Note 5). Our approach would allow one to catalogue solvated 2D and 3D pore isomers78, as illustrated in Fig. 5c, d, thus benefiting the Materials Genome Initiative79,80. Finally, by applying a travelling transverse (corrugated) wave one might be able to create a nanoscale ionic pump driving ions against the electrochemical gradient81 (see Fig. 6, Supplementary Movie 5 and Note 6). Thus, the proposed method paves the way to the optimisation and design of novel nanoionic devices in many contexts and applications82.

Fig. 5: Possible applications of our proposed method to DNA sequencing and the cataloguing of the solvated nanopore isomers.
figure 5

a Hydration pattern wetting the single-stranded DNA in a nanopore. Densities are obtained from Eq. (4) and are shown as grey isosurfaces (isovalue = 1.2) and coloured cross-sections. The DNA is represented by connected coloured spheres (oxygen (red), nitrogen (green), carbon (black), phosphorus (yellow)). Graphene lattice is shown in grey balls and sticks. b Three consecutive snapshots of the K+ ionic trajectory along the groove of the DNA, depicted by the red line. c, d Represent the nontrivial structure of the solvation patterns around one of the pore isomers, taken from ref. 78 (isomer for the case of N = 19, the frequency of observation of the isomer—4%—total number of possible isomers—48878). The isosurfaces are calculated at the isovalue of 1.2. Graphene lattice is shown in black balls and sticks.

Full size image
Fig. 6: Water density patterns in a (5,10) CNT with a mechanical wave travelling along with it.
figure 6

Panels a–d show frames taken at intervals \({t}_{m}=\left(\frac{0}{6},\frac{1}{6},\frac{2}{6},\frac{3}{6}\right)T\), where T is the period of the wave. This provides an illustration of how a transverse wave, travelling along the CNT axis, may lead to the directed motion of ions against the electrochemical gradient, thus providing a potential design for mechanically induced ion-pumping. See Supplementary Movie 5 for a visualisation.

Full size image

Conclusions

We have proposed and validated an analytic method for revealing and quantifying the hydration patterns around an ion near a sub-nanopore, using a product of the RDFs measured in the bulk electrolyte. The complexity and fragmentation of these patterns result from interference between ionic hydration shells and the hydration cloud around the nanopore. The method can be extended to account for the distribution of hydrogen atoms, and hence to obtain insight into local electrostatic interactions on the nanoscale. It is fast, and it allows one to take explicit account of the nanopore’s detailed geometry (size, length, shape, number of layers, offset eclipse), pore constituents (charge and atom type), and of the ion’s position, as well as implicitly to take account of the ion type, solvent, temperature, and pressure. Multi-parameter optimisation of the pore parameters and surroundings is expected to create a way of engineering and regulating the energy barrier in nanopores via control of the hydration pattern. This strategy opens a new avenue to designing, optimising, and controlling the permeability and selectivity properties of nanoionic devices, as well as regulating the translocation of DNA through solid-state nanopore sequencers.

Methods

MD simulation details

We used VMD83 to build the systems in question, and NAMD284 with the CHARMM27 force field for MD simulations at T = 300 K with a time step of 1 fs and the velocity Verlet algorithm. The Lorentz–Berthelot combining rule was applied. The TIP3P85 water model was used. A cutoff of 1.2 nm for non-bonded interactions (Lennard–Jones and Coulomb) with a switching distance of 10.0 Å (1 nm) was used, and the full electrostatic calculation was performed using the particle mesh Ewald (PME) scheme86. Periodic boundary conditions were imposed in all directions.

All lattice atoms were fixed by setting the beta parameter to 1. The system first underwent equilibration during the initial 1000 steps in the Nosé–Hoover thermostat at pressure p = 1 atm (101.325 kPa), with the remaining simulation running under NVT conditions. The frames were captured every 100 steps. No external electric field was applied. Unless otherwise stated, typical production runs took 10 ns and were sampled every 100 steps to yield statistically significant figures.

RDFs were measured for a free atom (C, K+, Na+) in the water surroundings. When an ion’s RDF was measured, the corresponding number of free counterions (Cl− for K+, and K+ for Cl−) was added to neutralise the system. Finally, RDFs were measured using the VMD plugin gofr. Molecular structures were visualised and rendered using VMD83.

Data analysis

MD trajectories were read by a readdcd.m and readpdb.m scripts MDToolbox87 for MATLAB88. A homemade script to build and analyse the histrogram distributions included built-in MATLAB functions and some external functions: histcn.m89, smooth2a.m90. All densities were normalised by the bulk values. A 3 × 3-point smoothing window was applied to suppress the pixelation of the MD histograms, unless otherwise stated.