Scientific paper
Effects of Translational and Rotational Degrees of Freedom on the Hydration of Ionic Solutes as Seen by the Popular Water Models
Tomaž Mohoric,1* Urban Bren23 and Vojko Vlachy1*
1	Chair of Physical Chemistry, Faculty of Chemistry and Chemical Technology, University of Ljubljana, Ve~na pot 113,
SI-1000 Ljubljana, Slovenia
2	Laboratory for Physical Chemistry and Chemical Thermodynamics, Faculty of Chemistry and Chemical Technology,
University of Maribor, Smetanova 17, SI-2000 Maribor, Slovenia
3 Laboratory for Molecular Modeling, National Institute of Chemistry, Hajdrihova 19, SI-1001 Ljubljana, Slovenia * Corresponding author: E-mail: vojko.vlachy@fkkt.uni-lj.si Received: 12-12-2014 Dedicated to Prof. Dr. Jože Koller on the occasion of his 70th birthday.
Abstract
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 sol-vation 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.
Keywords: Ionic hydration, water models, degrees of freedom, molecular dynamics
1. 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-4 Hydration of molecules having hydrophobic and ionic groups, such as polyelec-trolytes 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-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, TR, was different than the translational one TR.
In our previous study we used this approach14 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-Jo-nes particle representing a hydrophobe, (ii) the same Len-nard-Jones particle decorated by either positive or negative charge to represent a cation or an anion. An important
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 references14,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-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 Len-nard-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.
2. 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/cm3 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-of-mass positions) and rotational (orientations) degrees of freedom, which enabled us to couple them to separate thermostats. As before14,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 aSS equal to 3 À and eSS equal to eWW. Lo-rentz-Berthelot mixing rules were used to compute oWS and eWS 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.
3. Results and Discussion
The simulation results are presented in the form of solute-water RDF (gSO) 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.
3. 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-
Figure 1. Singly charged cation (+1)-water RDFs for water models of group I (left) and group II (right). For the sake of clarity, the results for different water models are shifted vertically and horizontally. Different lines correspond to different values of TR: 300 K (solid red), 500 K (dotted green) and 700 K (dashed blue). Arrows indicate the direction of the RDFs variation with a rotational temperature increase. Inset: angular distribution functions of OH bonds in the first hydration shell of the solute for group I (left) and group II (right) water models, i.e. SPC/E and TIP5P respectively. For TIP5P also the angular distribution of OL bonds is shown.
ce: an increase in TR 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 TR 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.
Figure 3. The same as Figure 1 but with the prescribed rotational temperature applied only to the water molecules in the cation's first hydration shell.
Figure 4. The same as Figure 1 but with the prescribed rotational temperature applied only to all water molecules outside of the cation's first hydration shell (the bulk water).
b) The indirect effect. In aqueous systems the orientation dependent interaction - the water-water interaction - is always present. An increase of TR 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 ro-tationally 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 gSO 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 8~0, indica-
Figure 5. The same as Figure 1 but for an anion(-l).
Figure 6. The same as Figure 2 but for an anion(-2).
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 TR. 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 gSO 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 gSO for both TIP5P and ST2 models. Analogously, for models of class I the first peaks of solute-water RDFs exhibit weaker dependence on TR, as the two opposing effects largely cancel each other.
3. 2. Solvation of Anions
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 solva-tion 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.
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 solute-water 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.
3. 3. 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 TR (Figure 10, top). Furthermore, the first peaks
- I	' TR [k]: 300 ' 500 \j 700 V V*
	1 i 1 1
■j	I 1 I I \ ojcr
	I I I I I
	! 1 1 ! 1
	
