Molecular Dynamics Study of the Structure and Dynamics of the

The observed order of residence times for the monovalent ions, τ(Li+) > τ(Na+) ...... According to Degrève and Quintal,12 solvent molecules inside ...
4 downloads 0 Views 445KB Size
Download PDF
J. Phys. Chem. 1996, 100, 15677-15687

15677

Molecular Dynamics Study of the Structure and Dynamics of the Hydration Shell of Alkaline and Alkaline-Earth Metal Cations Stefan Obst* and Hans Bradaczek Institut fu¨ r Kristallographie, Freie UniVersita¨ t Berlin, Takustrasse 6, 14195 Berlin, Germany ReceiVed: May 14, 1996X

Molecular dynamics simulations of hydrated alkaline and alkaline-earth metal cations at room temperature (T ) 300 K) were carried out using the CHARMM22 force field. Dynamic and static properties of systems containing one ion and 123 or 525 water molecules were investigated by analysis of trajectories of 1 ns duration and compared to experimental and theoretical results. In addition, the size and the direction of the elementary motions of both the ions and the water molecules were investigated on the scale of the integration time step of 1 fs. Comparison between systems of different size revealed that for the larger system the diffusion coefficient and the number of hydrogen bonds were increased. Radial pair distribution functions and coordination numbers are in good agreement with X-ray and neutron scattering data. The diffusion coefficient D of bulk TIP3P water in a system with 528 water molecules was by one-fourth higher than the experimental value. Minor differences of approximately 10% between experimental and simulated diffusion coefficients were found for Li+, Na+, K+, and Mg2+. On the other hand, D was underestimated by the simulation for Ca2+ and Sr2+ by as much as 30%. On the average, 2.9 hydrogen bonds per bulk water molecule were found. The observed order of residence times for the monovalent ions, τ(Li+) > τ(Na+) > τ(K+), is in good agreement with the literature. Although τ was expected to increase with decreasing mass of the ion, the exchange of water molecules from the solvation shell of Mg2+ occurred much faster than for Ca2+.

1. Introduction It is well-known that the simple model of ionic diffusion, known as the Stokes-Einstein relation1

D)

kBT 6πηrsol

(1)

with viscosity η and solute radius rsol, which relates the selfdiffusion coefficient D to the reciprocal of rsol, does not always strictly hold true, although it describes the diffusive motion of many liquids well.2 Especially for simple inorganic cations such as the mono- and divalent ions of the alkaline and alkalineearth metals, respectively, this has been shown by both experiments3 and molecular dynamics simulations.4,5 A theoretical attempt to improve the understanding of ionic mobility is the so-called “solventberg” concept.5 This model assumes that ions are strongly hydrated by the water molecules of the first hydration shell resulting in a comigration of the hydration shell with the central ion. This leads to an increased effective radius of the ion with the water molecules at fixed positions within the first hydration shell. In contrast to the static picture of the first hydration shell emerging from the solventberg concept, it is known from NMR measurements3 and computer simulations4,5 that there is motion of the water molecules within the first hydration shell and exchange with molecules of the bulk, too. For the alkaline ions, it has been shown that the binding strength of water molecules to the central ion increases with decreasing radius.4,5 Moreover, the reproduction of some of the well-known features of solutions of alkali and alkaline-earth metal ions should be a good benchmark for the ability of a molecular dynamics force field to produce reliable results. Such calcula* Corresponding author: e-mail, [email protected]; phone, +49 30 838 3457; FAX, +49 30 838 3464. X Abstract published in AdVance ACS Abstracts, September 1, 1996.

S0022-3654(96)01384-6 CCC: $12.00

tions have been made for different water models; however, to our best knowledge, no such data have been published up to now for the CHARMM22 force field. Recently, it was noted by Sansom et al.6 that the TIP3P water model used in CHARMM overestimates the diffusion coefficient by as much as 50%-75% at T ) 300 K as compared to the experimental value. In the present paper, we have calculated static and dynamic properties of alkaline and alkaline-earth metal ions in water using the CHARMM22 force field and compared them to results derived from both experiments and simulations. The static properties include radial distribution function g(r), coordination number ncoor, and number of hydrogen bonds n(HB); the dynamic properties investigated are diffusion coefficient D, hydration number nshell, and mean residence times τ. In addition, we analyzed the fluctuations of the hydration shell at the scale of the integration time step of 1 fs in order to improve our knowledge on the molecular basis of the different mobility of the investigated ions and their hydration shells. For the investigation of the dynamic properties of the hydrated ions, extended simulation times were necessary. Initially, small systems with 123 water molecules were used; however, the obvious size dependence of some system properties obliged us to extend the system size to 525 water molecules. As a consequence, the computing time becomes fairly long even on modern workstations, thus limiting the trajectories to 1 ns. The rest of the paper has the following structure. Section 2 contains some technical details of the simulations and of the force field used. In section 3, our results will be presented and discussed. In the last section (section 4), a summary of our work is given. 2. Computational Methods Molecular dynamics (MD) simulations were performed on SGI IRIS Indigo workstations (Silicon Graphics) using the program CHARMM7 with the CHARMM22 force field as © 1996 American Chemical Society

15678 J. Phys. Chem., Vol. 100, No. 39, 1996

Obst and Bradaczek

TABLE 1: Parameters Used for Setting up the Different Simulationsa simulation +

Li Na+ K+ Mg2+ Ca2+-123 Ca2+-525 Sr2+ H2O-125 H2O-528

n(H2O)

n(ion)

box size (Å)

t (ps)

525 525 525 525 123 525 525 125 528

1 1 1 1 1 1 1 0 0

25.08 25.08 25.08 25.08 15.50 25.08 25.08 15.50 25.08

1020.80 1021.30 1020.90 1023.90 1209.60 1024.50 1027.10 120.00 120.00

a n(H O): number of TIP3P water molecules in the system. n(ion): 2 number of ions in the system. Box size: size of the cubic simulation box. t: total simulation time, including heating (10 ps) and equilibration (10 ps) phases.

provided by MSI (Molecular Simulations Inc., Waltham, MA). Two different system sizes were used: 125 molecules of water in a cubic box of 15.50 Å length (small system) and 528 molecules of water in a box of 25.08 Å side length (large system). The density was 1.00 g/cm3. Water was represented by the TIP3P model,8 and standard parameters from CHARMM22 were used for the ions Li+, Na+, K+, Mg2+, Ca2+, and Sr2+. Periodic boundary conditions9 in x, y, and z directions with 26 images were used. The time step for integration of the leapfrogalgorithm was 1 fs. Energies, atomic coordinates, and velocities were saved every 100th time step (0.1 ps). Bonds containing hydrogen atoms were constrained by the SHAKE algorithm.10 Nonbonded interactions were cut off at rcut ) 12 Å by multiplying the energy function with the following shifting term:7

(

1-

)

2rij2

rij4

rcut

rcut4

+ 2

(2)

The electrostatic interactions were scaled by  ) 1, because the solvent was explicitly modeled. The water boxes were energy minimized using 1000 steps of steepest descent minimization followed by 4000 steps of adopted basis Newton Raphson minimization.7 The ion was located at the center of the simulation box, and all water molecules within a sphere of radius r ) 2.6 Å around the ion were deleted. The resulting ensemble was minimized by following the same protocol as was used for pure water. Next, an initial 10 ps phase of heating from 0 to 300 K and an equilibration phase of 10 ps at 300 K followed. The simulations of pure water were continued for another 100 ps, whereas the ionic solutions were simulated for at least another 1000 ps. Ca2+ was simulated for 1 ns both in the large and in the small box. Details are shown in Table 1. 3. Results and Discussion 3.1. Energies and Temperatures. In Table 2, information about energies and temperatures during the simulations is given.

Figure 1. Temperature dependence of the diffusion coefficient D of water. D is given in 10-5 cm2/s. Data taken from Weinga¨rtner11 (dots, experimental results; line, extrapolation).

