Ion clusters and ion-water potentials in MD-simulations

A well known, if little documented, problem in many molecular simulations of aqueous ionic solutions at finite concentrations is that unrealistic cation-cation associations are frequently found. One might suspect a defect in the ion-ion interaction potentials, about which not much is known. However, we show that this phenomenon can also be traced to the fact that, in the pair-potential approximation, the cation-water potentials are too deep compared with the other ones and we investigate this phenomenon in some detail. We then attempt to draw some general conclusions.


Introduction
It is a truism to stress the importance of aqueous salt solutions for many fields of science. Enormous efforts have consequently been made in the past decades to go beyond phenomenological considerations and to gain a fuller understanding of their structure and dynamics at the molecular level. A great many strategies in theoretical and computational physics and chemistry have been deployed for this purpose.
It turns out that while the hydration of almost all simple positive ions has by now been investigated in Molecular Dynamics (MD) simulations at many levels of approximation (many pair potentials (see e.g. [1]), three body potentials [2], various QM-MM methods (see e.g. [3] and references therein), Car-Parinello MD (see e.g. [4]), Brownian dynamics and related approaches [5]), studies of solutions at finite concentrations are much scarcer. This is of course due to the much more difficult modeling task that one faces in this case. At least six interaction potentials (ten pair potentials in this approximation) must be determined in a consistent way. Their relative contributions to the total energy will vary widely when the salt concentration is varied. Any small imbalance between various terms can then lead to artefacts. We will come back to this point below.
We have recently studied aqueous LiCl solutions in their entire concentration range at 300 K and normal densities [6] using a well established class of models that had been used for many salt solutions at low concentrations. Among the findings we observed at intermediate concentrations long-lived (with respect to the simulation times of a few 100 picoseconds) aggregations of Li + -ions. Figure 1 shows such a situation as it is found in a 10 m solution. Partly hydrated Li + -ions aggregate, leaving thus space for small, possibly interconnected, 'pools' of water with the weakly solvated Cl − -ions. The shape of these cationic aggregates is generally more or less linear and the ions are either in direct contact or barely solvent separated, as seen in the figure.
The structural details of very concentrated ionic solutions at room temperature are little known. There is some evidence for the existence of, albeit possibly small (down to 5 H 2 O molecules), nano-pools water, α = 1, see text. The Li + -ions are dark green, the Cl − -ions blue, water molecules closer than 2.5 Å to a Li + -ion (hydration shell) are bright red, other water molecules pink. Hydrogen atoms are omitted. The Li + -clusters together with the hydration water and Cl − -ions close to Li + -ions (ion pairs) have been joined together by a green surface, the 'bulk' water by a brownish one.
of water or worm-like structures in LiCl solutions at low temperatures [7]. Aggregations of some cations have also been found to influence the formation of certain complexes in aqueous solutions at room temperature [8]. Even though they can thus not be entirely excluded, we do not know of any convincing experimental evidence for the existence of such aggregates, concerning either the ions or the solvent, in bulk solutions at higher temperatures. We therefore, decided to consider the Li + aggregates mentioned above as a possible artefact of our model and to study them in some detail in an attempt to trace their existence to some feature of our interaction model. We report here the results of this analysis and new simulations of the same solutions where this phenomenon has been eliminated by various means.

Interaction models, general considerations
The interaction models used in simulations of salt solutions are usually built up by combining a model which has been developed and tested for pure water (e.g. ST2, MCY, SPC, SPCE, TIP4P, BJH) with ion-water potentials developed in various ways, mostly by fitting under different constraints (e.g., the charge of the ion and geometry and partial charges of the water model) empirical functions to ab-initio energies. Simple charged Lennard-Jones (LJ) spheres are mostly, also in this work, assumed for the ion-ion interactions.
Besides the problem discussed in the introductory section, several other ones have been identified. There are many critical discussions of the models for pure water [9][10][11][12][13], but we shall not enter this vast debate here. Pusztai et al. [14] have in particular concluded from comparisons of scattering experiments, reverse Monte Carlo (rmc), and MD simulations that, at higher salt concentrations, the water-water interactions developed for pure water may no longer be adequate.
There is indeed no guarantee whatsoever that interaction potentials of so different origins, often involving different approximations (effective potentials), can be combined in (almost) arbitrary proportions, depending on the salt concentration. Studying single ions, or solutions at low concentration, alleviates this problem. However, most 'real' solutions (biological fluids, sea water, in industry) are such that concentration effects cannot be neglected.
In particular, an imbalance between the cation-water and water-water interaction potentials could lead to artefacts above a certain concentration, i.e. when the ion-water energies are no longer small 23001-2 compared to e.g. the water-water ones. It can then happen that it becomes energetically more favorable to accept some ion-ion repulsion (especially between monovalent ions) and somehow allow these ions to associate so that they can share the remaining (scarce) water molecules to form a common 'hydration shell' for some sort of M n+ n aggregate. We explore this phenomenon (see e.g. in [6]) here by modifying the Li + -water interaction employed in [6] in two different ways: First by reducing the Coulomb-terms in the Li + -O and Li + -H pair potential, and secondly by rigidifying in different ways the water molecule.

