Instantaneous , parameter-free methods to define a solute ’ s hydration shell

A range of methods are presented to calculate a solute's hydration shell from computer simulations of dilute solutions of monatomic ions and noble gas atoms. The methods are designed to be parameter-free and instantaneous so as to make them more general, accurate, and consequently applicable to disordered systems. One method is a modified nearest-neighbor method, another considers solute-water Lennard-Jones overlap followed by hydrogen-bond rearrangement, while three methods compare various combinations of water-solute and water-water forces. The methods are tested on a series of monatomic ions and solutes and compared with the values from cutoffs in the radial distribution function, the nearest-neighbor distribution functions, and the strongest-acceptor hydrogen bond definition for anions. The Lennard-Jones overlap method and one of the force-comparison methods are found to give a hydration shell for cations which is in reasonable agreement with that using a cutoff in the radial distribution function. Further modifications would be required, though, to make them capture the neighboring water molecules of noble-gas solutes if these weakly interacting molecules are considered to constitute the hydration shell.


I. INTRODUCTION
2][3] The number of molecules in the hydration shell, termed as the coordination number n, is regarded as the number of water molecules adjacent to or perturbed by the solute.Despite the convenience and usefulness of the hydration shell, it has no rigorous definition in a liquid or solution.Many methods exist to define it which inevitably lead to a wide range of values of n and thus ambiguity in understanding solution structure.In a crystal, there is clarity about what constitutes a solute's hydration shell because molecules vibrate about discrete, well-defined, symmetrical positions, making the nearest neighbors (NNs), all approximately equidistant, obvious candidates for the hydration shell.In a liquid or solution, however, the molecules are continuously and non-uniformly distributed, making the demarcation of nearest neighbors unclear, 4 hence leading to many methods and values.
The main approach to determine n is based on the idea of proximity to the solute.It most commonly uses the radial distribution function g(r), available from X-ray diffraction, neutron diffraction, or computer simulation, an attribute that makes g(r)-based methods a powerful test of simulation models.However, there is no unique way to extract n from g(r).5][6][7][8][9] Alternatively, slightly smaller values are found by integrating over the symmetrized, first peak of rg(r) or r 2 g(r). 1,4,10,11In recent times, placing the cutoff at the first minimum in g(r) has become the predominant method.A number of variants are possible by computer simulation.With pre-computed cutoffs, one can determine values of n for each configuration and build up a histogram of the ensemble.Another approach considers the g(r)s of individual, ranked nearest neighbors to the solute and determines an integer value of n at the rank at which the g(r) moves from the first to the second minimum, so avoiding the need for a distance cutoff.A related method considers a cutoff in the force-field energy between two molecules, which depends on both distance and atom type 12,13 but this is rarely used for hydration shells.While conceptually simple, cut-off based methods suffer from a sensitivity of n to the cutoff, require well-converged g(r)s to produce a clear first minimum, have limited ability to resolve overlapping hydration shells arising from the variety of atomic arrangements in disordered systems, neglect correlations between molecules, and there is uncertainty if cutoffs are transferable between different atom types, thermodynamic conditions, and phases, versus the otherwise combinatorial number of cutoffs required.
Another class of methods infers n from the average perturbation from bulk behaviour of some property probed spectroscopically or thermodynamically such as dynamics, activity, density, or compressibility. 1 By these methods, n is often referred to as an effective solvation number 14 rather than coordination number.These approaches work better when solute-bound water has distinctly different properties to bulk water but are less useful when the perturbation is small or when there are multiple solutes contributing to the perturbation.
Implicit is the assumption that water can be partitioned into perturbed and bulk water, 15 a division often derived from the concentration or temperature dependence of the measured perturbation. 14,16However, there is growing evidence that all water molecules are perturbed by a solute, [17][18][19][20][21][22][23][24] less so at longer range than those near the solute, but never such that bulk behaviour is restored beyond some finite distance from the solute.Perturbation-based methods often yield quite different values of n compared to those based on proximity, the most well-known discrepancy being the n values of anions being close to zero when measured by nuclear magnetic resonance (NMR) proton shifts 25 or dielectric relaxation. 14 different strategy is to define the hydration shell in an instantaneous, parameter-free fashion that reflects each molecule's particular environment, accounts for correlation between particles, potentially gives greater understanding of why the hydration shell forms, and makes possible a general method applicable to a wide range of conditions.Only a few such methods exist.8][29][30][31][32][33][34][35] Atoms are defined as neighbors if their polyhedra share a common surface.While attractive in principle, the algorithm is computationally expensive, very sensitive to particle motion, and yields values of n much larger than those using g(r) cutoffs, being 10-20 for common liquids. 27,28,35,36 recent method is the solid-angle nearest-neighbor (SANN) algorithm 37 which sums up the solid angles of all nearest neighbors until the next-nearest particle lies beyond a cut-off distance at which the total solid angle equals 4π.Effectively, the cutoff is where a particle's distance exceeds the offsetted average distance of all closer particles, making the method similar to the nearest-neighbor method.
For the more specialized case of hydrogen-bond coordination, there is a parameter-free hydrogen-bond definition 38,39 whereby a hydrogen donates to the nearest acceptor as long as that is the shortest donor-acceptor distance between the two molecules involved.If there is a shorter distance, then the donor hydrogen is said to have a broken hydrogen bond.A more general form of this definition considers electrostatic force in place of distance. 39This takes into account both distance and charge and is useful if there are multiple types of acceptor.These two definitions had previously been referred to as "topological" because they are based on the idea of resolving hydrogen-bond arrangements using the separating transition states at which a donor experiences opposing, equal forces from two acceptors.However, the names "nearest acceptor" and "strongest acceptor" are adopted here because they are more informative and consistent with the naming scheme for the other definitions to be introduced later in this paper.The strongest-acceptor definition is readily applied to determine n for acceptors such as anions in water because water molecules donate a hydrogen to the anions.Effectively a nearest-neighbor method, it works well for hydrogen because a hydrogen usually donates to only one acceptor.However, it can neither detect hydrogen bonds with a hydrogen bifurcated between two oxygens 40,41 nor be applied to cations or neutral atoms with which water molecules do not form hydrogen bonds.Considering only a cation's nearest acceptor would lead to the extremely small value of n = 1 for a cation and n = 0 for a neutral atom.
In this paper, we explore five new methods to determine n from computer simulations of aqueous cations and noble gases.The methods are instantaneous, parameter-free and incorporate both proximity to the solute and perturbation by the solute.One method assesses whether ranked nearestneighbor water molecules are closer to the solute than another water.Another method considers a water molecule as lying in the hydration shell from the moment it approaches closer to the solute-water Lennard-Jones minimum until its acceptor hydrogen bonds change.Three related methods compare water-solute forces with various combinations of water-water forces to assess if each water interacts more with the solute than with its neighboring water molecules.We compare the n values by these methods with those derived by g(r) cutoffs, nearest neighbors, and the strongest-acceptor hydrogen-bond method for anions.

