Origin and control of ionic hydration patterns in nanopores

In order to permeate a nanopore, an ion must overcome a dehydration energy barrier caused by the redistribution of surrounding water molecules. The redistribution is inhomogeneous, anisotropic and strongly position-dependent, resulting in complex patterns that are routinely observed in molecular dynamics simulations. Here, we study the physical origin of these patterns and of how they can be predicted and controlled. We introduce an analytic model able to predict the patterns in a graphene nanopore in terms of experimentally accessible radial distribution functions, giving results that agree well with molecular dynamics simulations. The patterns are attributable to a complex interplay of ionic hydration shells with water layers adjacent to the graphene membrane and with the hydration cloud of the nanopore rim atoms, and we discuss ways of controlling them. Our findings pave the way to designing required transport properties into nanoionic devices by optimising the structure of the hydration patterns. The redistribution of water molecules when an ion passes through a nanopore is known to create complex patterns. Here, an analytical model accurately predicts the patterns when an ion passes through a graphene nanopore, and reveals the physical origins of the patterns.

I onic permeation through narrow water-filled channels and pores is of profound importance in both biophysics and technology 1,2 , finding applications in fuel cells 3 , biological ion channels 4 , water desalination 5-8 , gas 9 and isotope 10 separation, DNA sequencing 11,12 , and "blue energy" harvesting [12][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 nanopores 8,[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 research 8,[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 pore 18,22,29 , as has been confirmed by numerous molecular dynamics (MD) simulations 5,8,22,30,31 . The latter reveals that the patterns are highly inhomogeneous, anisotropic and spatially clustered near membranes 32 , inside nanotubes 33 , and in nanoslits 16,34 . This clustering determines the distribution of charge around the ion, the local dielectric permittivity 35,36 , and the strength of the ions' electrostatic interactions with their surroundings 37 . However, to understand, characterise, and predict the hydration patterns theoretically is a challenging problem 2 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 channels 24,38,39 , functionalised sub-nanopores 24,25,40,41 , nanotubes 42 and nanoslits 16,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 lattice 43 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.
For simplicity, we consider the carbons and the ion as being clamped. The spatial distribution of water oxygens is given by where H(⋅) is the all-atom Hamiltonian, comprising the water-ion U iw , water-lattice U ic and water-water U ww interactions; Z is the partition function, β = 1/k B T, k B is Boltzmann's constant and T is the absolute temperature. Integration over momenta p w n leads to a factor that subsequently cancels out due to normalisation. Here fr w n g, fr c m g, r i 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 r 1 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 as 47,p. 77 Wðrj r c 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 multiion 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 α-Hemolysin 48 and to describe DNA translocation through nanopores 49 . Importantly, the water PMF is related rigorously to the correlation function W (CF) via [44][45][46] W ðrÞ ¼ Àβ À1 ln gðrÞ; ð3Þ assuming that W ! 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 T 50 . The relations (2) and (3) allow one to rewrite the water spatial density ρ(r) in the form 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. 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.
This can clearly be seen near the graphene lattice in the lower panels of Fig. 2.
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 (g iw term) and the water layers surrounding the graphene membrane (∏g cw 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 approximation 50-52,p.186 . It has also been generalised to describe the ice-water interface 50 and solvated DNA 53 . 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 10 2 −10 4 speed-up in comparison with all-atom MD simulations 50 . 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-ray 54 or neutron 55 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 liquids 46,[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. MoS 2 , hBN, WS 2 , 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- 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.
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 ∏g iw 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 g iw 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 works 8,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,  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.
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.
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 radius 60 .
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 . Secondly, we propose adjustments of some intrinsic features: (e) the nanopore geometry, (f) rim charge, (g) atom type and (i) layer number (see Supplementary . 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.
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 CFs 50 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 hydration 62 and reorientation of the water dipoles in confinement which reduces 63,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 diffusivity 65 in biological channels, as well as the properties of the electric double layer 66,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) theory 68 . 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 research 35,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 eclipse 71 ), stretching, skewing, bending, and twisting of the lattice, the type of atom around the rim of the pore, and the nature of the 2D material 72 . This opens the way to multi-parameter optimisation 73 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 pathways [74][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 process 7 ); and secondly, energy harvesting, to maximise energy conversion by minimising Joule heating, thereby raising the output power density above the level of 5 W/m 2 needed to make harvesting economically viable 12,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 sequencing 11,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 isomers 78 , as illustrated in Fig. 5c, d, thus benefiting the Materials Genome Initiative 79,80 . Finally, by applying a travelling transverse (corrugated) wave one might be able to create a nanoscale ionic pump driving ions against the electrochemical gradient 81 (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 applications 82 .

Conclusions
We have proposed and validated an analytic method for revealing and quantifying the hydration patterns around an ion near a subnanopore, 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 VMD 83 to build the systems in question, and NAMD2 84 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 TIP3P 85 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) scheme 86 . 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 VMD 83 .

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Code availability
MATLAB codes for computing the ionic hydration patterns according to Eq. (4) and creating a movie analogous to Supplementary Movie 2, and transforming the crystal lattice are available as ZIP archives in Supplementary Data 1 and 2 (see Supplementary Note 7 for details). Other codes that support the results reported in this paper and other findings of this study are available from the corresponding author upon reasonable request.
Received: 15 March 2021; Accepted: 6 May 2021; Fig. 6 Water density patterns in a (5,10) CNT with a mechanical wave travelling along with it. Panels a-d show frames taken at intervals t m ¼ 0 6 ; 1 6 ; 2 6 ; 3 6 À Á 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.