Simulation details
The details of the simulation are as in previous work [6] except that either the Li + -water interactions are explicitly modified, or the geometry of the flexible water molecule is fixed, thus removing any mechanical polarization that may exist, see e.g. figure 12 in [6] for a distribution of the water molecular dipoles in LiCl solutions. Consequently, two options were pursued: a) to more or less arbitrarily reduce the Coulomb part of the Li + -oxygen and Li + -hydrogen pair potentials and b) to rigidify the water molecules in a reasonable geometry, see for comparison figure 7 in reference [15] and, as mentioned, figure 12 in [6] for the distributions of angles and the molecular dipole moments of flexible water in solutions. ∠ HOH = 104.5 • , red curves; r OH = 0.9721 Å, ∠ HOH = 102.3 • , green curves; r OH = 0.9872 Å, ∠ HOH = 101.7 • , blue curves; r OH = 1.0115 Å, ∠ HOH = 100.0 • , pink curves, the curves have been shifted by 2 Å in rdirection for better visibility.), for C 2v arrangements. The arrow shows that the potential at α = 0.9 for a deformed water molecule corresponds to the full potential for water in its approximate gas phase geometry.
In all other respects, the simulations were identical to the previous ones: 556 BJH [16] water molecules and, depending on concentration, 10 to 150 ion pairs at experimental densities and 300 K. We have simulated 1 m, 2 m, 4 m, 5 m, 8 m, 10 m, 12 m, and 15 m solutions. The simulations were run essentially under (NV E )-conditions, except that a Kast-type thermostat [17] was very loosely coupled (in order to minimize the perturbations of the dynamics, see below) to compensate an unavoidable numerical error during the ≈10 6 integration steps. Ewald summation was used throughout. Data for the analysis were collected, after extensive equilibration runs, for 125 ps or 250 ps. The average temperature of the runs was 〈T 〉 = (298±1.5) K. In one test case we increased the box size by a factor of 27 (15012 water molecules and 2700 ion pairs, 10 m solution). Starting with the end configuration of a smaller run we could extend this simulation only for about 6.75 ps. No changes in the structural data reported below could be detected in this run.
The reduction in the cation-water interaction was simply achieved by multiplying the Coulomb terms of the Li + -O and Li + -H pair potentials by factors α < 1, which is tantamount to multiplying for these interactions the partial charges of Li + , O and H by α. While this is inconsistent, e.g., from a dielectric standpoint (e.g.; which charges should be used to compute the dielectric constant ǫ?), this is sufficient for our purposes here. We have applied factors of 0.8 α 1.0. Figure 2 shows the energies for a Li + -water supermolecule with C 2v geometry obtained for various values of α and for rigidified flexible BJH water molecules. We have also conducted simulations with such rigidified water (always with full charges, α = 1) using two slightly different geometries: r OH = 0.9572 Å, ∠ HOH = 104.52 • and ∠ HOH = 109.43 • . Figure 3 shows for the flexible models some of the contributions to the total potential energies, divided by the salt concentration for better comparison, as a function of concentration, for the flexible water simulations, with α as a parameter. The middle panel shows how lowering α, i.e. diminishing the Coulomb interaction between the Li + -ion and the water, changes this energy component in the simulations as a function of concentration. The top and bottom panels show how the average water-water and Cl − -water energies react to this modification of the Li + -water interaction. It is seen how lowering (increasing the magnitude) of the cation-water terms 'pushes up' the water-water energies, and that they can even become positive at higher concentrations for the larger αs, as already noted previously [6]. This means that, loosely speaking, the water molecules are, on average, unfavorably oriented with respect to each other (i.e. not adopting hydrogen-bond-type configurations), their orientations being presumably controlled by the cations. The Cl − -water energies (with unmodified pair potentials) increase in magnitude (i.e. are lowered) when the Li + -water terms are diminished. Other energy terms are modified in a similar way, which will not be discussed in detail here. where rc and fc stand for 'reduced charges' (α=0.8, 0.85, 0.9) and 'full charges' (α = 1) in the Li + -water interactions, respectively. By analogy, rc will also be used for the rdf's obtained from rigid water simu-