A. Hydration-shell methods
The eight methods used to calculate n, the number of water molecules in a solute's hydration shell, are explained here, supported by matching illustrations in Figure 1.Each method is summarized in a brief one-sentence explanation, followed by the details of how n is calculated.A qualification may be included about whether the method is applicable to cations, anions, or noble-gas solutes.Hydrogen bonds are determined using the strongest-acceptor definition. 38,39 g(r) cutoff (GC): All water molecules within a distancecutoff are in the solute's hydration shell.The radial distribution function g(r XO ) is constructed from all trajectory frames, where r XO is the solute-oxygen distance.Integrating 4πr 2 g(r XO ) from zero to the first minimum in g(r XO ) and multiplying by the average number density yields n. 2. Nearest Neighbors (NNs): All nearest-neighbor water molecules up to an integer cutoff lie in the solute's hydration shell.For every snapshot, water molecules are ranked as nearest neighbors of the solute.A set of potentials of mean force (PMFs) PMF i = −RT ln g i (r) is constructed for each water with rank i from all simulation frames, where R is the universal gas constant, T is temperature, and g i (r) is the radial distribution function for water of rank i.For cations and noble gases, r is r XO , and for anions, r is r XH .The crossover point in the PMFs from small to large distance is defined at the point between the two highest, consecutive PMFs, PMF i , and PMF i+1 .In other words, the pair with the highest value of the PMF at the point of intersection PMF i = PMF i+1 defines the crossover point for being in the solute's hydration shell.The ranking i of the lower PMF in the pair becomes the integer cutoff for the hydration shell, and this i is set to n.This method is closely equivalent to the GC method because the crossover point corresponds to a point of low probability at the first minimum.3. Nearest Neighbor's Nearest Water (NNNW): This also assigns to the hydration shell all nearest-neighbor water molecules up to an integer cutoff.However, the PMFs considered are designed to test whether a water molecule interacts more strongly with the solute, in which case it is in, or with a nearby water molecule, in which case it is out.For cations and noble gases, PMF i uses the distance r OH between the water i's oxygen and the nearest hydrogen of another water.For anions, PMF i uses the distance r HO between the water i's hydrogen nearest the anion and the nearest oxygen of another water.For PMFs that are stabler at a large distance, this indicates that the water must be interacting strongly with the solute and thus in the hydration shell.If the PMF is stabler at a small distance, this indicates a water-water hydrogen bond, that the water interacts weakly with the solute, and that it is not in the hydration shell. A for the NN method, the integer cutoff, which defines n, is between the consecutive nearest neighbors whose PMF intersection is the highest.