All simulations reached the designated temperature of 300 K with the exception of the small system with Ca2+ (Ca-123), where the average temperature was only 290 K. The reason for this deviation remains unclear, as the same simulation setup has been used for all simulations of both large and small systems. Nevertheless, comparing results between Ca-123 and Ca-525, one should keep in mind the difference in mean simulation temperatures because dynamic properties such as diffusion coefficients D are known to be temperature dependent. For example, from extrapolation of experimental values (cf. Figure 1) given by Weinga¨rtner,11 it can be assumed that for water D(290 K) amounts to only about 80% of D(300 K). The total energy was well conserved for all simulations as indicated by the small value of the standard deviation (SD). The fluctuations in total energy for pure water were smaller than those of the ionic solutions by more than 1 order of magnitude. At first glance, this seems to indicate problems arising from the fact that the ionic charge of the cation was not counterbalanced. However, if this was true, one would expect that such effects are more pronounced for systems containing divalent cations, and since this was not observed, the reason for the high energy fluctuations remains unclear. As expected, the magnitudes of the fluctuations of the potential and the kinetic energies are comparable. 3.2. Radial Pair Distribution Functions and Coordination Numbers. The radial distribution functions (rdf’s) g(r) between the oxygen atoms of the water molecules and the ions were averaged from the first 100 ps after equilibration of each simulation. For pure water, the rdf was calculated between the oxygen atoms of all water molecules. The maximum distance rmax ) lb/2 for evaluation of the rdf is determined by the simulation box length lb due to the minimum image convention.9 The radius of the first hydration shell was defined by the first minimum of the rdf. The coordination numbers ncoor of the ions and of pure water were obtained from the rdf g(r)4πr2 dr by

TABLE 2: Simulation Times and Energies for the MD Simulationsa +

Li Na+ K+ Mg2+ Ca2+-123 Ca2+-528 Sr2+ H2O-125 H2O-528

∆t (ps)