Results and discussion
lations. Figure 4 also shows ∆n(r ) (right ordinate). The overall effect of α and concentration on ∆n(r ) is studied in a more systematic way in figures 6.     Figure 7 shows that rigidifying the water molecules, i.e., reducing their average intermolecular energies (due to the absence of the mechanical polarization of the flexible model) leads to similar, but smaller changes in the Li + -Li + rdf as explicitly reducing the Li + -water interaction. Similar observations are made for the other g -functions. Figure 8 shows the spectral densities of motion d(ν) for the water molecule oxygens and the Li + -ions in the solutions of lowest and highest concentrations for the four values of α. The spectral density d(ν) is the Fourier-cosine transform of the normalized velocity autocorrelation functionc v v (t ):

Further considerations
The fact that the polarizabilities of the ion and of the water are not, or only in an average fashion and not explicitly, taken into account in the interaction models has often been criticized, especially if such a model is to be used in inhomogeneous environments, e.g. at surfaces [18,19]. However, another effect may also be important and of comparable, if not larger, magnitude: A minimal amount of charge transfer δ|e| (of the order of less than ≈0.03 |e| per water molecule in the first hydration shell of the central ion, say) will result in a discharged ion (in our example Li (1−6δ)+ , on the average) interacting with water molecules charged by the additional δ, which should be apportioned in some way to the partial charges of oxygen (−0.6597|e| in the BJH model) and hydrogen (0.32985|e|). Of course, it is a priori not clear to which site of the model water the charge is best transferred. As figure 9 shows, such a transfer leads in a few plausible cases (e.g. δ|e| transferred entirely to the oxygen site of the water molecule, or δ|e| distributed equally among the three water sites) to markedly shallower ion-water potential wells, comparable to what is seen in figure 2 for different values of α. Note that this plot does not take into account the fact that if several hydration water molecules are present, the ion will be even more discharged, as described above.
In contrast to the present approach, the water-water interactions between first shell water molecules and the other water molecules would also be altered in case of a transfer of charge. The consequences of this will be explored in future work. In any case, these simple considerations show how effective a transfer of charge between an ion and its first shell of solvent molecules might be in modifying the interactions.

Summary
We have studied aqueous LiCl solutions in the entire solubility range of this salt at room temperature and the densities corresponding to ambient pressure. Such systems are particularly interesting in terms of evaluating the underlying models because varying the concentration in a wide range means varying the contributions originating (in the pair potential approach) from various kinds of particle pairs (waterwater, ion-water, ion-ion) in an equally wide range. It was the purpose of this work to show how small modifications of these interactions affect the results of the simulation. We chose first to vary the cationwater interaction since this is where we suspected in earlier work a bias leading to an artefact. Besides, we fixed the geometry of the flexible water model, thus eliminating polarization effects, which occur mainly in the ionic hydration shell, but also elsewhere.
We have shown that decreasing the Li + -water infraction, in an admittedly arbitrary fashion, does indeed lead to qualitatively different results concerning the structure of the solutions, in particular at 23001-7 Ph.A. Bopp, K. Ibuki intermediate concentrations. Rigidifying the water, which leads to a reduction both in the ion-water and (less) the water-water energies, leads to similar, but smaller effects.
This shows how sensitive the results of such simulations can be to details of the interaction model and to the balance between the components of such a 'global' model that may have been taken from very different sources and incorporating different approximations: effective potentials and ab-initio potentials, rigid vs. flexible models, polarization, charge transfer, etc.. The theoretical chemists and solution chemists have their work cut out.