Strongest-Acceptor Hydrogen Bond (SAHB, anions only):
The solute's hydration shell comprises all water molecules forming a SAHB to the anion.A water lies in the anion's hydration shell if at least one of its hydrogens has the anion as the acceptor with which it has the most attractive force.Summing all such waters per trajectory frame and averaging over all frames gives n.

Lennard-Jones Overlap/Strongest-Acceptor Hydrogen
Bond (LJO/SAHB, cations and noble gas solutes only): This assumes that a water remains in the solute's hydration shell from the moment it collides with the solute until there is a change in the hydrogen bonds it accepts.A water is in the hydration shell whenever it lies closer than the minimum of the solute-oxygen Lennard-Jones energy function.It remains in the shell for all subsequent structures in time until there is a change in the hydrogen bonds that this water accepts.Summing all such waters per trajectory frame and averaging over all frames gives n.
Anions are not considered because the SAHB definition is already adequate for them, and because water's hydrogens for the SPC/E model used here have zero Lennard-Jones radius.6. Weakest Donor (WD, cations and noble gases only): This puts a water in the hydration shell if the water's solutewater interaction is stronger than its weakest water-water hydrogen-bond interaction.A water lies in the hydration shell if the force between the solute and its oxygen, excluding Lennard-Jones repulsion, is more attractive than the force between the water's oxygen and the hydrogen donating most weakly to that oxygen.Summing all such waters per trajectory frame and averaging over all frames gives n. 7. Strongest Non-Donor (SND, cations and noble gases only): This is more relaxed than the WD definition and puts a water in the hydration shell if the water's solutewater interaction is stronger than its strongest water-water interaction with a water not interacting via a hydrogen bond.A water lies in the hydration shell if the force between the solute and the water's oxygen, excluding Lennard-Jones repulsion, is more attractive than the force between the water's oxygen and the hydrogen with the most attractive force to that water but not donating to it.Summing all such waters per trajectory frame and averaging over all frames gives n. 8. Weakest Donor/Strongest Non-Donor (WD/SND, cations and noble gases only): This is intermediate between the WD and SND definitions and examines if a water's solutewater interaction is weaker than its weakest water-water hydrogen-bond interaction but stronger than the strongest water-water non-hydrogen-bonded interaction.A water lies in the hydration shell if the force between the solute and the water's oxygen, excluding Lennard-Jones repulsion, is closer to the force between the water's oxygen and the hydrogen donating most weakly to that water than the corresponding force between the water's oxygen and the hydrogen with the most attractive force to that water but not donating to it.Summing all such waters per trajectory frame and averaging over all frames gives n.

