Chemical Physics Letters 538 (2012) 60–66
Contents lists available at SciVerse ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Car–Parrinello molecular dynamics simulations of Na+ solvation in water, methanol and ethanol Yongping Zeng a,⇑, Junmei Hu a, Yu Yuan a, Xiaobin Zhang b, Shengui Ju c a
College of Chemistry and Chemical Engineering, Yangzhou University, Yangzhou 225002, China Information Engineering College of Yangzhou University, Yangzhou 225002, China c College of Chemistry and Chemical Engineering, Nanjing University of Technology, Nanjing 210009, China b
a r t i c l e
i n f o
Article history: Received 10 January 2012 In final form 19 April 2012 Available online 27 April 2012
a b s t r a c t Car–Parrinello molecular dynamics simulations have been performed on Na+ in water, methanol and ethanol. The structure of the first solvation shell around the Na+ is in a good agreement with the experimental data. It has been found the Na+ has a coordination number (5.13) in water greater than that in methanol and ethanol. However, the decay of orientational correlations for solvation shell methanol and ethanol molecules shows slower relaxation compared to that of bulk. Translational dynamics of Na+ in these systems is affected in a small degree by molecular size. It produces a significant difference on their relative relevance bands. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction Solvation is a universal phenomenon playing an important role in a wide range of chemical and biochemical processes. The solvation structure and dynamics of ions in polar solvent influence a wide range of chemical physics, with implications in medical, biochemical and technological research. To study the ion solvation, especially in water, many experimental have been applied including neutron [1,2], X-ray scattering [3], infrared and Raman spectroscopy [4,5] and, more recently, femtosecond spectroscopy [6,7]. The polarization-selective IR pump–probe and two dimensional infrared (2DIR) spectroscopy have been used to study the H-bond dynamics and the spectral diffusion dynamics [8]. While, even for some fundamental properties such as the average distance and the average coordination numbers for the first solvation shell, the results widely depend on the experimental and theoretical method of investigation. Characteristic features of the solvation shell, an integral part of a hydrated ion, have been intensely studied experimentally and theoretically [9]. However, experimental determination of a detailed microscopic structure of solvation shells of ions is very difficult [10]. Classical computer simulation techniques, such as molecular dynamics or Monte Carlo can provide detailed and statistically reliable information on structure and dynamical properties of solvent molecules around the ions. While, these classic simulation methods rely on the quality of the interaction force field parameters. These
⇑ Corresponding author. E-mail address:
[email protected] (Y. Zeng). 0009-2614/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2012.04.035
force field functions may be parameterized by fitting the results from available experimental data [11] or quantum chemical calculations [12]. However, the validation of empirical potentials in the liquid phase condition is difficult. On the other hand, the transferable force field parameters are still difficult for different simulation systems. Car–Parrinello molecular dynamics (CPMD) simulations have been shown to be a powerful tool to obtain microscopic insight into some solution systems [13–16]. In recent years, the Car–Parrinello molecular dynamics based on the density-functional theory (DFT), is one of the more popular alternatives. The CPMD technique has been successfully applied in simulations of both pure water [17,18] and ions in water [19–21] as well as many other systems of chemical interest [22–24]. The essence of ab initio molecular dynamics consists in calculating the forces between the nuclei are determined ‘on the fly’ from quantum chemical calculations at each time step, rather than fixing the potential at the outset of the simulation. The CPMD could, in principle, provide information on the dynamics and structure of molecules on a theoretically solid basis. In addition, this information can be used for derivation of more reliable effective potential functions for classical molecular simulations [25]. In this Letter, dilute solutions of Na+ in three different solvents, water, methanol and ethanol are studied by CPMD simulations. Water represents a strongly polar molecule while methanol and ethanol are intermediate between nonpolar and strongly polar nonaqueous molecules. For methanol, the smallest molecule characterized by both a hydrophobic and a hydrophilic group, is widely used as a solvent, but the number of computational studies on the structural and dynamical properties of the liquid [26,27] and its
Y. Zeng et al. / Chemical Physics Letters 538 (2012) 60–66
ionic solutions [28,29] is still limited. In the past few years, the interest in methanol has also grown as a possible fuel cell component [30,31]. Therefore, the comprehension of the liquid methanol properties and its interactions with ions and simple molecules is becoming important. Comparison with methanol, the ethanol molecule has a stronger hydrophobic group which maybe has a significant effect on the structure of solution. Ethanol has distinct solvation properties, and it is easily soluble as its polar hydroxyl group can participate in the hydrogen bonded network and the hydrophobic ethyl group is relatively bulky. Consequently, the understanding of the ethanol properties and its interactions with Na+ ion is an interesting topic. At present, we report results from Car–Parrinello molecular dynamics simulations of Na+ in water, methanol and ethanol. Particular attention has been paid to the structure of the first solvation shell and to the dynamical properties of the ions and of the solvent molecules. The details of simulation methods are described in Section 2 and our results are discussed in Section 3. This section is subdivided into various subsections on structural and dynamical properties. Finally, our conclusions are summarized in Section 4.
2. Computational methods The simulations have been performed with the CPMD code [32] for the present work. Constant volume Car–Parrinello molecular dynamics simulations have been performed on solutions containing a single Na+ ion and either 32 water molecules or 25 methanol or 25 ethanol molecules in cube boxes applied with periodic boundary conditions in all three directions. The size of cubic boxes is 9.89 Å for water, 11.93 Å for methanol and 13.60 Å for ethanol, respectively. In this work, Troullier and Martins [33] pseudopotentials have been used along with Kleinman and Bylander [34] decomposition for the C and O atoms, whereas a Von Barth–Car pseudopotential for the H atoms [35]. The choice has been shown to give good results for pure liquid methanol [36]. Calculations with a Martins–Troullier pseudopotential have been also performed for comparison on the structure properties and the results are in the expected agreement for this kind of simulations [17]. The metal ions are taken into account using Goedecker pseudopotential [37] which gives both structure and electronic satisfactory results. The electronic states are expanded in a plane-wave basis have been truncated at 70 Ry. The time step in the numerical integration of the equation of motions is 5 a.u. (0.12 fs). A fictitious mass of 800 a.u., is assigned to the electronic degrees of freedom. The Kohn–Sham formulation of DFT is performed to calculate the electronic structure. In the DFT calculations, we use the gradient-corrected BLYP exchange correlation functional. The type of the functional has been used previously to study the liquid water [38] which the structure and dynamics properties were described very well. In order to be able to use a relatively larger time step, the hydrogen atoms of solvent are fully substituted by deuterium. The ionic motions are less affected by this choice [39]. In the liquid methanol and ethanol, the hydrogen bonds are the dominant interactions as in liquid water. The initial configuration for the ab initio simulations is taken from a classical molecular dynamics simulation that was run for 1 ns using empirical OPLS potentials [40]. We perform 3 ps of thermalization of this system at 300 K by velocity scaling, and then, 10 ps of free dynamics to make the systems equilibrated. Finally, 32 ps in the NVE ensemble is run, storing the cell dipole every five steps and the atomic coordinates and velocities every step during this production phase of the simulation run for our structural and dynamical analyses.
61
3. Results and discussion 3.1. Structure properties The radial distribution functions (RDFs) for Na O are shown in Figure 1, along with their running integration numbers.
nðrÞ ¼
4pN V
Z
r
gðrÞr 2 dr 0
where N is the number of oxygen atoms, and V is the volume of the simulation cell. The distribution functions show at least two well resolved peaks in all three systems. The most distinguished feature is the sharp first maximum of Na+ O RDF. In addition, the amplitude of the Na+ O RDF in the first minimum dose not decay to zero in a rather large interval for the three systems. Thus the first solvation shell of the Na+ ion is not very stable. To leave the first solvation shell for solvent molecules is not very difficult. This is why exchanges of solvent molecules between the first solvation shell and the bulk are frequent. Simultaneously, we find that the intensity of the RDF increases from water to methanol and ethanol. This suggests that stability of the first solvation shell is ethanol > methanol > water. The higher first peak in the Na O RDF for ethanol state that there is more structuring of the first shell in ethanol relative to bulk ethanol than in water or methanol. However, this results largely from the lower bulk density of oxygen atom in ethanol since the molar volumes are in the order water < methanol < ethanol. So to incorporate a molecule of ethanol is the most difficult in these three systems. This may be attributed to the steric hindrance of the bulky –C2H5 groups. Examination of gNa–O(r) for water in Figure 1 reveals two prominent peaks centered at 2.45 and 4.43 Å, which correspond to the first and second solvation shells around the sodium ion. The position of this secondary peak, at 4.43 Å, suggests that the second solvation shell is hydrogen bonded to the first shell. There is also evidence of subsidiary structure around 5.5 Å, which suggests that the influence of the Na+ ion extends quite far into the solvent. The first peak of Na+ O for the methanol and ethanol are 2.37 and 2.29 Å, respectively. The first peak positions for the three systems had minor differences. The coordination number can directly be obtained from the integration over the first maximum of the Na+ O RDF of water, methanol and ethanol. The first Na+ O peak of the three systems is respectively at 2.45, 2.37, and 2.29 Å. This indicates that the groups (H, methyl, ethyl) have influence on the first peak position. It can be seen that the cutoff distance has an effect on the coordination
Figure 1. Radial distribution functions and relative integration numbers for Na+ O.
62
Y. Zeng et al. / Chemical Physics Letters 538 (2012) 60–66
Figure 3. Na+ O distances refer to solvent molecule in the first coordination shell for Na+ in ethanol along the simulation run, (a) Na+ in water; (b) Na+ in methanol; and (c) Na+ in ethanol (Only distances corresponding to the solvent molecules that belong to the first solvation shell at any time along the simulation are displayed).
Figure 2. The coordination number along the simulation runs (a): water; (b): methanol; and (c): ethanol.
number due to the fact that the first minimum in the g(r) is widespread for all three solutions. For Na+ in water, X-ray diffraction
studies [41] locate the first peak position at 2.443 Å with a ‘relatively rigid octahedral’ cage thus proposing a 6-fold coordination. Rempe and Pratt [42] observed that the coordination number is 4.6 on average which is slightly smaller than 5.13(this work). White et al. applied the first principles simulation to obtain the first peak position for Na+ solutions is found at a slightly larger (2.49 Å) [43] The results of the present simulation are in agreement with experiments [41] and simulations [43]. In this work, the first peak
Y. Zeng et al. / Chemical Physics Letters 538 (2012) 60–66
63
Figure 4. Spatial distribution of Na+ ion relative to a solvent molecules in the first solvation shell (top panel) and oxygen atoms of the first solvation shell relative to Na+ ion (bottom panel). (a): water (red to intensity 80; blue to intensity 30); (b): methanol (red to intensity 100; blue to intensity 60); (c): ethanol (red to intensity 100; blue to intensity 40); (d): water (red to intensity 50; blue to intensity 20); (e): methanol (red to intensity 100; blue to intensity 50); and (f): ethanol (red to intensity 120; blue to intensity 80). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
position of Na+ O RDF is similar to the Na O distances of 2.38 Å and 2.4 Å found for sodium ion in water. For Na+ in methanol, the first peak position of Na+ O RDF is similar to the Na+ O distances of 2.35 Å found for sodium in methanol by Jorgensen et al. [44] They observed that the coordination number is 6, which is greater than the present work. While Dauster et al. [45] reported a 5-fold coordinated Na+ in the methanol, which is very close to this work. This is very different to water and methanol probably related to the more bulky –C2H5 alkyl group for ethanol. We note in this study that the 5- and 6-fold coordinated Na+ in the water and methanol cases were not conclusive for the ethanol solution. Forck et al. [46] found that 5- and 6fold coordinated Na+ could not compete with the optimization hydrogen bond network interactions resulting in minima with a 4-fold coordinated Na+. These observations are in good agreement with our results of coordination number. It is important to note that the coordination number obtained with CPMD simulations for Na+ ions in methanol and ethanol is lower than the corresponding value in water, essentially due to the greater steric effects of the methanol and ethanol molecules. To obtain an insight of solvent exchange, coordination numbers in the first solvation shell of Na+ ion have been calculated from the trajectories. Figure 2 shows as a function of time the molecules bound to Na+ in the first solvation shell. The molecule is defined to belong to the first solvation shell if the distance of its oxygen atom from the Na+ ion is smaller than the first peak position in the RDF (3.3 Å). The number of first-shell solvent molecules fluctuated frequently for water and methanol during 12 ps time scale. This finding confirms the slow exchange of ethanol molecules in the primary solvation shell. The exchange of waters between solvation shells can be analyzed in detail by monitoring the Na+ O distance. Figure 3 shows a plot of the distance of the oxygen atoms of water molecules from the Na+ ion as a function of time. Only those solvent molecules which move to within a distance of 3.3 Å of the Na+ are included
in the figure. Throughout the period of the simulation, we find that there are some exchanges for the solvent molecules of the first solvation shell with outer solvent molecules during our simulation time scale. The results show that the solvation shell of Na+ in water and methanol have a short-lived structures. A molecule is defined to belong to the first solvation shell if the distance of its oxygen atom from the Na+ ion is smaller than the first peak position in the g(r) function. It can be seen that some of the molecules are permanently in the shell; the other molecules are almost constantly in the shell and other molecules alternate as first neighbors. The overall balance leads to a coordination number of 5.13, 4.4 and 4.15 for water, methanol and ethanol respectively. This can be attributed to the larger ionic radius of Na+ that allows, on average, a higher number of first neighbor solvent molecules. It can be seen that ethanol molecules in the shell and other molecules alternate as first neighbor molecules. The overall balance leads to a coordination number which is slightly more than 4. This can be attributed to the larger size of ethanol molecules that results in, on average, a lower number of first neighbor molecules. Inspection of Figure 3 reveals that four water molecules clearly remain within the first coordination shell, while two molecules make exchanges between the first and second shells. One of the molecules involved in the exchange resides most in the second shell. The simulations exhibit both associative as well as dissociative interchange mechanisms, i.e., exchange mechanisms that proceed through intermediate states where the ion–water coordination number is either larger or smaller than the average coordination. For methanol molecules, four of them remain within the first coordination shell, while three molecules make exchanges between the first and second shells. The structure observed at the beginning and end of the simulation scale shows that some significant fluctuations in the first shell distances suggesting a destabilization of the solvation structure. It is obvious that two molecules make exchanges at about 9.5 ps in our simulation time scale. For ethanol molecules, there are three molecules within the first coordination shell, while
64
Y. Zeng et al. / Chemical Physics Letters 538 (2012) 60–66
Figure 5. The first and second rank orientational correlation functions of the vector of OH bonds of water, methanol and ethanol molecules. Results are presented for water, methanol and ethanol molecules in the first solvation shell of Na+ ion (upper panel), bulk water, methanol and ethanol molecules (lower panel).
three molecules make exchanges between the first and second shells. The figure shows that two ethanol molecules make an exchange at about 3.2 ps for this time scale. However, in our simulation time scale, it still is not enough to find the rate of exchange and its mechanisms. So, these exchanges happened on the picosecond time scale. Characteristic orientation of solvent molecules in the Na+ solvation shell can be studied from spatial distribution functions (SDFs). A three-dimensional description of the distribution of the Na+ ions around the solvent molecules is shown in Figure 4 (top panel). The isosurfaces close to the solvent molecules represent the probability to find the Na+ in a certain position, giving a clear representation of the mobility of the solvation cage. Interestingly, we found that each of them corresponds to a different picture shown in Figure 4. The results show that the ordering observed in these solutions. It can clearly show that the maximal of the Na+ distribution lies along the same directions as maximal of acceptor-hydrogen bonded water and methanol molecules. Since the oxygen atom of ethanol is not H-bonded perfectly along Na O bond, there is a broad configurational phase space of Na+. We also show a spatial distribution function, where the local coordinate system is fixed by the Na+ and oxygen atoms of water, methanol and ethanol molecules in the first hydration shell. This SDF is displayed in Figure 4 (bottom panel). The SDF shown in this figure clearly shows that the first solvation shell around Na+ contains radially and angularly different subshells. The results show that the positions of the water molecules in the first hydration shell become very well defined in their positions. The Na+ has a very well defined first hydration shell of water molecules in a
Table 1 First and second rank reorientational relaxation time for water, methanol and ethanol. Species
s1 (ps)
s2 (ps)
Water in the solvation shell of Na+ Methanol in the solvation shell of Na+ Ethanol in the solvation shell of Na+ Water in bulk Methanol in bulk Ethanol in bulk
7.69 10.42 10.11 8.64 7.81 7.63
6.85 7.83 7.46 7.19 5.48 4.41
octahedral symmetry. However, methanol and ethanol molecules in the first hydration shell are rather flexible. 3.2. Dynamics properties 3.2.1. Orientational relaxation of solvent molecules The time scales of reorientation of solvent molecules in these solution systems are crucial to insight its ability to solvate other species. Compared to water, there is limited experimental data on the reorientational dynamics of liquid alcohols. The orientational correlation function Cl(t) can provide valuable information on the reorientation dynamics for a molecular system. The rotational dynamics of water, methanol and ethanol molecules can be studied by calculating the orientational time correlation functions Cl(t) of lth order:
C l ðtÞ ¼
hP l ðv^ ð0Þ v^ ðtÞi hPl ðv^ ð0Þ v^ ð0Þi
Y. Zeng et al. / Chemical Physics Letters 538 (2012) 60–66
65
Figure 6. Nomalized translational velocity autocorrelation function of the Na+ ion in water, methanol and ethanol (left panel) and the corresponding power spectra of the velocity autocorrelation function of the Na+ ion in the water, methanol and ethanol (right panel).
^ ðtÞ denotes the where Pl is the lth order Legendre polynomial and v unit vector along the OH bond of a solvation molecule and Cl(t) is averaged over all the vector bonds that are present in the solution. The time correlation function for l = 1 is related to the infrared spectrum and that l = 2 relates to the depolarized Raman spectrum of these systems, which are usually measured by experiment [47,7]. The calculated C1(t) and C2(t) of the first solvation shells and bulk water, methanol and ethanol molecules are shown in Figure 5. The Cl(t) functions have been fitted to exponential function at long times. The corresponding orientation correlation time, based on an exponential fitting of the long-time decay of the orientational correlation functions and these values are given in Table 1. It is found from the figures that there is a fast inertial decay followed by a slower relaxation for Cl(t). An oscillation is also observed in the relaxation at short times due to librational contribution to the rotational correlation functions [48]. The short-time relaxation is faster for water and slower for the ethanol compared to that of methanol. The decay profiles of orientational correlations of the water molecules in the first solvation shell of Na+ ion show faster relaxation for both first rank and second rank correlations compared to bulk water. The interactions of the Na+ ion force the water molecules to be densely packed in its surroundings in the first solvation shell. Such compact arrangement restricts the orientational motion of the water molecules in the first solvation shell. However, methanol and ethanol has relatively bigger molecule size, the molecular size has a significant effect on the orientational relaxation. So the decay of orientational correlations for solvation shell molecules show somewhat slower relaxation (both first rank and second rank) compared to that of bulk methanol and ethanol. 3.2.2. Translational dynamics of Na+ The velocity autocorrelation functions (VACFs) for the Na+ in water, methanol and ethanol are plotted in Figure 6. Vibrational spectra are one of the main tools for investigating the dynamical properties of condensed systems. It can be conveniently evaluated from the Fourier transform of time correlation functions. Their corresponding power spectra have been calculated and shown in Figure 6.
In Figure 6, the VACF of Na+ is different in the three systems. Na+ in water has a faster initial decay and backscattering. Translational dynamics of Na+ in these systems is affected in a small degree by solvent molecular size. It results in a slightly deeper backscattering of the VACF function for water. It produces a significant difference on their relative relevance bands. In the Na+-methanol and Na+-ethanol systems, two relevant bands can be fitted appear at 150 and 250 cm1, which is in good agreement with the bands from 125 to 300 cm1 obtained experimentally in NaI methanolic solutions [49]. There is only one relevant band in Na+-water system. Their frequencies are in the range of the ones corresponding to the first bands observed in the ion’s spectra for methanol and ethanol. The lower frequency band can be assigned to the ion–solvent complex vibrations, whereas the higher frequency band corresponds to vibrational motions of the solute in the cage of its first solvation shells.
4. Conclusions In the present study, Car–Parrinello molecular dynamics of Na+ in liquid water, methanol and ethanol are reported. The results of structure of the solvation shell are in an agreement with previous some experimental data. The coordination number of the Na+ has been computed, showing that the first solvation shell is hex-coordinated as in water, whereas for methanol and ethanol the average of the coordination number is slightly lower than in water. Orientational decay profiles of the solvent molecules in the first solvation shell of Na+ ion is quite different from those of bulk molecules. The decay profiles of orientational correlations of the water molecules in the first solvation shell of Na+ show faster relaxation for both first rank and second rank correlations compared to bulk water. However, the decay of orientational correlations for solvation shell methanol and ethanol molecules showed somewhat slower relaxation (both first rank and second rank) compared to that of bulk for the sake of the bigger molecular size. Translational dynamics of Na+ in water, methanol and ethanol are affected in a small degree by molecular size. It results in a
66
Y. Zeng et al. / Chemical Physics Letters 538 (2012) 60–66
slightly deeper backscattering of the VACF for water. It produces a significant difference on their relative relevance bands. Their frequencies are in the range of the ones corresponding to the first bands observed in the ion’s spectra for methanol and ethanol. The lower frequency band can be assigned to the ion–solvent complex vibrations, whereas the higher frequency band corresponds to vibrational motions of the solute in the cage of its first solvation shells. Acknowledgements This work was supported by Natural Science Foundation of China (Grant No. 20806064). Generous allocations of computer time by the ScGrid plan of Supercomputing center of CAS. References [1] T.W.N. Bieze, R.H. Tromp, J.R.C. van der Maarel, M.H.J.M. van Strien, M.C. Bellissent-Funel, G.W. Neilson, J.C. Leyte, J. Phys. Chem. 98 (1994) 4454. [2] Y. Kameda, K. Sugawara, T. Usuki, O. Uemura, Bull. Chem. Soc. Jpn. 71 (1998) 2769. [3] T. Megyes, S. Balint, T. Grosz, T. Radnai, I. Bako, J. Chem. Phys. 128 (2008) 044501. [4] J.M. Alia, H.G.M. Edwards, J. Solution Chem. 29 (2000) 781. [5] F. Rull, C. Balarew, J.L. Alvarez, F. Sobron, A. Rodriguez, J. Raman Spectrosc. 25 (2005) 933. [6] A.W. Omta, M.F. Kropman, H.J. Bakker, Science 301 (2003) 347. [7] K.J. Tielrooij, C. Petersen, Y.L.A. Rezus, H.J. Bakker, Chem. Phys. Lett. 471 (2009) 71. [8] S. Park, M. Odelius, K.J. Gaffney, J. Phys. Chem. B 113 (2009) 7825. [9] D. Marx, A. Chandra, M.E. Tuckerman, Chem. Rev. 110 (2010) 2174. [10] A.K. Soper, K. Wecstom, Biophys. Chem. 124 (2006) 180. [11] L.X. Dang, J. Chem. Phys. 96 (1992) 6970. [12] M. Migliore, G. Corongiu, E. Clementi, G.C. Lie, J. Chem. Phys. 88 (1988) 7766. [13] K. Leung, S.B. Rempe, O.A. von Lilienfeld, J. Chem. Phys. 130 (2009) 204507. [14] M. Gaigeot, M. Sprik, J. Phys. Chem. 107 (2003) 10344. [15] Y. Liu, H. Lu, Y. Wu, T. Hu, Q. Li, J. Chem. Phys. 132 (2010) 124503.
[16] K. Hermansson, P.A. Bopp, D. Spangberg, L. Pejov, I. Bako, P.D. Mitev, Chem. Phys. Lett. 514 (2011) 1. [17] J. VandeVondele, F. Mohamed, M. Krack, J. Hutter, M. Sprik, M. Parrinello, J Chem. Phys. 122 (2005) 14515. [18] T. Todorova, A.P. Seitsonen, J. Hutter, I.-F.W. Kuo, C.J. Mundy, J. Phys. Chem. B 110 (2006) 3685. [19] D. Marx, M. Sprik, M. Parrinello, Chem. Phys. Lett. 273 (1997) 360. [20] M.M. Naor, K.V. Nostrand, C. Dellago, Chem. Phys. Lett. 369 (2003) 159. [21] L. Petit, R. Vuilleumier, P. Maldivi, C. Adamo, J. Chem. Theory Comput. 4 (2008) 1040. [22] E.S. Shamay, V. Buch, M. Parrinallo, G.L. Richmond, J. Am. Chem. Soc. 129 (2007) 12910. [23] H. Gunaydin, K.N. Houk, J. Am. Chem. Soc. 130 (2008) 15232. [24] A. Rodriguez-Fortea, M. Iannuzzi, M. Parrinello, J. Phys. Chem. B 110 (2006) 3477. [25] A.P. Lyubartsev, A. Laaksonen, Chem. Phys. Lett. 325 (2000) 15. [26] J.A. Morrone, M.E. Tuckerman, Chem. Phys. Lett. 370 (2003) 406. [27] J.W. Handgraaf, E.J. Meijer, J. Chem. Phys. 121 (2004) 10111. [28] M. Pagliai, G. Cardini, V. Schettino, J. Phys. Chem. B 109 (2005) 7475. [29] C. Fafalli, M. Pagliai, G. Cardini, V. Schettino, Theor. Chem. Acc. 118 (2007) 417. [30] M. Seo, Y. Yun, J. Lee, Y. Tak, J. Power Sources 159 (2006) 59. [31] Y. Yang, Y.C. Liang, J. Power Sources 165 (2007) 185. [32] CPMD V 3.11; copyright IBM Corp. 1990–2007, copyright MPI fur Festkorperforschung Stuttgart 1997–2001. [33] N. Troullier, J.L. Martins, Phys. Rev. B 43 (1991) 1993. [34] L. Kleinman, D.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. [35] R. Vuilleumier, M. Sprik, J. Chem. Phys. 115 (2001) 3454. [36] J. Handgraaf, T.S. van Erp, E.J. Meijer, Chem. Phys. Lett. 367 (2003) 617. [37] S. Goedecker, M. Teter, J. Hutter, Phys. Rev. B 54 (1996) 1703. [38] M. Sprik, J. Hutter, M. Parrinello, J. Chem. Phys. 105 (1996) 1142. [39] K. Laasonen, M. Sprik, M. Parrinello, R. Car, J. Chem. Phys. 99 (1993) 9080. [40] W.L. Jorgensen, D.S. Maxwell, J. Tirado-Rives, J. Am. Chem. Soc. 118 (1996) 11225. [41] R. Caminiti, G. Licheri, G. Paschina, G. Piccaluga, G. Pinna, J. Chem. Phys. 72 (1980) 4522. [42] S.B. Rempe, L. Pratt, Fluid Phase Equilib. 183–184 (2001) 121. [43] J.A. White, E. Schwegler, G. Galli, F. Gygi, J. Chem. Phys. 113 (2000) 4668. [44] W.L. Jorgensen, B. Bigot, J. Chandrasekhar, J. Am. Chem. Soc. 104 (1982) 4584. [45] I. Dauster, M.A. Suhm, U. Buck, T. Zeuch, Phys. Chem. Chem. Phys. 10 (2008) 83. [46] R.M. Forck, I. Dauster, U. Buck, T. Zeuch, J. Phys. Chem. A 115 (2011) 6068. [47] S. Park, M.D. Fayer, Proc. Natl. Acad. Sci. 104 (2007) 16731. [48] L.H. de la Pena, P.G. Kusalik, J. Chem. Phys. 121 (2004) 5992. [49] B. Guillot, P. Marteau, J. Obriot, J. Chem. Phys. 93 (1990) 6148.