T ( SD (K)

Etot ( SD (kJ/mol)

Ekin ( SD (kJ/mol)

Epot ( SD (kJ/mol)

1001.1 1001.3 1000.9 1003.8 1189.6 1004.5 1007.1 100.0 100.0

299.96 ( 6.25 299.99 ( 6.25 300.03 ( 6.11 299.98 ( 6.19 289.51 ( 12.53 299.99 ( 6.10 300.09 ( 6.23 303.72 ( 12.61 299.96 ( 6.12

-21011.7 ( 29.9 -20854.1 ( 33.5 -20780.4 ( 29.7 -22072.3 ( 27.8 -6121.5 ( 6.0 -21893.6 ( 31.2 -21718.2 ( 34.1 -4892.2 ( 0.1 -20670.5 ( 0.5

3930.9 ( 81.7 3928.5 ( 81.9 3928.6 ( 80.1 3931.2 ( 81.0 888.2 ( 38.5 3928.6 ( 79.9 3932.5 ( 81.4 943.2 ( 39.1 3946.8 ( 80.5

-24929.5 ( 86.5 -24775.2 ( 87.7 -24702.2 ( 85.3 -25990.4 ( 85.9 -7009.7 ( 38.7 -25814.3 ( 85.8 -25637.6 ( 88.2 -5835.4 ( 39.2 -24617.2 ( 80.5

a All values given are the mean ( standard deviation (SD) and are representative of the simulation phase of duration ∆t. ∆t: simulation time (ps) used for data analysis excluding heating and equilibration phases. T: average temperature (K) during the simulation time ∆t. Etot (Ekin, Epot): total (kinetic, potential) energy in kJ/mol, averaged over the simulation time ∆t.

MD Simulation of Hydrated Cations

J. Phys. Chem., Vol. 100, No. 39, 1996 15679

Figure 2. Comparison of the radial distribution function gion-oxygen(r) of monovalent (left: solid, Li+; dot, Na+; dash dot, K+) and divalent ions (right: solid, Mg2+; dot, Ca2+; dash dot, Sr2+).

integrating from 0 up to the first minimum of g(r).9 ncoor represents a static property which should not be confused with the dynamic hydration number nhydr introduced by Impey et al.,4 which is given by the ratio of residence times of bulk water and hydration water of the ions. Degre`ve and Quintal12 argued that the definition of ncoor based on the integration of the rdf is close to most of the experimental data. No dependency of the rdf on the system size was observed neither for H2O nor for Ca2+ when comparing rdf data for distances from 0 up to one-half the simulation box’s size. This is in accordance to results published by Andrea et al.,13 who found that distribution functions of water are insensitive to the size of the system, provided the cube’s size is equal to or larger than twice the range of the potential, using system sizes of 64 and 512 molecules and potentials truncated at 6 and 12 Å, respectively. On the other hand, Brooks14 found large effects in ion-water distribution functions for truncation distances between 7 and 12 Å. Our simulations were carried out using a 12 Å cut-off radius, thus minimizing artifacts induced by the truncation scheme. For Ca2+ and Sr2+, up to three hydration shells can be identified as peaks of the rdf, whereas for the other ions and water, only two distinct hydration shells were observed. The height of the first peak of the rdf varies with the charge and the size of the ion simulated, with values ranging from 4.3 to 12.3 for the monovalent ions K+ and Li+ and from 10.3 to 17.2 for the divalent ions Sr2+ and Mg2+, respectively (cf. Figure 2). These results are in good agreement with trends observed by Probst et al.15 for Mg2+ and Ca2+ using a central force water model. For water, a peak value of 3.2 was observed. Although the exact values of peak height are not very meaningful since they are extremely sensitive to the interval of sampling distance data r and to smoothing of the curve (cf. Figure 3), the following comparison can be made: The height of second peaks of the rdf decreases in the following order: Mg2+ ≈ Ca2+ ≈ Sr2+ ≈ Li+ > Na+ ≈ K+ (cf. Figure 4). In general, in experimental studies, the hydration shell structure was found to be weakened and to be less well defined with increasing ionic radius for both the alkaline and the alkaline-earth metal cations.16 In a simulation study, Probst et al.15 noticed a somewhat better defined second hydration shell for Ca2+ than for Mg2+. Our observation of a highly ordered hydration shell for Sr2+ is inconsistent with the results of Neilson and Broadbent.17 They concluded from neutron diffraction experiments that Sr2+ is only weakly coordinated by water molecules and shows a less defined hydration shell structure than Mg2+ and Ca2+. On the other hand, using X-ray diffraction, Caminiti et al.18 earlier found a “fairly ordered” second hydration shell for Sr2+. Due to technical limitations, most of

Figure 3. Dependence of apparent height of peaks of the rdf on the interval used for data sampling and on subsequent data processing. Smoothing procedures applied to the data result in a decrease of the peak height which is comparable to a reduction of the sampling resolution. On the abscissa, the resolution at which data were sampled is given in Å: dots, harmonic smoothing, ys(i) ) 1/4(y(i - 1) + 2y(i) + y(i + 1)); squares, different sampling intervals; diamonds, sliding i+j average, yj(i) ) (1/(2j + 1))∑k)i-j y(k) with 2j + 1 ) number of data points used in the averaging procedure. In this example, data from the first maximum of the rdf of Mg2+ were used.

the experimental studies cited here were carried out using solutions with minimum concentration of 1 mol/L. At a concentration of 1 mol/L, less than 55 water molecules surround one cation, if the hydration of the anion is taken into account. At this molar ratio, more than two hydration shells can hardly be expected to be observed. On the other hand, the concentration in the simulated large system is approximately 0.1 mol/L (525 water molecules per ion), which might enable the ions to extend their hydration shell structure farther into the water. The second peak is in general by 1 order of magnitude less pronounced than the first peak, and it is of almost equal height for all simulations, i.e., for pure water as well as for the ionic solutions. Our results for TIP3P water (cf. Table 3) are in good agreement with the values given by Lee and Rasaiah5 for TIP4P water. For Li+ and Na+, the r value of the first minimum of g(r) is by 0.3 Å smaller than the values given by Lee and Rasaiah, whereas for K+, our result is by 0.1 Å higher. Taking into account the difficulty of determining the exact position of the minimum due to the shallowness of the rdf curve, the observed difference is tolerable. Moreover, Lee and Rasaiah used the center-of-mass of the water molecules instead of the oxygen atom as a reference point for the calculation of the rdf. Our values for the position of the first maximum, rio(Li+) ) 2.0 Å, rio(Na+) ) 2.5 Å, and rio(K+) ) 2.9 Å, are in good agreement with experimental results of 1.95, 2.4, and 2.7 Å compiled by Enderby19 and with the values given by Impey et

15680 J. Phys. Chem., Vol. 100, No. 39, 1996

Obst and Bradaczek

Figure 4. Number nHB(r) of hydrogen bonds accepted and rdf g(r) versus distance r from the central ion for monovalent ions (left) and divalent ions (right). Left: Li+ (upper), Na+ (middle), K+ (lower). Right: Mg2+ (upper), Ca2+ (middle), Sr2+ (lower). Solid: nHB(r). Dot: g(r). Notice the different scale factors on the right scale (g(r)) of each figure.

TABLE 3: Minima and Maxima of Radial Distribution Functions g(r) of Ionic Solutions and Pure Water at 300 Ka this work first max ion Li+

Na+ K+ Mg2+ Ca2+ Sr2+ H2O

rio (Å) 2.0 2.5 2.9 2.2 2.5 2.9 2.8

literature first min

gio 12.3 6.7 4.3 17.2 13.3 10.3 3.1

rio (Å)

first max gio

2.8 3.2 3.8 3.1 3.5 3.7 3.4

0.0 0.2 0.5 0.0 0.0 0.0 0.9

roo (Å) 2.2(IX);

1.95(VI)

2.5(IX); 2.4(VI) 2.7(IX); 2.7(VI) 2.1(II);2.1(VIII) 2.4(II); 2.4(IV) 2.6(II); 2.6(VII) 2.9(I); 2.8(V)

al.4 of 1.98 , 2.29, and 2.76 Å for Li+, Na+, and K+, respectively. Other studies using various water models yielded comparable results (Petitt and Rossky20 and references therein). For the alkaline-earth cations Mg2+, Ca2+, and Sr2+, we determined the position of the first minimum to be at rio ) 2.2 Å, rio ) 2.5 Å, and rio ) 2.9 Å, respectively, whereas Neilson and Enderby16,21 give experimental values of rio ) 2.06, rio )

24

9.6 7.0 5.6 19.2(III) 14.0(IV) 3.1(I); 2.8(V)

roo (Å)

goo

3.1 3.5 3.7

0.03 0.22 0.35

3.2-3.7(IV)

0.1(IV)

3.3(I); 3.3(V)

0.7(I); 0.9(V)

rio: radius between water oxygen o and ion i (I): Soper and Phillips. (II): Data compiled by Neilson and Enderby. (III): Probst et al.15 (IV): Dietz et al.45 (V): Jorgensen et al.8 (VI): Enderby.19 (VII): Neilson and Broadbent.17 (VIII): Skipper et al.;23 Powell et al.46 (IX): Lee and Rasaiah.5 (I, II, VI, VII, VIII): Experimental data. (III, IV, V, IX): Simulation. a

55

first min goo

21

2.46, and rio ) 2.64 Å. Hewish et al.22 obtained rio ) 2.4 Å for the first peak of the rdf of Ca2+ from neutron diffraction. The position of the second peak of the rdf of Mg2+ agrees well with the result obtained from neutron scattering experiments of aqueous solutions of NiCl2. In terms of hydration properties, Ni2+(aq) and Mg2+(aq) can be regarded as structurally isomorphic species.23

MD Simulation of Hydrated Cations

J. Phys. Chem., Vol. 100, No. 39, 1996 15681

TABLE 4: Coordination Number ncoor of Different Ions and Water and Their Ionic Radiia r (Å) ncoor ncoor(lit.) ncoor(lit.) ncoor(lit.) ncoor(exptl) ncoor(exptl) a

Li+

Na+

K+

0.73 4.5 6(I) 5.3(II) 3.9(XI) 4-6(IX) 5.5(IV)

1.16 6.0 6.6(I) 6.0(II) 7.0(III) 4-6(IX) 5(XV)

1.65 7.6 8.0 7.6(VII) 5.5(XVI)

Mg2+

Ca2+

Sr2+

0.86 6.2

1.26 8.0

1.32 9.1

6(XIII) 6-8(IX) 6(XVII)

9(XII) 6-8(IX) 6-10(VIII)

6-8(IX) 10(XVIII)

H2O 1.57(X); 1.44(XIV) 4.9 4.9(VI) 5.0(III) 4.4(V)

radii47,48

r: crystallographic ionic for the coordination number given under ncoor. ncoor: coordination number, determined by integration of the rdf up to the first minimum. ncoor(lit.): coordination numbers from literature determined by simulation. ncoor(exptl): experimentally determined coordination numbers. (I): Lee and Rasaiah,5 Tip4P-water. (II): Impey et al.,4 MCY-water. (III): Garcia and Stiller,27 TIPS3-water. (IV): Newsome et al.,41 exptl. (V): Narten et al.,25 exptl. (VI): Svishchev and Kusalik,49 SPC/E-water. (VII): Forester et al.,50 SPC-water. (VIII): Hewish et al.,22 exptl (ncoor was found to be concentration dependent, increasing from 6 to 10 as the molality decreases from 4.5 to 1). (IX): exptl., cited from Neilson and Enderby.16,21 (X): r(H2O) calculated from V(H2O) ) 16.18 Å3.51 (XI): Degre`ve and Quintale,12 MC simulation, TIP4P-water. (XII): Probst et al.,15 CF water model. (XIII): Dietz et al., 45 CF water model. (XIV): Agmon.52 (XV): Skipper and Neilson,53 exptl. (XVI): Neilson and Skipper, 54 exptl. (XVII): Skipper et al.,23 exptl. (XVIII): Caminiti et al.,18 exptl. An extensive review on solvation numbers of ions is given by Hinton and Amis.55

Figure 5. Radial distribution function g(r) (solid) and coordination number (dotted) of TIP3P water.

In Table 4 it can be seen, that the coordination number ncoor increases with increasing radii of monovalent and divalent ions. For water, ncoor ) 4.9 was obtained. This value is slightly higher than the experimental value of ncoor ) 4.5 given by Soper and Phillips24 or ncoor ) 4.4 given by Narten et al.25 and Gaballa and Neilson26 but is close to the value obtained by Garcia and Stiller,27 who used a TIPS3 water model. The difference between these results might again be explained in part by the difficulties arising from the shallowness of the rdf (cf. Figure 5). Furthermore, Sciortino et al.28 showed for ST2-water a strong dependence of the number of water molecules in the first hydration shell on the density of the solution and, to a lower extent, on the simulation temperature. Our simulations were carried out at 300 K, whereas Sciortino et al. conducted simulations at 235 and 273 K with n(HB) ) 4.6 at F ) 1.0 g/cm3. For the monovalent ions, we found that ncoor ranges from 4.5 for Li+ to 7.6 for K+, whereas for the divalent ions, we measured values between 6.2 for Mg2+ and 9.1 for Sr2+. Our values of ncoor for the ions are in good agreement with the results reported in literature (cf. Table 4). Only in the case of Ca2+, is the coordination number ncoor ) 8 determined by our MD simulation smaller than the value ncoor ) 10 measured by neutron diffraction for a solution of molality 1 and smaller than ncoor ) 9, obtained by Probst et al.15 and by Heinzinger29 by MD simulation using a central force water model. Unfortunately, to our knowledge, no simulation data using the CHARMM force field for the divalent alkaline-earth metals have been published to date. The increase of ncoor with increasing ionic radius might in part be explained by geometrical considerations, i.e. the larger the ion the greater the available space for hydrating water molecules in the first sphere. Comparing mono- and divalent ions, the

Figure 6. Fraction p of water molecules participating in a given number of hydrogen bonds n(HB) for pure water and for Ca2+: black dots, H2O-125; diamonds, H2O-528; black squares, SPC-water (Kowall and Geiger33); open squares, Ca-525. The lines connecting the data points are for clarity only and have no physical meaning.

higher charge density of the divalent ions might be the reason for a higher coordination number and a stronger hydration (as can be seen by the increase in residence times τ when going from the mono- to the divalent ions, cf. section 3.5) as supposed by Balbuena et al.30 3.3. Hydrogen Bonds. For the analysis of hydrogen bonds between the water molecules, the following common geometric definition of a hydrogen bond was used: Two water molecules are called hydrogen bonded, if the distance between their two oxygen atoms is less than 3.4 Å and the O-H‚‚‚O angle is larger than 135°. Similar definitions were used by other authors in the past.27,31,32,33 The number of hydrogen bonds n(HB) per water molecule was analyzed using the first 100 ps of each simulation after equilibration. n(HB) is given as the time average of the ratio of the number of hydrogen bonds both donated and accepted to the number of water molecules found at the same distance from the ion or a randomly chosen water molecule. In our terminology, the hydrogen bond acceptor is the water oxygen which is not covalently bound to the hydrogen atom participating in the hydrogen bond. To check for a time dependence, the 1 ns trajectory of Ca525 was divided into 10 fragments of 100 ps duration each which were analyzed separately. Comparison between n(HB) calculated from the different fragments yielded equivalent results of n(HB) ) 2.9. There was a significant difference in n(HB) between the large and the small system of pure water. In Figure 6, the population of water molecules participating in a given total number of hydrogen bonds n(HB) (both accepted and donated) is shown

15682 J. Phys. Chem., Vol. 100, No. 39, 1996

Figure 7. Distance dependence of nHB(r) for the small (H2O-125, dot) and the large systems (H2O-528, dash) of pure water compared to its rdf g(r) (solid). nHB(r): number of hydrogen bonds accepted per water molecule. The distance r from a central water molecule is given in Å.

for both systems. In the smaller system, more water molecules participate in zero, one, or two hydrogen bonds and less have three, four, or five hydrogen bonds than in the larger system. A number of five hydrogen bonds per water molecule can be explained in terms of bifurcated bonds, which are know to occur in hydrate crystals34 and which were also found in MD simulations.28,35 Kowall and Geiger33 found a distribution of n(HB) which was shifted to higher n(HB) as compared to our results. They give a value n(HB) ) 3.34 for bulk SPC-water. Depending on the lifetime of the hydrogen bonds investigated, Sciortino and Fornili36 found values of n(HB) between 0.5 and 4.6 for minimum lifetimes between 3.95 and 0 ps, respectively. The average number of hydrogen bonds in our simulation is n(HB) ) 2.6 for the small system versus n(HB) ) 2.9 for the large system, corresponding to hydrogen bond lifetimes of 0.35 and 0.75 ps in Sciortino’s work. Depending on the energy cutoff used for hydrogen bond assignment, Geiger et al.37 found by MD simulation n(HB) ) 2.4 and n(HB) ) 3.1 for VHB ) 2.1 × 10-20 J and VHB ) 2.5 × 10-20 J, respectively. The same effect of the system size on n(HB) as for pure water was observed in the two systems containing Ca2+. For the small system (Ca-125), n(HB) ) 2.5 was obtained, vs n(HB) ) 2.9 for the large system (Ca-525). No difference in the average n(HB) was detected for water in the large systems, no matter which ion it contained. In other words, the water molecules in the large ionic solutions show a bulklike n(HB). This might be due to the increasing amount of water molecules which are to far away from the ion to be affected by its water-structuring properties as compared to the small systems. These water molecules dominate the n(HB) statistics over the small number of water molecules in the first hydration shell with a distorted hydrogen bond pattern. Thus the chosen system size with more than 500 water molecules seems to ensure the correct representation of these properties. To improve our knowledge of the effects of the mono- and divalent ions and of water on the hydration water, we investigated the distance dependence of n(HB), nHB(r). r is the distance between the water oxygen atom and the central ion or a central water molecule. One would expect to observe a strong orientational order imposed on the water of molecules of at least the first hydration shell by the central cation,38 resulting in a first hydration shell which is mainly a donator but not an acceptor of hydrogen bonds, because the negatively charged water oxygen atoms will be pointing toward the cation. Therefore, we decided to omit the contribution of the donated hydrogen bonds and focused on accepted hydrogen bonds only. In Figure 7, nHB(r) is plotted for both the small and the large systems of pure water together with the rdf. Obviously, there

Obst and Bradaczek is a significant difference between H2O-125 and H2O-528 in the number hydrogen bonds. On the other hand, it can be clearly seen that these curves have similar structures. Furthermore, the structure of nHB(r) follows the undulation of the rdf g(r). The rise of nHB(r) slightly precedes the rise of g(r) and reaches its maximum at the same distance where the rdf has its first peak. At the very radius rmin where the rdf reaches its first minimum, nHB(r) has a local minimum. Interestingly, the undulation of nHB(r) and g(r) does not run parallel for r > rmin. The difference in nHB(r) of approximately 20% between H2O-125 and H2O528 must be the outcome of the different size of the simulation box, thus showing the limits of the periodic boundary approach. The undulation of nHB(r) can be interpreted as follows. From the rdf, two hydration shells can be identified. Assuming that the members of each hydration shell will preferentially adopt a structure which enables them to participate in a maximum number of hydrogen bonds, one might expect that in these layers specific hydrogen bond patterns evolve. Every water molecule which travels between these two shells will loose some of its hydrogen bonds during its voyage because the hydrogen bond pattern is optimized for the molecules in the first and second hydration shell but not for those molecules which are between them. Therefore, the minimum in nHB(r) is correlated to rmin of the rdf. This hypothesis is supported by the finding of Degre`ve and Quintal12 that water molecules in the first hydration shell of Li+ ions preferentially adopt configurations which enable them to form hydrogen bonds with molecules outside the shell. Going from the second to the third hydration shell, the influence of the ion (or the central water molecule) vanishes and bulklike properties evolve, which leads to a plateau value of nHB(r). Vaisman et al.39 noticed a reduced number of hydrogen bonds in the first layer as compared to other layers further away from the solute when simulating an ammonium ion solvated in 499 molecules of SPC-water. The effect was less pronounced for methane as solute. The average number of hydrogen bonds is 3.25 and 3.29 for the first and the other shells in the case of methane and 2.67 and 3.02 hydrogen bonds for the solvated ammonia ion, respectively. Jorgensen et al.8 reported a value of n(HB) ) 3.50 for the TIP3P-water model. When comparing these values with our results, one should keep in mind the difference between counting only accepted hydrogen bonds and counting every hydrogen bond a water molecule participates in. It is not clear whether this is the reason for the observation of Vaisman et al. that n(HB) is lower in the first hydration shell than in layers farther away from the solute. Within the ionic solutions, the relative number of hydrogen bonds accepted, nHB(r), is zero from the center of the ion up to the beginning of the first hydration shell, which is indicated by the rise of the rdf. This is followed by a steep ascent to a maximum value of nHB(r) ≈ 1.6 which is reached at a distance which matches the second peak of the hydration shell. With increasing distance from the center of the ion, a slight decrease in the value of nHB(r) was observed after the first maximum. For Ca-123, this decrease of nHB(r) was much stronger than for the large system Ca-525 (cf. Figure 8). This might be explained by the limited size of the simulation box of (15.5 Å),3 although one would expect serious problems only for distances exceeding rmax ) 7.75 Å, with 2rmax ) box length lb. For the monovalent ions Li+ and K+, nHB(r) rises steeply to n(HB) ) 1.5 and remains stable at greater distances. For Na+, a maximum value of nHB(r) ) 1.7 and a decrease at distances of >5 Å were found. For the divalent ions, the hydrogen bond profile nHB(r) showed a steep ascent to nHB(r) ) 1.7 with a plateau value of

MD Simulation of Hydrated Cations

J. Phys. Chem., Vol. 100, No. 39, 1996 15683

Figure 8. The number nHB(r) of intra-water hydrogen bonds accepted depends on the system size: solid, Ca-525; dash, Ca-123.

nHB(r) ) 1.6. No significant differences between Mg2+, Ca2+, and Sr2+ were found. Compared to the monovalent ions, the ascent was at a greater distance and it was steeper. Only Na+ showed a behavior which was similar to that of the divalent ions. For the monovalent ions, the onset of nHB(r) correlated with the ionic radius, whereas for all of the divalent ions, the rise started at the same distance r. Although the shape of nHB(r) for the ionic solutions resembles nHB(r) for water, an important difference in the peak position must be noted. For water, the maximum values of nHB(r) coincide with the first peak of the rdf g(r), whereas for the ionic solutions, they appear at r values which correlate with the second peak of g(r), where r ) r2nd max. Taking into account that nHB(r) gives the number of hydrogen bonds accepted by a water molecule at distance r from the central molecule (water or cation), it is obvious that any water molecule in the first hydration shell of the cations can accept hydrogen bonds only from those members of the second hydration shell which are oriented with their hydrogen atoms toward the ion. The stronger the hydration of the ion (i.e., the farther the hydration shell extends into the water), the less probable such a disorder of a water molecule will occur. Thus, one would expect nHB(r) ) 0 for r < r2nd max, especially for the divalent ions. The small but nonzero number of hydrogen bonds which, in fact, were observed for r < r2nd max might be explained by the formation of hydrogen bonds within the first hydration shell and by water molecules exchanging between second and first hydration shell. On the other hand, in the case of pure water, molecules of the first hydration shell can accept hydrogen bonds from the hydrogen atoms of the central water molecule. Therefore, nHB(r) reaches its peak within the first hydration shell. 3.4. Diffusion Coefficient D. The diffusion constant was obtained from the simulations using the well-known Einstein relation:1

〈[zi(t) - zi(0)]2〉 tf∞ 2t

D ) lim

(3)

where zi(t) is the z coordinate of particle i at time t and zi(0) is the z coordinate of particle i at time t ) 0. Mean square displacements for each spatial direction x, y, and z were analyzed in segments of 5 ps to improve the statistical accuracy of the results and averaged afterward to give the mean value and its standard deviation (SD). In the following, all D values are given in 10-5 cm2/s. Time dependency was observed neither for water nor for Ca2+ when comparing the values of D obtained for subsequent segments of 100 ps taken from the 1 ns simulation (cf. Figure 9).

Figure 9. Diffusion coefficient D of water and Ca2+ is independent of the simulation time: gray, H2O; black, Ca2+. Error bars give the standard deviation from the average of 20 D values within the 100 ps time interval. Time interval indicates the 100 ps fragment evaluated from the simulation Ca-525.

D is larger for the monovalent ions than for the divalent ions. The diffusion coefficient increases with increasing molar mass of the monovalent ions, whereas the D values of the divalent ions appeared to be almost independent of their molar mass. This finding is in good agreement with experimental results for the diffusion of alkaline-earth metal ions compiled by Hertz.3 Water has the largest D value (cf. Table 5). The diffusion rate of water in the small system of pure water D(H2O-125) ) 2.74 is by 5% lower than the value for the large system, D(H2O528) ) 2.87. This finding is in accordance with the observation of Svishev and Kusalik,40 who simulated two SPC/E-water systems with 108 and 256 molecules with D108(H2O) ) 2.15 and D256(H2O) ) 2.24. In the large system Ca-525, the water molecules outside the first hydration shell of the solute show a bulklike diffusion rate of DCa-525(H2O) ) 2.83. In comparison, the diffusion rate of these water molecules calculated for Ca-123, DCa-123(H2O) ) 1.99, is reduced by one-third. This reduction might in part be explained by the reduced average simulation temperature TCa-123 ) 290 K as compared to TCa-525 ) 300 K (cf. Figure 1 and section 3.1). On the other hand, an influence of the difference in concentrations between Ca-123 and Ca-525 cannot be ruled out on the basis of these simulations. For Ca2+, a strong dependency of the value of D from the box size was observed: for the small system, DCa-123(Ca2+) ) 0.29 reached only 55% of the large system’s value DCa-525(Ca2+) ) 0.55. These differences cannot be explained solely by different simulation temperaturessthis would account only for 50% of the effect. Obviously, a small system size gives rise to results which are not reliable, as can be judged by comparing these values with the experimental results reported by Hertz3 and Weinga¨rtner.11 It must be noted though that, even with the large ionic systems containing 525 water molecules, our simulation results for the divalent ions slightly underestimate the experimental value of D(ion). The strong system size dependence of the values of D for the divalent Ca2+ ion can be interpreted as arising from a strong influence of this divalent ion on water molecules in the first and second hydration shell. Even outside the first hydration shell, the water molecules in Ca-123 do not reach bulk properties, as mentioned above. The complex consisting of the central ion and its tight hydration shell, the so-called solventberg, almost completely fills up the simulation box. Due to the limited size of the small system, this solventberg is not surrounded by a sufficiently large number of bulk water molecules which could enable a free diffusive motion of the ion-solvent complex.

15684 J. Phys. Chem., Vol. 100, No. 39, 1996

Obst and Bradaczek

TABLE 5: Diffusion Coefficient D of Water and Ions in Simulation Boxes of Different Sizes and Comparison to Literaturea Mr (g/mol)

D(Sim-525) ( SD (10-5 cm2/s)

D(Sim-125) ( SD (10-5 cm2/s)

D(lit.) (10-5 cm2/s)

D(exptl) (10-5 cm2/s)

r(cryst) (Å)

r(calcd) (Å)

6.94 22.99 39.10 24.31 40.08 87.62 18.02

1.06 ( 0.06 1.22 ( 0.05 1.80 ( 0.08 0.62 ( 0.09 0.55 ( 0.04 0.54 ( 0.07 2.87 ( 0.05

nd nd nd nd 0.29 ( 0.06 nd 2.74 ( 0.11

0.66 ( 0.18(II) 0.78 ( 0.13(II) 0.84 ( 0.28(II)

1.03 1.33 1.96 0.70 0.79 0.77 2.30

0.73 1.16 1.65 0.86 1.26 1.32 1.44

2.54 2.11 1.29 4.11 4.48 4.65 0.86

Li+ Na+ K+ Mg2+ Ca2+ Sr2+ H2O

2.67(I)

a Mr: molecular mass.56 Sim-525: V ) (25.08 Å)3, 525 water molecules and one ion (Li+, Na+, K+, Mg2+, Ca2+, Sr2+) or 528 water molecules (H2O). Sim-125: V ) (15.50 Å)3, 123 water molecules and one Ca2+ or 125 water molecules (H2O). SD ) standard deviation of the mean value of D from the analysis of 1000 ps simulation time for the ionic solutions and of 100 ps for pure water. D(exptl): Hertz,3 Weinga¨rtner.11 D(lit.): (I) Zhang et al.;57 (II) Lee and Rasaiah.5 r(cryst): crystallographic effective ionic radii.47,48 r(cryst)H2O: Agmon.52 r(calcd): hydrodynamic radii rsol, calculated from D using Stokes’ law rsol ) kBT/(6πηD), with T ) temperature (K), kB ) Boltzmann’s constant, and η ) viscosity (ηwater(298 K) ) 0.891 gm-1 s-1).

TABLE 6: Mean Residence Time τ and Number nlongtime of Water Molecules Not Leaving the Hydration Sphere of Ions and Water within the Simulation Timea molecule t (ps) τImpey (ps) τLee/Rasaiah (ps) nlongtime (1000 ps) nShell tav (ps) tmax (ps)

Li+ 41.4 33.3 184 0 175 25.4 199.1

Na+ 14.7 9.9 38 0 351 17.6 137.3

K+ 8.3 4.8 15 0 441 17.3 145.7

Mg2+

Ca2+-123

Sr2+

H2O-125

421.9

nd

(I)

Ca2+-525 nd(I)

50.7

0 164 38.0 936.3

4 44 120.8 1000

6 47 169.7 1000

0 217 41.8 729.9

4.7 4.5 5.3 0 89 4.5 49.6

H2O-528 4.6 0 108 5.5 55.3

a (I): no value due to persistence of the hydration shell throughout the entire simulation phase. τ Impey: mean residence times in 64 (125) MCY water molecules at T ≈ 280 K, taken from Impey et al.4 τLee/Rasaiah: mean residence times in 216 TIP4P water molecules at T ) 298 K, taken from Lee and Rasaiah.5 nlongtime(1000 ps): number of molecules not leaving the hydration sphere within 1 ns simulation time (for the initial hydration number, the reader is referred to nhydr, given in Table 4). nShell: total number of water molecules which were found at least once inside the first hydration shell during the 1 ns simulation time. tav: average time water molecules spent in the hydration shell of the ion. tmax: maximum hydration time observed during the simulation; values of 1000 indicate an incomplete decay of the hydration shell during the 1 ns simulation phase.

There is quite a large difference between our results for D(monovalent ion) and those reported by Lee and Rasaiah,5 but our values are in good agreement with the experimental results compiled by Hertz.3 Nevertheless, all results agree in the finding that D(monovalent ion) increases with increasing molar mass, which accords with the simple picture of Stokes’ law predicting a reciprocal proportionality to the radius only if one takes not ionic radii but larger hydrodynamic radii.38 The calculated radii r(calcd) given in Table 5 do not correspond to the radii of the ionic hydration shells which do increase with increasing molar mass. For K+ and for H2O, the calculated radii are smaller than the crystallographic radii. Assuming the validity of Stokes’ law, this implies that the observed D values are too high. On the other hand, the radii calculated for the divalent ions and for Li+ and Na+ are larger than the crystallographic radii. This is interpreted as the effect of tightly bound solvent molecules of the first hydration shell which increases the effective hydrodynamic radii of these ions. However, no direct correlation between the size of the first hydration shell derived from the rdf and the effective hydrodynamic radii was found. Obviously, the dynamics of the hydration shell have to be taken into account. 3.5. Mean Residence Time. Residence times were evaluated in terms of maximum and average residence times of water molecules in the first hydration shell of the ions or a central water molecule. Mean residence times τ for water molecules in the first hydration shell of different ions were computed by following the description given by Impey et al.4 The mean residence time reflects the decay of the hydration shell of an ion in the course of time. A tolerance time t* ) 2 ps was used in accordance with the time chosen by Impey et al. Any molecule that does not leave the hydration shell for a time longer than t* without returning into it in the meantime is treated as not having left the hydration shell at all. This definition was

introduced by Impey et al. and has been used since that time with different values of t* by other authors.27,31,33 A large number of very short hydration times do result from solvent molecules located close to the border between the first and second hydration shells, which are oscillating between the inside and outside of the hydration shell. Any value derived from such a procedure reflects not only the strength of the hydration but also the choice of the time interval of coordinate saving and of the chosen radius of the hydration shell. The introduction of t* reduces the problems arising from the rattling motions of water molecules of the first hydration shell. In Table 6, mean residence times τ and average and maximum residence times are given. Mean residence times τ could not be evaluated for both systems containing Ca2+. The hydration shells of the calcium cations did not decay completely within the simulation time of 1000 ps. Even after more than 1 ns, six (four) water molecules, which were part of the initial hydration shell at t ) t0, still were present for Ca2+ in the large (small) system. Moreover, the high stability of the hydration shell of Ca2+ results in a significantly reduced number nShell of water molecules visiting the hydration shell within the 1 ns simulation time (cf. Table 6). The minimum values, nShell(Ca-123) ) 44 and nShell(Ca-123) ) 47, amount to only 10% of the maximum value, nShell(K+) ) 441. For Li+, τ is longer than τ(H2O) by 1 order of magnitude, whereas τ(K+) and τ(Na+) are only 2-3 times longer than τ(H2O). Our result for lithium is in good agreement with τ(Li+) ≈ 10-11 s estimated by Newsome et al.41 For both systems containing Ca2+, a very slow decay of the hydration shell was observed, but the observation of the complete decay is without the reach of our current computational resources, as one might expect that it would take more than another 2 ns to occur. Burgess38 gives mean residence times for water molecules in the primary hydration shell of τ(Ca2+) ≈ 1 ns and τ(Mg2+) ≈ 1

MD Simulation of Hydrated Cations

J. Phys. Chem., Vol. 100, No. 39, 1996 15685

Figure 10. Lennard-Jones interaction energies between water-oxygen and the alkaline-earth metal cations Mg2+ (dash dot), Ca2+ (dot). and Sr2+ (dash).

µs, the latter being definitely out of reach by means of conventional MD simulation at the present time for us. For pure water, Bellissent-Funel and Teixera42 derived a residence time of 1.25 ps at 293K from neutron scattering experiments, which is close to our result from the MD simulation of 4.5 ps. The average hydration time is not significantly affected by the system size. Whereas a decrease in average hydration times is observed with increasing molar mass of the monovalent ions, the average hydration time of the divalent ions is highest for the medium-sized Ca2+ and much lower for the smaller Mg2+ and for the bigger Sr2+ ion. This is in contradiction to experimental findings43 which show the following decreasing order in residence times: Sr2+ < Ca2+ , Mg2+. For water, the average hydration time is lower than for any of the ions simulated. For the monovalent ions, the maximum hydration time decreases in the same way the average hydration time does. For Ca2+, the maximum residence time could not be determined exactly, as it exceeds the simulation time. This is also indicated by the nonzero number nlongtime of water molecules not having left the hydration shell of Ca2+ at all. The residence times for the alkaline cations found by Impey et al.4 at T ≈ 280 K and by Lee and Rasaiah5 at T ) 298 K and the value of τ(Na+) ) 20.0 ps reported by Garcia and Stiller27 for TIPS3 water are in good agreement with our results at T ) 300 K (cf. Table 6). From Raman scattering experiments of aqueous solutions of LiCl, NaCl, KCl, MgCl2, and CaCl2, Wang and Tominaga44 concluded that the influence of monovalent ions on the dynamic structure of water is weaker than the effects of the divalent ions. They found the following order of water relaxation times: τ(Mg2+) > τ(Ca2+) > τ(Li+) > τ(Na+) > τ(K+). This is in accord with the data given by Burgess,38 who gives residence times for Mg2+ which are by several orders of magnitude larger than those of Ca2+. On the contrary, in our simulations, less exchange of water molecules from the first hydration shell of Ca2+ than of Mg2+ occurred. This might be due to inconsistencies of the CHARMM22 force field used. According to Degre`ve and Quintal,12 solvent molecules inside the hydration shell of a cation are energetically in a more stable state than bulk water molecules. They have to pass an energetic barrier to reach the bulk, whereas the kinetic energy of bulk water molecules usually is sufficient for passing the barrier and reaching the hydration shell. Therefore, the parameters of the ion-water interaction in the CHARMM22 force field might give rise to the observed effects for Ca.2+ and Mg2+. From comparison of the profiles of the Lennard-Jones interaction energies between a water-oxygen atom and the cations Mg2+, Ca2+, and Sr2+ (cf. Figure 10), the increased stability of the hydration shell of Ca2+ cannot be explained. Although the energy minimum for the Sr2+-water interaction is deeper than

Figure 11. Snapshots of the hydration shell structure of the alkalineearth metal cations Mg2+ (top), Ca2+ (middle), and Sr2+ (bottom). The stereo drawings were generated using SCHAKAL (E. Keller, Freiburg).

for the Ca2+-water interaction, the hydration shell of Ca2+ was more persistent. A high degree of order was found by Probst et al.15 for the water molecules in the first hydration shell of Ca2+ (distorted hexahedron) and Mg2+ (octahedron), which might give rise to the extraordinary stability of these hydration shells. In Figure 11, the alkaline-earth metal cations are shown together with their first hydration shells. The solvation shell of Mg2+ is arranged as an octahedron, whereas the water molecules surrounding Ca2+ form a square antiprism. Sr2+ is hydrated by nine water molecules and their arrangement can be described as a trigonal prism with an extra ligand above the face of each rectangular prism. 3.6. Jump Length and Angle. In a very simple approach, the macroscopic diffusion coefficient D and the average jump angle and length should be related. To gain deeper insight into the molecular mechanism leading to the discrepancy between the simple Stokes’ law and the observed diffusion coefficients, we decided to take a closer look at the elementary steps of the simulation. Therefore, a new simulation was started for each ion using coordinates and velocities from a restart file written after 500 ps of simulation. It consisted of a 3 ps simulation time with an interval of saving energies, coordinates, and velocities which was equal to the simulation time step of 1 fs. This allowed us to investigate the distance traveled and the change of direction between two subsequent simulation time steps. For the water molecules, data refer to the center-of-mass motion.

15686 J. Phys. Chem., Vol. 100, No. 39, 1996

Obst and Bradaczek

TABLE 7: Average Angle θ between Subsequent Velocity Vectors Calculated at tn and tn+1a molecule

Li+

Na+

K+

Mg2+

Ca2+-123

Ca2+-525

Sr2+

H2O-125

H2O-528

D (10-5 cm2/s) θ(ion) (deg) θ(H2O-in) (deg) θ(H2O-out) (deg)

1.06 4.28 1.67 1.34

1.22 1.77 1.40 1.33

1.80 1.01 1.31 1.33

0.62 3.92 2.13 1.34

0.29 2.17 1.74 1.34

0.55 2.05 1.66 1.32

0.54 1.53 1.45 1.33

2.74 1.32 1.37 1.36

2.87 1.09 1.36 1.35

a The width of a time step is 1 fs. D: diffusion coefficient determined from simulation data through the Einstein relation (cf. section 3.4). θ(ion): angle θ, determined for the ion denoted under molecule. θ(H2O-in): θ for all water molecules H2O-in which are located within the first hydration shell (cf. Table 3) of the specific molecule at a given time step tn+1 without taking into account its location at tn. θ(H2O-out): θ for all water molecules H2O-out, which are located outside the first hydration shell, i.e. which do not belong to H2O-in.

TABLE 8: Average Distance Traveled between Subsequent Simulation Time Steps tn, tn+1 of 1 fsa molecule

Li+

Na+

K+

Mg2+

Ca2+-123

Ca2+-525

Sr2+

H2O-125

H2O-528

distance(ion) (Å) distance(H2O-in) (Å) distance(H2O-out) (Å) D (10-5 cm2/s)

0.43 0.28 0.28 1.06

0.21 0.27 0.28 1.22

0.19 0.29 0.28 1.80

0.22 0.27 0.28 0.62

0.20 0.28 0.28 0.29

0.17 0.28 0.28 0.55

0.12 0.28 0.28 0.54

0.31 0.28 0.28 2.74

0.26 0.27 0.28 2.87

adistance(ion): average distance traveled by the ion during 1 fs. distance(H O-in): average distance traveled during 1 fs by water molecules 2 H2O-in belonging to the first hydration shell (cf. Table 7). distance(H2O-out): average distance traveled during 1 fs by water molecules H2O-out not belonging to the first hydration shell. D: diffusion coefficient.

Figure 12. Relation between the average jump angle θ and the molecular mass Mr. Mr56 is given in g/mol; θ is given in degrees. The lines connecting the data points are for visibility only and have no physical meaning.

In Table 7, the average angle θ between the velocity vectors computed at two subsequent simulation steps of 1 fs is compared to the diffusion coefficient D. No direct relation between the jump angle and diffusion coefficient can be seen. Instead, for both mono- and divalent ions, θ decreases with increasing atomic mass (cf. Figure 12). This effect can be interpreted in terms of classical mechanics as arising from the encounter between two particles of different mass. Therefore, a more complex relationship between the elementary events of subsequent jumps and macroscopic diffusion must exist, probably involving effects arising from the presence of water molecules in the first hydration shell. A small difference between the θ values determined for the large (Ca-525) and the small system (Ca-123) containing Ca2+ was noticed both for the ion and for water within the first hydration shell. On the other hand, no such difference was found between both systems of different size comprising of pure water only. For the water molecules outside the first hydration shell, no significant difference can be seen between the water molecules in different ionic solutions and pure water. Water molecules inside the first hydration shell show a non-bulklike behavior for all ions except K+. The effect is more pronounced for the divalent ions than for the monovalent ions. This finding is in accord with the order of relaxation times found by Wang and Tominaga,44 τ(Mg2+) > τ(Ca2+) > τ(Li+) > τ(Na+) > τ(K+) (cf. section 3.5). The average jump length, i.e., the distance traveled during the 1 fs time step, decreases with increasing molar mass of the

ion (cf. Table 8). The average jump length of water molecules is the same both inside and outside the first hydration shell. For Ca2+, the jump distance was smaller for the large system than for the small system simulated, whereas for water, no size dependence of the average distance was detected. Nevertheless, the possibility cannot be ruled out that these differences do not arise from effects of the system size but from statistical inaccuracy. The uncertainty of the results for the ions may be estimated by comparing the values for bulk water H2O-out and H2O-in (averaged over the complete ensembles of 528 and 125 molecules, respectively) with those obtained for one randomly chosen molecule of water (instead of an ion) in the pure water simulations. Whereas the error for the distances is small, the large difference between the θ values for bulk water and for an individual water molecule is dissatisfying. Nevertheless, the differences in θ of the different ions are significantly larger than the estimated error. The average jump length of water molecules of the first hydration shell is not affected by the type of ion in the center, whereas the average jump angle decreases with increasing molar mass of the ion (cf. Figure 12) and also depends on the charge state of the ion. Divalent ions lead to larger deflection angles than monovalent ions. In these terms, a water molecule in the center gives rise to results comparable to those of Na+ or K+. 4. Summary We have investigated static and dynamic properties of solutions of alkaline and alkaline-earth metal cations in water by means of molecular dynamics simulations. The radial distribution function of water was not influenced by the size of the simulated system. The local density of water in the vicinity of the cations in terms of rdf peak heights depends on the charge and the radius of the ion. Divalent cations and the small Li+ ion lead to a higher local density than the larger monovalent ions Na+ and K+. In general, good agreement was found between experimental results and our results, although in the simulations, the solvent shell structure extends farther into the water than in the experiments. The coordination number increases with increasing ionic radius. Compared to literature data, ncoor(Ca2+) seems slightly underestimated. The hydrogen bond distribution reaches its equilibrium quickly, as can be judged from comparison between values calculated from different time segments. For pure water and Ca2+, n(HB) varies with the system size, reaching values which are approximately 10% higher in the large systems than in the

MD Simulation of Hydrated Cations small ones, which seem to give unreasonably low results. On the average, water in the large ionic solution systems adopts bulk water properties. The structure of the hydration shell as indicated by the rdf influences the spatial distribution of hydrogen bonds, nHB(r). Water molecules can form hydrogen bonds to other water molecules, whereas the cations are not capable of doing so, thus leading to a distortion of the hydrogen bond pattern which gives rise to the observed differences in nHB(r). Diffusion coefficients D from our simulations were in good agreement with experimental results. D decreases with decreasing system size for both pure water and Ca2+. For the monovalent ions, D increases with increasing molar mass, whereas for the divalent ions, D slightly decreases with increasing molar mass. The values of D are not directly related to neither the average jump length nor the average jump angle θ between subsequent elementary (1 fs) steps of the simulation. The following order of mean residence times τ was found: τ(Mg2+) . τ(Sr2+) > τ(Li+) > τ(Na+) > τ(K+) > τ(H2O). Mean residence times τ were not determined for Ca2+ because the hydration shell of this cation did not decay completely within the simulation time of 1 ns. This was probably caused by the highly ordered structure of the Ca2+ solvent shell. The mean residence time for pure water, τ(H2O), increases with increasing system size. Nevertheless, in both systems, τ(H2O) is smaller than the residence times τ observed for the alkaline metal cations. The observed effects of the system size on the molecular dynamics simulation results strongly suggest the use of larger systems for the determination of the properties of ionic solutions. To investigate the dynamic properties of the alkaline-earth metal cations more thoroughly, simulation times should be increased by at least 1 order of magnitude. The good overall agreement of our simulation results with experimental data underlines the reliability of the CHARMM22 force field for simulations involving cations, i.e., simulations of charged biological macromolecules such as DNA or proteins with their counterions. Only the dynamic properties of the Ca2+ ion were not reproduced in a completely satisfactory manner. Acknowledgment. The authors thank Dr. T. Gutberlet, Dr. M. Kastowsky, P.-J. Koch, and Dr. A. Sabisch for valuable discussions. References and Notes (1) Einstein, A. Ann. Phys. 1905, 17, 549. (2) Atkins, P. W. Physical Chemistry, 3rd ed.; Oxford University Press: Oxford, U.K.,1988; Chapter 27. (3) Hertz, H. G. In Water: A ComprehensiVe Treatise; Franks, F., Ed.; Plenum Press: New York, 1973; Vol. 3, Chapter 7. (4) Impey, R. W.; Madden, P. A.; McDonald, I. R. J. Phys. Chem. 1983, 87, 5071. (5) Lee, S. H.; Rasaiah, J. C. J. Chem. Phys. 1994, 101, 6964. (6) Sansom, M. S. P.; Kerr, I. D.; Breed, J.; Sankararamakrishnan, R. Biophys. J., 1996, 70, 693. (7) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (8) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (9) Allen M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K.,1987; Chapter 1. (10) van Gunsteren, W. F.; Berendsen, H. J. C. Mol. Phys. 1977, 34, 1311.

J. Phys. Chem., Vol. 100, No. 39, 1996 15687 (11) Weinga¨rtner, H. Z. Phys. Chem. NF 1982, 132, 129. (12) Degre`ve, L.; Quintale, C., Jr. J. Chem. Phys. 1994, 101, 2319. (13) Andrea T. A.; Swope W. C.; Andersen, H. C. J. Chem. Phys. 1983, 79, 4567. (14) Brooks, C. L., III. J. Chem. Phys. 1987, 86, 5156. (15) Probst, M. M.; Radnai, T.; Heinzinger, K.; Bopp, P.; Rode, B. M. J. Phys. Chem., 1985, 89, 753. (16) Neilson, G. W.; Enderby, J. E. AdV. Inorg. Chem. 1989, 34, 195. (17) Neilson, G. W.; Broadbent, R. D. Chem. Phys. Lett. 1990, 167, 429. (18) Caminiti, R.; Musinu, A.; Paschina, G.; Pinna, G. J. Appl. Crystallogr. 1982, 15, 482. (19) Enderby, J. E. Chem. Soc. ReV. 1995, 23, 159. (20) Petitt, B. M.; Rossky, P. J. J. Chem. Phys. 1986, 84, 5836. (21) Neilson, G. W.; Enderby, J. E. In Water: A ComprehensiVe Treatise; Franks, F., Ed.; Plenum Press: New York, 1979; Vol. 6, Chapter 1. (22) Hewish, N. A.; Neilson, G. W.; Enderby, J. E. Nature 1982, 297, 138. (23) Skipper, N. T.; Neilson, G. W.; Cummings, S. C. J. Phys.: Condens. Matter, 1989, 1, 3489. (24) Soper, A. K.; Phillips, M. G. Chem. Phys. 1986, 107, 47. (25) Narten, A. H.; Levy, H. A. In Water: A ComprehensiVe Treatise; Franks, F., Ed.; Plenum Press: New York, 1972; Vol. 1, Chapter 1. (26) Gaballa, G. A.; Neilson, G. W. Mol. Phys. 1983, 50, 97. (27) Garcia, A. E.; Stiller, L. J. Comput. Chem. 1993, 14, 1396. (28) Sciortino, F.; Geiger, A.; Stanley, H. E. Nature 1991, 354, 218. (29) Heinzinger, K. Pure Appl. Chem. 1985, 57, 1031. (30) Balbuena, P. B.; Johnston, K. P.; Rossky, P. J. J. Phys. Chem. 1996, 100, 2706. (31) Kowall, T. Molekulardynamische Simulation der Hydratation des Kronenethers 18-Krone-6 und der Komplexierung eines Kalium-Kations. Thesis, Universita¨t Dortmund, Germany, 1992. (32) Kru¨ger, P.; Strassburger, W.; Wollmer, A.; van Gunsteren, W. F. Eur. Biophys. J. 1985, 13, 77. (33) Kowall, T.; Geiger, A. J. Phys. Chem. 1994, 98, 6216. (34) Falk, M.; Knop, O. In Water: A ComprehensiVe Treatise; Franks, F., Ed.; Plenum Press: New York, 1973; Vol. 2, Chapter 2. (35) Ohmine, I.; Tanaka, H. Chem. ReV. 1993, 93, 2545. (36) Sciortino, F.; Fornili, S. L. J. Chem. Phys. 1989, 90, 2786. (37) Geiger, A.; Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1979, 70, 4185. (38) Burgess, J. Ions in Solution: Basic Principles of Chemical Interactions; Ellis Horwood: Chichester, U.K., 1988. (39) Vaisman. I. I.; Brown, F. K.; Tropsha, A. J. Phys. Chem. 1994, 98, 5559. (40) Svishev, I. M.; Kusalik, P. G. J. Phys. Chem. 1994, 98, 728. (41) Newsome, J. R.; Neilson, G. W.; Enderby, J. E. J. Phys. C 1980, 13, L923. (42) Bellissent-Funel, M. C.; Teixera, J. J. Mol. Struct. 1991, 250, 213. (43) Frey, U.; Merbach, A. E.; Powell, D. H. In Dynamics of Solutions and Fluid Mixtures by NMR; Delpuech, J.-J., Ed.; John Wiley & Sons: Chichester, U.K., 1995; Chapter 6. (44) Wang, Y.; Tominaga, Y. J. Chem. Phys. 1994, 101, 3453. (45) Dietz, W.; Riede, W. O.; Heinzinger, K. Z. Naturforsch. 1982, 37a, 1038. (46) Powell, D. H.; Neilson, G. W.; Enderby, J. E. J. Phys.: Condens. Matter 1989, 1, 8721. (47) Shannon, R. D. Acta Crystallogr. 1976, A32, 751. (48) Shannon, R. D.; Prewitt, C. T. Acta Crystallogr. 1969, B25, 925. (49) Svishev, I. M.; Kusalik, P. G. J. Chem. Phys. 1993, 99, 3049. (50) Forester, T. R.; Smith, W.; Clarke, J. H. R. J. Phys. Chem. 1995, 99, 14418. (51) Gogonea, V.; Osawa, E. J. Comput. Chem. 1995, 16, 817. (52) Agmon, N. J. Phys. Chem. 1996, 100, 1072. (53) Skipper, N. T.; Neilson, G. W. J. Phys.: Condens. Matter 1989, 1, 4141. (54) Neilson, G. W.; Skipper, N. T. Chem. Phys. Lett. 1985, 114, 35. (55) Hinton, J. F.; Amis, E. S. Chem. ReV. 1971, 71, 627. (56) IUPAC. Pure Appl. Chem. 1991, 64, 1520. (57) Zhang, L.; Ted Davis, H.; Kroll, D. M.; White, H. S. J. Phys. Chem. 1995, 99, 2878.

JP961384B