B. Simulation protocol
All simulations were performed using the GROMACS 4.5.5 molecular dynamics package.The Lennard-Jones nonbonded parameters for Li + , Na + , K + , Rb + , Cs + , F − , Cl − , Br − , and I − are those of Joung and Cheatham. 42For the noble gas solutes Ne, Ar, Kr, and Xe, the parameters are those of Guillot and Guissani. 43Each system consists of a single ion or atom in a cubic box of 1000 SPC/E water molecules with length ∼3.1-3.3 nm.GROMACS automatically neutralizes any net charge of the system by providing a net background charge through dipole and charge corrections.Initial coordinate generation and elimination of short range repulsive interactions were carried out using Packmol. 44This was followed by NPT equilibration for 2 ns using a Berendsen barostat at 1 bar and velocity rescaling thermostat at 298 K. Data were then collected at 1 ps intervals for a further simulation 1 ns long.All simulations use LINCS to constrain bonds, a 2 fs time step, a 1 nm non-bonded cutoff, and particle mesh Ewald.All methods to calculate hydration shells from the GROMACS coordinate were carried out using our own code in Fortran 77 which can be found at http:// personalpages.manchester.ac.uk/staff/henchman/software/.

III. RESULTS AND DISCUSSION
The coordination numbers n calculated by all eight methods for all thirteen aqueous solutions are presented in Table I.The cutoffs used here to calculate the GC values, ranked from smallest solute to largest, are 0.27, 0.32, 0.36, 0.37, and 0.39 nm for the cations, 0.33, 0.38, 0.39, and 0.41 nm for the anions, and 0.49, 0.53, 0.54, and 0.55 nm for the noble gases.0][51][52] Some spread arises due to different simulation force field parameters, experimental accuracy, concentration dependence, and sensitivity to the cutoff.][55][56][57][58][59][60][61] The GC and NN values calculated here compare well with those from the literature.The purpose of this study, however, is to compare n from the new methods with GC values.
The PMFs required for the NN method, up to the fourteenth nearest water for ions and twenty-fifth nearest water for noble gases, are plotted in Figure 2 for all solutes except Cs + whose plots are in Figure S1 of the supplementary material, 62 being similar to those of Rb + .][65][66] The n values also closely mirror the GC values.This is expected, given that the switchover from a water being in the hydration shell to outside occurs near the first minimum in g(r).They are not exactly the same, though, especially for the noble gases, presumably because the additional r 2 weighting in g(r) brings about a slightly different position in g(r)'s minimum compared to the PMFs.The advantages of the NN method compared to GC are that there is no cut-off distance and the hydration shell's boundary is smoother.However, the method is still a post-processing technique rather than instantaneous because the PMFs must be compiled over the full structural ensemble.There are cases when a PMF straddles both hydration shells.For example, the fifth water of Li + and the seventh water of F − sometimes intrude from outside into the hydration shell while Na + 's sixth water is not always within the hydration shell.A more refined, non-integer value of n could be derived that includes the relative contributions in the two hydration shells.The switch-over from Na + to K + , noted by Collins, 67 is clear in that Na + 's PMFs display a clear gap, indicating a more defined hydration shell, whereas those of K + are spaced more uniformly.A smaller switch-over 67 is seen for anions between F − and Cl − .
The PMFs required for the NNNW method are plotted in Figure 3 for all solutes except Cs + in Figure S2. 62In the case of cations, the minima in Figures 3(a)-3(d) at r OH ∼ 1.8 Å represent a water molecule that donates to the water under consideration and therefore is assumed to exclude that water from the cation's hydration shell.The minima on the right indicate that no water molecule donates to that water because its oxygen is occluded by the cation.For the anions, the left minima indicate a water molecule to which the water under consideration donates, preventing it from donating to the anion and thus excluding it from the anion's hydration shell, while the right minima are water molecules that do not donate to another water and must therefore be donating  FIG. 3. Nearest-neighbor's nearest-water PMFs ranked as a function of water-solute distance r .For cations and noble gases, r = r OH is the distance between the ranked water's oxygen and its nearest hydrogen on another water.For anions, r = r HO is the distance between the ranked water's hydrogen nearest the anion and the nearest oxygen on another water.Plots are coloured as in Figure 2.
the hydrogen-bond forming minima on the left are usually stronger.Evidently, the NNNW method gives values of n based more on perturbation than proximity.While appropriate for the anions and strongly interacting cations, the NNNW method is clearly overly strict, and while parameter-free, it is not instantaneous.The SAHB method, only applicable to anions, gives n at around 6, peaking at 6.4 for Cl − and slightly less for the other anions.Its n values agree well with those from other simulations that use either the SAHB method but a different ion force field, 39 for which F − to I − are 6.1, 6.9, 7.5, and 7.6, or that use the GC method based on g(r XH ) 60,68 whose n values are 6.2-6.5 for F − and 7.0-7.4for Cl − .GC n values obtained from g(r XH ) by neutron diffraction 8,45,54,69 are also close at 6.7 for F − , 5.5-6.4 for Cl − , 5.5-5.9 for Br − , and 5.8 for I − .Evidently, n values derived from r XH are slightly lower than the values derived from r XO , especially for the larger anions.Ions of larger size have more neighbors, and yet their lower charge density weakens their ability to accept hydrogen bonds.
The LJO/SAHB n values give a remarkably close match to the GC n values for the cations, rising from 3.8 for Li + to 7.8 for Cs + .Thus, the method's workings provide an informative, although not unexpected, understanding of what constitutes the hydration shell.All water molecules hydrating cations contact the cation and only leave when other water molecules donate to them to replace the cation.The noble gas n values are improved compared to the NN and NNNW methods but are still much smaller than GC n values, being 8.1 for Ne and ∼6 for Ar, Kr, and Xe.Thus, it would appear that the LJO/SAHB method reasonably captures the inner subshell of water molecules hydrating non-polar solutes 70 but fails to capture the outer subshell because the water molecules here do not contact the solute.One could either regard the hydration shell as only the inner subshell or the cutoff used could be expanded beyond the Lennard-Jones radius, for example, to the point of inflection.It should be noted that the LJO/SAHB method may miss Lennard-Jones points of contact if a water molecule is adjacent to the solute but happens to be non-overlapping with the solute in every trajectory frame considered.The larger the time interval and fewer trajectory frames there are, the more likely that a point of contact would be missed.This could lead to a water not being in the hydration shell when it would have been had frames been saved more often.This may be responsible for the slightly lower n values by LJO/SAHB than GC.These results demonstrate that hydration shells can be defined instantaneously and without additional cut-off parameters.However, the Lennard-Jones radius could be regarded as a cut-off parameter, albeit a force-field proscribed one.
The WD method yields n values that increase from 4.1 for Li + to 5.5 for Na + but then decrease with increasing cation size to only 1.6 for Cs + .This mirrors the results from the NNNW method.The n values lie fairly close to those derived by the degree of perturbation either by temperature-dependent NMR proton shifts of 3-3.5 for cations 1,25 or by dielectric spectroscopy of 4.2 for Na + . 14Given the weak water-noble gas interactions compared to the water-water interactions, whether with the water's oxygen as presented here or with water's hydrogen, their n are all zero.Comparing the solutewater interaction to the weakest donor's interaction with the candidate water is likely too strict.The solute may be adjacent to the water but its interaction with it happens to be weaker than that of all the donors to that water.This possibility is addressed by the SND method discussed next.
SND n values vary weakly with cation size, ranging only from 8 to 9. Clearly, the SND method has a more generous criterion than the WD method for a water to be in the hydration shell.The n values agree well with those from the first four methods for the larger cations but they are too large for the smaller ions for which more distant water molecules are included.While the hydration shell may be observed to extend further away from the solute for ions of higher charge density, 3,71 the SND method appears to be not strict enough.The threshold interaction of water molecules not donating to the candidate water may be quite weak.The more generous criterion makes no difference to the noble gases, though, whose n are still all zero.
The WD/SND method is intermediate between the WD and SND methods because the threshold force is the average of the weakest-donor and strongest non-donor forces.The n values reflect this.There is a rise in n from 4.4 for Li + to 6.1 for K + and only slightly down again to 5.2 for Cs + .This improves the agreement with the GC and NN methods, especially for the smaller ions, but still leaves a discrepancy of 2-2.5 for the larger cations with respect to the GC method.Having a force cutoff closer to the strongest non-donor force could boost n of the larger cations, while still keeping similar n of the smaller cations, but at the cost of introducing a cut-off parameter.
Overall, the mixed results from comparing all methods make clear the difficulty in creating an instantaneous, parameter-free method to determine the hydration shell of a cation or noble gas solute in an analogous way to the successful force-based SAHB method for anions.Much of this can be attributed to the asymmetry of the water molecule and the preference of a hydrogen to donate to only one acceptor versus the oxygen's preference to accept more donors.The LJO/SAHB method yields n in reasonable agreement with GC values but with a weak dependence on the frequency of frame sampling and being low for noble gases.The three methods based on a comparison of solutewater and water-water forces, WD, SND, and WD/SND, yield a range of n, with the intermediate WD/SND having reasonable agreement with GC values.This agreement could be improved by shifting the force cutoff by an arbitrary amount but at the cost of introducing an arbitrary parameter.The n values of zero for the noble gas solutes could be boosted by only comparing the solute-oxygen and oxygen-oxygen Lennard-Jones forces instead of the combined Lennard-Jones and electrostatic forces, although this would fail to account for the electrostatics of weakly charged non-polar atoms.Alternatively, because the WD, SND, WD/SND, and SAHB methods depend more on perturbation than proximity, one would not have to assume that the GC method, being based on proximity, is the absolute reference for what constitutes a hydration shell.While the WD/SND and SAHB methods share the feature of comparing solute-water and water-water forces, ultimately, it is unsatisfactory to have different methods for solutes of different charges.The desirable method to pursue is one that is not only instantaneous and parameter-free but also applicable to any solute regardless of the charge or atom type.

