Effects of Translational and Rotational Degrees of Freedom on the Hydration of Ionic Solutes as Seen by the Popular Water Models

We employed molecular dynamics simulations with separate thermostats for translational and rotational temperatures in order to study the effects of these degrees of freedom on the hydration of ions. In this work we examine how water models, differing in charge distribution, respond to the rise of rotational temperature. The study shows that, with respect to the distribution of negative charge, popular water models lead to different responses upon an increase of the rotational temperature. The differences arise in hydration of cations, as the negative charge distribution on the model solvent represents the determining factor in such cases. The cation-water correlation increases with the increasing rotational temperature if negative charge is placed in (or close to) the centre of the water molecule (a typical example is the SPC water model) and decreases, when the negative charge is shifted from the centre (as in the TIP5P model of water). Because all the water models examined here have similar distributions of positive charge, they all exhibit similar trends in solvation of anions. In contrast to above, the effect of translational temperature variation is similar for all water-solute pairs; any increase in translational temperature decreases the solute-water correlations.


Introduction
Properties of solutions are not determined solely by solute-solute interactions but also by the ability of solvent to solvate the solute particles. [1][2][3][4] Hydration of molecules having hydrophobic and ionic groups, such as polyelectrolytes and proteins, is relevant for electrochemistry, geochemistry and especially for biology. Interesting questions regarding the role of different fundamental degrees of freedom (translation, rotation and vibration) in hydration have been addressed recently. A relevant physical situation may occur upon interaction of microwaves with aqueous solutions and has been studied by the non-equilibrium molecular dynamics simulation with the electric field modeled explicitly. [5][6][7][8][9][10][11][12][13] Bren and Jane`i~ considered another approach to this problem, 14 employing the steady-state molecular dynamics simulations with separate thermostats for different degrees of freedom. They examined the hydration of hydrophobes and cations under conditions where the rotational temperature, T R , was different than the translational one T R .
In our previous study we used this approach 14 to simulate the solvation of cations and anions (the latter were not studied before) for a range of conditions. 15 We used the molecular dynamics simulations with separate thermostats to examine effects of translational and rotational degrees of freedom on the structure of water in solution.
To describe water molecules we used the standard SPC/E model, while the model solutes were (i) the Lennard-Jones particle representing a hydrophobe, (ii) the same Lennard-Jones particle decorated by either positive or negative charge to represent a cation or an anion. An important Mohori~ et al.: Effects of Translational and Rotational Degrees ... conclusion of our study was that, while an increase of the translational temperature always decreases solute-water correlations, higher rotational temperature has varying effects, depending on nature of the solute. In summary, increased rotational temperature affected hydration of cations in an opposite way to anions. Our study also indicated that charge distribution on a solvent molecule could be important. This prompted us to explore to what degree the results obtained in references 14,15 depend on the chosen model of water. The question we wish to answer here is the following: can we expect the same trends as obtained for SPC/E also for other water models, or the results are model-dependent? For this reason we examined the behaviour of six popular water models, most of them are frequently used in the large-scale molecular-dynamics (MD) simulations of aqueous solutions, [16][17][18] in situations where the rotational temperature of the system is different than the translational one.
The models mimicking water molecules can be divided into two distinct groups. In most water models the negative charge coincides with (or is put very close to) the Lennard-Jones centre (e.g. SPC, SPC/E, TIP3P, TIP4P) 19-21 -let us call it group I. However, some water models try to mimic the oxygen's electron lone pairs by putting two negative charges further away from the Lennard-Jones centre (e.g. TIP5P, ST2) 22,23 -this is our group II. Distribution of the positive charge (mimicking the two hydrogen atoms) is roughly the same for all water models: the two positive charges are placed around 1 Å from the Lennard-Jones centre with the angle between them ranging from 104° to 109° degrees.
The differences between the water models of types I and II are expected to be reflected in hydration of polar and ionic species, as the solvent charge distribution is crucial in such cases. 24,25 In this work we simulate the hydration of singly and doubly charged ions and of a polar molecule (represented by a model water molecule), by applying six popular water models mentioned above. We have calculated the solute-water radial distribution functions (RDFs) and angular distribution functions (ADFs) for all these examples. The paper is organized as follows: after this Introduction we provide the model and simulation details, followed by the presentation and discussion of results. A brief summary of the most interesting results is placed at the end.

Model and Simulation Details
Molecular dynamics simulations were performed in a canonical (N, V, T) ensemble. The system consists of a fixed solute, placed in the centre of the simulation cell and many surrounding water molecules. When the solute was a singly charged ion or a polar molecule, 256 water molecules were placed in the system. For the doubly charged ions, 512 water molecules were used in the calculations. A simple spherical cut-off with a distance 9 Å was used for Lennard-Jones, as well as, for electrostatic interactions. 15 We have verified, that results did not change upon an increase of the cut-off distance.
The density of water molecules corresponded to 1.0 g/cm 3 and the time step used was 1 fs. The solutions were simulated for 5 ns, of which the first 0.1 ns were an equilibration period, followed by the 4.9 ns long production run, over which the statistics was collected. Equations of motion were solved separately for translational (centre-ofmass positions) and rotational (orientations) degrees of freedom, which enabled us to couple them to separate thermostats. As before 14,15 we used the velocity rescaling procedure at every time step to keep the two temperatures fixed. We carefully verified that results were not affected by the way how we thermostatted the system: temperature can be rescaled in every time step, or it can be held constant on the average. Because the velocity rescaling in every time step provides better control over the temperature values, we adopted the latter approach. 15 We applied the same protocol for the equilibration and production periods.
Ions were represented by the Lennard-Jones (LJ) particles with σ SS equal to 3 Å and ε SS equal to ε WW . Lorentz-Berthelot mixing rules were used to compute σ WS and ε WS values. The charges were located in the centres of ions. Polar solutes, studied here, were represented as the respective (model) water molecules. The water models considered here included the SPC, SPC/E, TIP3P, and TIP4P models as members of the group I, and TIP5P and ST2 models as the members of the group II. In simulations we kept one set of degrees of freedom at constant temperature of 300 K, while the other one assumed values of 300, 500, 600, and 700 K. For the sake of clarity the RDFs and ADFs for 600 K, which always fall between the results for 500 K and 700 K, are omitted from the graphs shown here.

Results and Discussion
The simulation results are presented in the form of solute-water RDF (g SO ) and ADF (dP/dθ) for the OH bonds (OL bonds in the case of the group II models, where L refers to the oxygen's electron lone pair) lying in the first hydration shell of a solute. In accordance with previous observations, 15 we find that faster translational motion is always associated with a decrease in solute-water correlation. Because the figures, showing the effects of the translational temperature do not provide new insights, we decided not to include them here.

1. Solvation of Cations
Effects of the rotational temperature variations on the structure of water solvating +1 and +2 charged cations are shown in Figures 1 and 2. General trends are easy to noti-  ce: an increase in T R causes for cations to become better solvated (the first peak of the solute-water RDF increases) in water models of group I and less solvated (the first peak of the solute-water RDF decreases), if water models of group II are applied. We can identify two opposing effects of an increase of T R to the solute-water correlation: 15 a) The direct effect. The decrease of the average solute-water correlation takes place when the solute-water interaction is orientation dependent. In case of a strong orientation dependence, faster rotational motion causes weaker hydration of a solute.  b) The indirect effect. In aqueous systems the orientation dependent interaction -the water-water interaction -is always present. An increase of T R makes the average water-water interaction weaker, improving water's ability to hydrate the solute. We attempted to partly separate one effect from the other by heating up only the water molecules in the solute's first hydration layer (mostly the direct effect) or all the others (the indirect effect) waters. Notice that every such separation is necessarily arbitrary and is here proposed merely as an attempt to better illustrate the effects of an increase in the rotational temperature. As revealed by Figure 3, the water models of group I and II are clearly distinguishable by their response to the rise of the rotational temperature in the cation's first hydration shell. While rotationally excited group I water molecules tend to accumulate in the cation's first hydration layer (the first peak of the cation-water RDF increases), the group II water molecules are escaping from it (the first peak of the cation-water RDF decreases). On the other hand, Figure 4 shows that an increase of rotational temperature of the bulk water molecules has the same effect for both groups of water models -causing a slight increase of the first peak in the cation-water RDF. In our approximate picture the stronger of the two effects determines the response to rise of the rotational temperature.
The insets in Figure 1 show ADFs for typical water models of group I (e.g. SPC/E), and group II (e.g. TIP5P). Angular distribution of OH bonds is relatively broad in the first case, and an increase of rotational temperature has only a marginal effect on it. The water-water interaction exhibits strong orientation dependence, so the indirect effect wins and the first peak in g SO increases, even when heating up only the first hydration shell water molecules (see Figure 3). Different angular distributions are observed for the group II water models, i.e. the one belonging to the OL bonds. In this case the distribution has a strong and narrow peak at θ∼0, indica-     ting that water favours orientations with lone pairs pointing directly to the cation. The OL bond distribution simultaneously decreases and widens substantially upon an increase in T R . This effect is strengthening the direct influence. Around singly charged cations (+1) both effects seem to largely cancel each other in TIP5P water, while in ST2 water the direct influence prevails and the first peak of g SO decreases. In ST2 water model the negative charge is placed at a distance of 0.8 Å, while in TIP5P at a distance of 0.7 Å from the oxygen centre.
This rationalizes the observation that the direct influence is stronger in the ST2 water. In the case of doubly charged (+2) cations, the angular distributions of group II water models become even more pronounced and narrower, thus strengthening the direct effect further. This results in a decrease of the first peak of g SO for both TIP5P and ST2 models. Analogously, for models of class I the first peaks of solute-water RDFs exhibit weaker dependence on T R , as the two opposing effects largely cancel each other. Figures 5 and 6 indicate that hydration of anions yields qualitatively the same behaviour for all the models examined here. In the same manner as for the cations above, we can try to separate the direct and indirect effects by applying the elevated rotational temperatures to the water molecules in the anion's first hydration layer or, alternatively, to all the others. Figures 7 and 8 reveal the same qualitative trends for both groups of water models: a) a substantial decrease in the first peak of the anion-water RDF when increasing the rotational temperature of the anion's first hydration shell water molecules (due to the direct influence) and b) a slight increase in the first peak when applying higher rotational temperature to all the other water molecules (the indirect influence). This is to be expected as all the water models have very similar distributions of the positive charge (mimicking the two hydrogens), which represents the most important factor in the solvation of anions. The angular distributions of OH bonds (see insets of Figures 3 and 4) posses a strong and narrow peak at 0°. Higher rotational temperatures decrease this peak significantly, causing the direct influence to be the dominant factor.

2. Solvation of Anions
The quantity often used for characterization of hydration of solutes is the hydration number. In our case this is the number of water molecules that lie within the solute's first hydration shell. The outer boundary of the shell is given by the position of the first minimum in the solutewater RDF. The changes in hydration numbers observed in our simulations are rather small, not much bigger than the estimated numerical errors. Upon an increase of the rotational temperature from 300 to 700 K, the hydration numbers for cations (+1) as well as for anions (-1), increase from less than 1% (TIP3P and ST2 water models) to about 4% (TIP5P and TIP4P water models). Interestingly, the qualitatively different response of the first peak in the solute-water RDF for cations and anions (or for different water models) upon increase of the rotational temperature is not reproduced in changes of hydration numbers. The latter, namely, slightly but systematically increase in all the examples considered here. An exception is an anion solvation by the ST2 water model, where the hydration number decreases by less than 1% or it remains unchanged.

Solvation of Polar Molecules
In aqueous solutions, water molecules have an option to solvate either the solute or other waters. In this subsection we consider the hydration of a model polar compound, represented by a water molecule typical of either group I or II. We define the cationic and the anionic part of the solute in the following way: the water molecule is fixed in the centre of the simulation box with its symmetry axis pointing in the z-direction. To obtain the solute-water RDF for the cationic part we consider only the distances with the molecules in the half of the box, which includes solute's hydrogen atoms. In opposite to this, we use the molecules from the other half of the box to evaluate the anionic RDF (see Inset in Figure 9). Figure 9 shows the solute-water RDFs as a function of rotational temperature for the SPC water model. Trends are similar as in our study of the SPC/E water model. 15 Total RDF (Figure 9, top) does not change significantly, which can be attributed to the opposing direct and indirect effects discussed before, that to a large extent cancel each other out. With an increase in rotational temperature, the first peak of RDF increases next to the cationic part (Figure 9, middle) of the solute molecule, and decreases in the vicinity of the anionic part (Figure 9, bottom).
As expected, the trends are different when we consider the hydration of a fixed TIP5P molecule in an environment of other TIP5P water molecules. Here the first peak of the total RDF monotonically decreases with the increasing T R (Figure 10, top). Furthermore, the first peaks around both the cationic and anionic parts of the solute decrease with the increasing T R (Figure 10, middle and bottom). The SPC and the TIP5P water molecules solvate the positive and negative part of the polar solute in the same qualitative way, as above seen for cations and anions. charge distribution on the solvent molecule, determines, which of the two effects dominates. At this point we also wish to mention that a variation of the translational temperature may screen the subtle effects of the rotational temperature changes, as it exhibits a stronger effect on the solute hydration, regardless of the water model. 14,15 Different behaviour of the two different groups of water models is also reflected in the hydration of polar solutes, which is in accordance with the solvation behaviour of these models for simple ions.
Rotationally "excited" water molecules may hydrate the solute better for some water models and worse for the others. In other words: the MD results, obtained at elevated rotational temperature, may be model dependent. This confirms previous observations (see, for example 24,25 ) that different water models can perceive the same thermodynamic and/or structural property differently. The experiments that would, as much as possible, mimic the rotationally excited water conditions of this study, would therefore provide a better understanding of the water physics and the hydration phenomenon. To the best of our knowledge no such experiments have been published so far.