1	i i i i i
2	4	6	8
rlA]
Figure 9. The fixed water-water RDFs at different rotational temperatures for the SPC water model. Top: the total RDF. Middle: RDF around the cationic part of the fixed water. Bottom: RDF around the anionic part of the fixed water. Arrows indicate the direction of the RDFs variations with the rotational temperature increase; the colour code as for Figure 1. Inset shows how the solute molecule was divided into the cationic and the anionic part.
around both the cationic and anionic parts of the solute decrease with the increasing TR (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.
2 -
'c та Ö L/l
Cn
1 -0 -
2	4	6	8
r [A]
Figure 10. The same as Figure 9 but for the TIP5P water model.
4. Conclusions
We investigated how different water models respond to independent variations in strengths of the translational and rotational degrees of freedom. We studied hydration of simple ionic (cations and anions) and polar solutes using six popular water models, which can be, with respect to the charge distribution, classified into two groups. We showed that an increase in the rotational temperature yielded different effects on the cation-water correlation, depending on the chosen water model. As far as the anion-water correlations were concerned, all the models exhibited qualitatively similar behaviour. We attempted to explain these results by the competing (direct and indirect) effects caused by an increase in TR. The solute-water angular distribution function, strongly influenced by the
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 example24'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 rotatio-nally 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.
5. Acknowledgement
This study was supported by the Slovenian Research Agency (ARRS) funds through the Program 0103-0201, Projects J1-4148 and J1-5448, and the Young Researchers Program (T.M.) as well as by the NIH research grant (GM063592). T.M. acknowledges Daan Frenkel for providing access to computer clusters at the Department of Chemistry, University of Cambridge.
6. References
1.	J. M. G. Barthel, H. Krienke, W. Kunz, in: H. Baumgartel, E. U. Franck, W. Grunbein (Eds.): Physical Chemistry of Electrolyte Solutions: Modern Aspects, Topics in Physical Chemistry Vol. 5, Springer, New York, 1998, pp. 1-53.
2.	Y. Marcus, G. Hefter, Chem. Rev. 2006, 106, 4585-4621. http://dx.doi.org/10.1021/cr040087x
3.	Y. Marcus, Chem. Rev. 2009, 109, 1346-1370. http://dx.doi.org/10.1021/cr8003828
4.	A. Ben-Naim, in: Molecular Theory of Water and Aqueous Solutions. Part I: Understanding Water, World Scientific Publishing Co. Pte. Ltd., Singapore, 2010, pp. 280-564.
5.	N. J. English, J. M. D. MacElroy, J. Chem. Phys. 2003, 118, 1589-1592. http://dx.doi.org/10.1063/L1538595
6.	N. J. English, J. M. D. MacElroy, J. Chem. Phys. 2003, 119, 11806-11813. http://dx.doi.org/10.1063/L1624363
7.	N. J. English, D. A. Mooney, J. Phys. Chem. B. 2009, 113, 10128-10134. http://dx.doi.org/10.1021/jp902500m
8.	N. J. English, D. A. Mooney, S. O'Brien, Mol. Phys., 2011, 109, 625-638.
http://dx.doi.org/10.1080/00268976.2010.544263 9. N. J. English, G. Y. Solomentsev, S. O'Brien, J. Chem. Phys. 2009, 131, 035106-1-10.
10.	N. J. English, G. Y. Solomentsev, D. A. Mooney, J. Chem. Phys. 2010, 133, 235102-1-7.
11.	J.-A. Garate, N. J. English, J. M. D. MacElroy, Mol. Simul. 2009, 35, 3-12.
http://dx.doi.org/10.1080/08927020802353491
12.	J.-A. Garate, N. J. English, J. M. D. MacElroy, J. Chem. Phys. 2011, 136, 055110-1-12.
13.	M. Tanaka, M. Sato, J. Chem. Phys. 2007, 126, 034509-1-9. http://dx.doi.org/10.1063/L2403870
14.	U. Bren, D. Janežič, J. Chem. Phys. 2012, 137, 024108-1-11. http://dx.doi.org/10.1063/L4732514
15.	T. Mohorič, B. Hribar-Lee, V. Vlachy, J. Chem. Phys. 2014, 140, 184510-1-7.
16.	M. O. Jensen, V. Jogini, D. W. Borhani, A. E. Leffler, R. O. Dror, D. E. Shaw, Science 2012, 336, 229-233. http://dx.doi.org/10.1126/science.1216533
17.	R. O. Dror, et al., Nature 2013, 503, 295-299.
18.	Y. Shan, et al., Nature Structural & Molecular Biology 2014, 21, 579-584. http://dx.doi.org/10.1038/nsmb.2849
19.	H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, J. Hermans, in: B. Pullman (Ed.): Intermolecular forces, Rei-del, Dordrecht, 1981, pp. 331-342. http://dx.doi.org/10.1007/978-94-015-7658-1_21
20.	W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Im-pey, M. L. Klein, J. Chem. Phys. 1983, 79, 926-935. http://dx.doi.org/10.1063/L445869
21.	H. J. C. Berendsen, J. R. Grigera, T. P. Straatsma, J. Phys. Chem. 1987, 91, 6269-6271. http://dx.doi.org/10.1021/j100308a038
22.	F. H. Stillinger, A. Rahman, J. Chem. Phys. 1974, 60, 1545-1557. http://dx.doi.org/10.1063/L1681229
23.	M. W. Mahoney, W. L. Jorgensen, J. Chem. Phys. 2000, 112, 8910-8922. http://dx.doi.org/10.1063/L481505
24.	S. Rajamani, T. Ghosh, S. Garde, J. Chem. Phys. 2004, 120, 4457-4466. http://dx.doi.org/10.1063/L1644536
25.	D. L. Mobley, A. E. Barber II, C. J. Fennel, K. A. Dill, J. Phys. Chem. B. 2008, 112, 2405-2414. http://dx.doi.org/10.1021/jp709958f
Povzetek
S pomočjo simulacije gibanja molekul v sistemu, kjer smo posameznim prostostnim stopnjam (translacija, rotacija) predpisali različne temperature, smo raziskali vpliv teh prostostnih stopenj na hidratacijo enostavnih ionov. Zanimalo nas je, kako se, glede na izbrani model vode, hidratacija enostavnih ionov odziva na povišanje rotacijske temperature. Pokazali smo, da se modeli vode, skladno s porazdelitvijo negativnega naboja na molekuli, ločijo v dve skupini. Do različnega odziva zato prihaja predvsem pri hidrataciji kationa. Za 3- in 4-točkovne modele vode (npr. SPC model), kjer je negativni naboj v središču molekule vode, hidratacija kationov z naraščujočo temperaturo rotacije narašča. V nasprotju s tem pa pri 5-točkovnih modelih vode (npr. TIP5P), hidratacija kationov z višanjem temperature rotacije slabi. Pri vseh obravnavanih modelih vode se hidratacija anionov zmanjšuje s povišanjem rotacijske temperature. Razlog je podobna porazdelitev pozitivnega naboja pri vseh modelih vode. Vpliv translacijske temperature je pri vseh topljencih in za vse modele vode kvalitativno enak - vsako povišanje temperature slabi korelacijo med topljencem in vodo