IV. CONCLUSIONS
A series of methods have been presented to define the hydration shell of monatomic ions and solutes from computer simulation based on distances and forces.In addition to the existing g(r) cutoff, nearest-neighbor, and strongest-acceptor hydrogen-bond methods, these include a nearest neighbor's nearest water method, a method based on penetration of the Lennard-Jones minimum combined with rearrangement of strongest-acceptor hydrogen bonds, and three methods based on a comparison of solute-water and water-water forces.The varied results reinforce the existing picture of hydration that the hydration shell is ambiguous and can be defined in many ways.The most suitable method depends on what can be measured or calculated and how the hydration shells are to be used.Our focus has been on deriving new methods that give coordination numbers n close to GC values but are parameter-free and instantaneous so as to enable the necessary generality and accuracy for disordered systems.No method is found to meet all these specifications for all solutes.The NN method is parameter-free and well reproduces GC values but is not instantaneous.The WD/SND and LJO/SAHB methods for cations and SAHB method for anions meet the requirements of being parameter-free, instantaneous, and agreeing well with GC n values for ions.However, this is not the case for noble gases whose weak interactions with water bring about n = 0 for the WD/SND method and n about half the GC values for the LJO/SAHB method.It is unclear if GC values are the most appropriate benchmark, given their shortcomings such as sensitivity and transferability.We are undertaking more studies to assess this as well as consider alternative parameter-free approaches.

234501- 5 Chatterjee,FIG. 2 .
FIG. 2. Nearest neighbor PMFs ranked as a function of water-solute distance r .For cations and noble gases, r = r XO and for anions, r = r XH .The plots are shown up to the fourteenth nearest water for ions and twenty-fifth nearest water for noble gases.Even-numbered lines for cations are blue (a) Li + , (b) Na + , (c) K + , and (d) Rb + ; for anions are red (e) F − , (f) Cl − , (g) Br − , and (h) I − ; and for noble gas solutes are green (i) Ne, (j) Ar, (k) Kr, and (l) Xe.Odd-numbered lines are coloured black for the hydration shell, fading to light gray outside it, to aid legibility.