J . Phys. Chem. 1993,97, 4 5 0 0 4 5 0 4
4500
Recovering the Concept of the Hydrated Ion for Modeling Ionic Solutions: A Monte Carlo Study of Zn2+in Water Rafael R. Pappalardo and Enrique Shchez Marcos' Departamento de Quimica Fisica, Facultad de Quimica, Universidad de Sevilla, 41 01 2 Sevilla, Spain Received: September 30, 1992
A [Zn(H20)6]2+-H20 interaction potential has been developed by means of ab initio calculations using the 3-21G* basis sets. Monte Carlo simulations of zinc ion in water have been performed by using the potential for the zinc hexahydrate and for comparative purposes that of Clementi developed for the single ion. Both simulations have used the same conditions. Structural results derived from simulation of zinc hexahydrate give a more disrupted second shell, with two maxima for the Zn-H radial distribution function at the second shell. Results for the hydration energy of Zn2+obtained from the hydrated ion approach (-5 18 kcal/mol) are much closer to the experimental estimation (-491 kcal/mol) than results obtained with the potential for the single ion (-845 kcal/mol). As a matter of fact, the main conclusion of this study is that the M C and MD simulations of ions in solution using the concept of a species beyond the single ion is not unrealistic and can be, under some conditions, an interesting alternative to be considered in the case of conflicting simulations of highly charged or transition metal ions. Advantages and limitations of this approximation are discussed.
TABLE I: Energetic and Geometrical Data of the
Introduction
[Zn(H20),#+ Cluster and H20
Several methods for the construction of interaction potentials between the components of a solution have been proposed. The most general method is that which uses the interaction energies derived from quantum mechanical calculations, principly for its independence from known experimental data.' Such a method allows the obtention of energetic information in the complete potential hypersurface. The case of highly charged metallic cations, Mn+-H20, usually presents several problems due to the fact, first pointed out by Clementi et that the correct dissociation products are M('+-l)+ H20+ instead of Mn++ H20, owing to the relative magnitude of the ionization potentials for the metallic cation and water. More recently, Cossi and Persico) have examined in great detail this problem of charge transfer and curve crossing in the [Be(H20)l2+system. Transition metal cation-water potential surfaces present crossing between electronic states at cation-water distancesshort enough to precludethe determination of the interaction potential. For the [Fe(H20)]2+ the crossing occurs at distances greater than 6 A, whereas that for the [Fe(H20)]3+ appears at 1.85 A, much closer to the Fe3+-H20 m i n i m ~ m .On ~ the other hand, the normal way to determine the interaction potential is using the pairwise approximation that seems to be particularly crude for transition-metal cations.5 This has motivated the inclusion of three-body corrections in the intermolecular potential.G8 In this paper, we want to investigate a different approach to the problem by using the hydrated cation as the entity to be considered in the determination of the ion-water interaction potential. The employment of the cation with its first hydration shell prevents the problems of charge transfer and state crossing, since this hydrated cation is stable enough against the interaction with a water molecule to avoid a charge transfer or a change of electronic state. In addition, the explicit consideration of several solvent molecules in the solute introduces implicitly some terms other than thoseof two-body interactions. Due to the preliminary character of the present study, which is devoted to the exploration of the above-mentioned approximation, the selected system has been the Zn2+cation, which experimentally is known to have a well-defined first hydration shell9 formed by six water molecules, and in a previous study it was shown that for the [Zn(H20)l2+ potential surfaces there is a charge transfer at short distances.1° Likewise, Rode et al." have extensively studied solutions of this
+
ET (ua) dwi (A) LHOH(deg) dzn-0
(A)
[Zn(Hz0)6I2+
H20
-2222.519 083 0.9690 108.5 2.067
-75.585 387 0.9572 104.5
cation by means of Monte Carlo simulationswhere intermolecular potentials including threebody effects have been considered. Construction of the [Zn(H20)#+-H20 IntermolecularPotential. The structure of the cluster [Zn(H20),l2+was determined by optimization of the geometry inside a spherical cavity using a SCRF method described elsewhere12to take into account average solvent effects. The dielectric permittivity of water was chosen for the dielectric continuum and a multipolar development up to 6th order was used. The basis set was 3-21G for H and 0 and 3-21G* for Zn in view of the large number of calculationsneeded to determine the intermolecular potential. The final geometric and energetic parameters for the cluster and the water molecule are in Table I. The reference system to determine the ion-solvent interaction potential will be zinc hexahydrate at the geometry optimized inside the cavity and a water molecule at the experimental geometry, that is, the geometry used for the MCY potential13 employed for water-water interactions. Due to the high symmetry of cluster, Th, the space to be scanned when constructing the interaction hypersurface can be reduced to a quadrant. In addition, certain privileged directions on the Zn-Ol axes and OIZnOz planes are equal by symmetry. Thus, the following parameters have been used (Scheme I): the distance Zn-07 between 2.0 and 9.0 A with variable increments to a total of 10 distances; the Zn0107 angle which was changed between 0' and 90' with increments of 30'; the dihedral angle between the 01ZnO2 and OIZn07 planes between Oo and 90° with intervals of 30'; the angle between the dipole moment of the seventh water molecule and the axis ZnOl between 0' and 180' with intervals of 45'; and the dihedral angle between the plane of the seventh water moleculeand theplaneOIZn07 which was changed between 0' and 180' with intervals of 45'. In principle, this sum implies 4000 points. Using the symmetry and fixing some parameters for distances Zn-07 greater than 6 A, the number of points is reduced to 1600. This involves about 700 CPU hours of a CONVEX 240.
0022-3654/93/2091-4500S04.00/0 0 1993 American Chemical Society
Monte Carlo Study of Zn2+in Water
The Journal of Physical Chemistry, Vol. 97, No. 17, 1993 4501
SCHEME I: Only the Oxygens of the Hexahydrated Ion Represented, and Thus a Oh Symmetry Shown. However, When the Hydrogens of the Water Molecules Are Considered the Symmetry Reduces to T,,,
CHART I
z
t
0 0 -Zn /01 X-
I’
I /O-0,
-Y
/ I
H 0
7(
H
The fitting of these values to an analytical function was performed by means of a least-squares procedure. As the most time-consuming part of a Monte Carlo simulation is the computation of the internal energy, the use of exponential functions was avoided and only terms in rnwere considered, where r is the distance between an atom of the cluster and an atom of the water molecule. There is a total of 57 of such distances. The explicit form of the potential employed after some trials is given in Chart I. The terms with n < 3 have been grouped in such a way that the integral J4?rr2E(r)dr vanishes as r m to avoid size dependence.I4 The values for the adjustable coefficients, C,, are collected in Table 11. With these Ci)s, the interaction energy is obtained in hartrees when distances are expressed in bohrs. For the rl terms, atomic charges have not been used as coefficients since they can be considered formally included in the fitted parameters Ci. The standard deviation of the complete set is 8.8 kcal/mol due to the high dispersion of the repulsive points. However, when only the points with energieslower than 15 kcal/ mol are considered, the standard deviation reduces to 1.8 kcal/ mol. Importance of the highly repulsive points during the simulationsis minimal, but these points are necessary during the fitting in order to avoid holes in the fitted potential, particularly in the central region of the quadrants. Figure 1 shows a plot of fitted vs computed energies. Another way to test the validity of the fitted functions is through the trace of isoenergetic contour maps. In Figure 2 is shown a set of such maps for a quadrant, demonstrating that there are no anomalous features. Monte Carlo Simulationsof Zn2+in Water. The Monte Carlo simulation has been performed in the canonical ensemble (constant NVT) using the algorithm of Metropolis et al.ls as implemented by Corongiu and Clementi,I6and modified by Lluch and Galera.” The solvent-solvent interaction energy was calculated by means of the MCY pair potential developed for the H20-H20 system.13 Analysis of the results was done using programs developed by Lluch and Galera.” We have performed two different simulations. In the first the potential developed by Clementi et was employed for the Zn2+-H20 system with 216 water molecules in a cubic box with adensityof 1 gcm-3. Fortheother simulation the [Zn(HzO)#+H20 intermolecular potential previously described was used, the basic box contained the zinc hexahydrateand 210 water molecules. Periodic boundary conditions were employed. Equilibration of the system at 298 K needed 3000K configurations for the simulation using thesingle ion potential, and 8000K configurations for that of the hydrated ion. Afterwards, 500K and lOOOK configurations were employed for analysis, respectively. The radial distribution functions (RDF) for Zn-0 and Zn-H and their integrations for the simulation where the single zinc ion potential proposed by Clementi is used are quite similar to those obtained by Clementi.13J8The gz,-(r) function shows a sharp peak at 2.0 A corresponding to the first hydration shell, whose integration leadsto a number of six water molecules. ThegZ,_H(r) function also shows a sharp peak centered at 2.7 A, whose integration gives 12 hydrogen atoms. The distance between the
-
/
2
12
2
2
6
2
6
main peaks of gzn-o(r) and gZn-H(r) indicates that the water molecules of this shell present a clear dipole orientation. The fact that both RDFs vanish once the first shell is formed means that this shell has quite a rigid structure. From the Zn-0 and Zn-H RDF a second hydration shell is also deduced, although, as expected, the peaks are broader than those of the first shell. The maximum for Zn-0 is at 4.0 A, and the integration value reaches 22 at the minimum after the second shell, corresponding to an average of 16 second shell water molecules. The maximum for Zn-H RDF is centered at 4.6 A. Rode et al.llahave recently critically reexamined this Clementi potential and have concluded that the reasonable results obtained are rather the consequence of some error compensation based on
Pappalardo and Sanchez Marcos
4502 The Journal of Physical Chemistry, Vol. 97, No. 17, 1993 300
I
I
I
I
a 200
,
100
p"'.. . ..
Z=O
..**t.
2.4 ..:.
I
I
I
I
i
0
100
200
300
0
Y
,E ,,
h
I
J
/ k c a l mol-'
Figure 1. Plot of fitted vs SCF-calculated energy points of the [Zn(H20)6I2+/H20 energy surface.
TABLE II: Adjusted Values of the Coefficients, C, for the Analytical Intermolecular Potential of the [Zn(H20),#+-H20 System coefficient
adiusted value
-4.541 I47 867 12.231 961 07 -14.376 274 55 8.211 474 469 -19.940 122 32 -46.878 472 61 -6.563 574 938 13.082 100 43 202.375 147 1 7.589 379 516 -19.348 410 1 1
coefficient CI 2 CI 3 cl4
CI 5 CI 6 CI 7 Cl8 C19 c 2 0 c 2I
adiusted value 23.530484 22 76 473.515 98 -0.762 834 951 0.076 097 195 23.649 773 83 -516.343 457 1 -1572.157 284 230 952.2185 628.148 879 8 8775.933471
the too simple pair-potential function and the minimal basis sets employed in the quantum chemical calculations to build this potential. These authors have developed a more involved interaction potential, including three-body corrections and calculating the interaction energy from quantum chemical calculations with DZV basis sets and pseudopotentials. They obtained RDF for Zn-O and Zn-H where two hydration shells are present, the first Zn-0 peak at 2.05 A, accompanied by a smaller spike at 2.25 A whose integration leads to seven molecules, and the first Zn-H peak at 3.0 A which integrates at 14-15 atoms. A second hydration shell clearer than that obtained from the Clementi potential appears, with maxima located at 4.20and 5.0 A for Zn-O and Zn-H, respectively. Integration leads to an average of 16-18 second-shell water molecules. Figure 3 plots the Zn-0 and Zn-H RDF and their integration resulting from the simulation, which uses the potential developed for the zinc hydrate. Obviously, since the first hydration shell is implicit in the solute, our RDFs deal with the second shell. The Zn-0 RDF shows a well-defined peak centered at 3.8-3.9 A, which integrates to 12water molecules. Contrary to the previous results, the Zn-H RDF presents two peaks at lower distances (3.7A) and greater (4.7A) than the maximum for Zn-0, whose integrations lead to 7.3 and 20.4 hydrogen atoms, respectively. This result seems to indicate that the picture for the second hydration shell corresponds to a structure where the ordered hydrogen bonding with water molecules of the first hydration shell is partially lost. Moreover, since the RDF does not go to zero at the minima, mixing between both types of hydrogens as well as with hydrogensof outer shells is allowed. This is obviously a consequenceof our potential and can be understood by realizing that the cluster-water interaction includes a detailed description of how water molecules can interact with the doubly-charged
z=2A
Figure 2. Isoenergetic contour maps built from the fitted function for the [Zn(Hz0)6]2+-H20 interaction. At each point the energy value for the most favorable orientation of the water molecule was chosen. (I = 0 is the OIZn02plane in Scheme I).
cation compromised simultaneously with the water molecules directly attached to this cation. In fact, as a simplified vision, our interaction potential is less directional and less attractive than those raising from consideration of the single ion. OrtegaBlake et al.19 have suggested similar ideas about the loss of hydrogen bond directionality in the Ca2+and Mg2+ hydration when nonadditivity was considered. Therefore, it is to be expected that water molecules in the second shell are allowed to get a more compromised orientation with water-water interactions. Comparing the twoRDFs it seems that thiseffect influencesthe relative orientation of water molecules rather than the average distance of these molecules to the central ion. Table I11 presents the water-solute, the water-water, and the total interaction energy per water moleculefor the two simulations performed. As is expected, the solutesolvent interaction energy is much smaller in the simulation of the zinc hydrate than in the simulation of the single zinc cation. Likewise,the solvent-solvent interaction energy is greater for the simulation of the hydrated
Monte Carlo Study of Zn2+ in Water
The Journal of Physical Chemistry, Vol. 97, No. 17, 1993 4503 26
1.5
Thus, the hydration energy of Zn2+ becomes
24
Mhrdr(ZnZ+) = -239.4 - 331.7 + 6(8.92) kcal/mol =
22 20
-5 17.6 kcal/mol
18
1.0
16
14
v I
12
p"
-
5 zo
10
0.5
8 6
4
0.0
'
2
'
'0 8
6
4
10
r (A) 3.5
2.5
.5
3'0 2.0
'
70
I
/
I
1
60 50
40 v
'z
u
30 1.0
1j
Concluding Remarks
i
/'
20
Figure 3. Radial distribution functions (solid lines) and their integration (dashed lines) from Monte Carlo testruns for [ Z I I ( H ~ O ) ~ ] ~in+ water.
TABLE III: Cation-Water, Water-Water, and Total Interaction Energy (kcal/mol) per Water Molecule Obtained from Simulations for Single Ion (216 HzO) and Hydrated Ion (210 H20) ~~
~
~~~
~
~~
~~~~
system
E(Zn2+-H20)
E(H20-H10)
Zn2+ 216 H,O Zn(H20)62+ 4 2 1 0 H20
-5.32 -1.83
-7.51 -8.23
+
This result is quite satisfactory since the experimental estimation is -491 kcal/mol.22 The dramatic difference to previous simulations should be again ascribed to the interaction potential used. When the interaction energy is computed for the single cationwater system, the strong interaction energy at distances of the first hydration shell are the same for all water molecules surrounding the cation. During the simulation this becomes an artifact because, for instance, for a first hydration shell octahedrally coordinated such as Zn2+, six times the same strong interaction energy is accounted, when it is well-known that the total interaction energy per water molecule becomes less and less attractive as the cluster size is enlarged.23Thus, when simulations deal with high charge solute, theoverestimation of the interaction energy will be particularly high, as in fact has been the case. On thecontrary, if thestrong interaction of the first water molecules is considered by quantum-chemical techniques as a whole, the previous deficiences are eliminated. Although, probably the nonadditivity of interaction remains on the second shell, its influence on the final results seems to be much smaller than that on the first shell.
-12.83 -10.06
ion than for that of the single ion. Taking into account that the evolution of this energy in simulation of ions with an increasing number of water molecules is toward more and more negative values,20pure water simulation being the limit (with the MCY potential,-8.92 kcal/mol),onecandeduce that the water structure in the simulation of the hydrated ions is more similar to that of the pure water than the water structure in the simulation with a single ion. Again, this is not surprising as in our simulation the water molecules most affected by the ion do not contribute to this value. For simulation of the single ion, the hydration energy of Zn2+ can be calculated as the difference between the total interaction energy of the [Zn(H20)2,6l2+system and the interaction energy of 2 16 water molecules. The value obtained is -844.6 kcal/mol. This result is similar to that obtained by Rode et al,I1awhich is -848 kcal/mol, using the potential with three-body corrections. For simulation of the zinc hexahydrate the hydration energy is -239.4 kcal/mol. To calculate the corresponding value for Zn2+ hydration we have to add two additional terms: the first is the energy to transfer six water molecules of the pure water solution to gas phase to form the cluster, which in our model will be, with a change of sign, six times the water-water interaction energy in the pure water simulation; the second is the formation energy of the cluster at 298 K, which has already been calculated in a previous paper21 for the same basis set, Le., -331.7 kcal/mol.
In this work we have examined the possibility of developing simulations of ionic solutions based on the concept of hydrated ion. Although construction of the hydrated ion-water potential is more complex than for single ion-water, when cases of high charge or transition-metal cations are studied, this technique represents an alternative to the many-body corrections. In the case of open-shell metallic cations, the stabilization of a given electronic state by coordination to solvent molecules can easily be reached, opening the possibility of extracting an interaction potential with another solvent molecule from a well-defined potential surface. This is not the case if the potential is developed from a single ion-solvent molecule system. It has been found that the structure of the second hydration shell does not follow a hydrogen bond directionality as strong as that found in the simulationsemployingthe singleion-water potentials. This result agree with the conclusions of Torrie and PateyZ4 about water structure in the proximity of the electrode, where beyond the first layer the solvent-solvent intermolecularforces are quite important in determining the interfacial structure. Likewise, it has been found that the theoretical hydration energy obtained as a sum of statistical and quantum-mechanical contributions agrees fairly well with the experimental value, contrarily to simulations which use the single ion as solute. Obviously,the main limitations rise from thedefinition a priori of the first hydration shell. Thus, dynamics of the exchange of solvent molecules of this shell cannot be obtained with our potential. The unimolecular rate constant for water release from the first hydration shell of Zn2+ is 3.2 X lo7S - I , ~ ~ which represents a relatively short half-life time(ca. 20 ns). To some extent the ridigity of the molecules of the first hydration shell may also exaggerate the structure of the second shell; however, this limitation could be overcome by means of a flexible potential, where the first shell water molecules are allowed to move, basically following the rotation of the water plane around the Zn-O axis and, in addition, considering libration of these water molecules. Apart from these limitations,we believe that for thermodynamical properties of ionic solutions this type of approximation should be quite useful. For instance, properties depending on the ionic mobility could be better reproduced under this approximation. In fact, the hydrated ions are generally more clearly defined
4504 The Journal of Physical Chemistry, Vol. 97,No. 17,I993
particularly for the type of ions which present more difficulties in getting a correct intermolecularpotential and use to have large half-life time of water molecule residence in the first hydration shell. In any case more effort should be made in this direction, examining the solvation of more critical cations, such as Cr3+, Fe3+, Rh3+, etc.
Acknowledgment. We are grateful to Prof. Lluch and Drs. Gondlez-Lafont and Galera, University Autonoma de Barcelona, for supplying us with a copy of their program and Drs.Andreu and Molero, University of Sevilla, for useful discussions. We thank the Direcci6n General Cientlfica y Tknica of Spain (PB890642) and the Junta de Andalucfa (Group 1039) for financial support. References and Notes (1) Clementi, E. In Modern TEchniques in Computational Chemistry (MOTECO; Clementi, E.; Ed.; Escom: Leiden, 1990; Chapter 1. (2) Corongiu. G.; Clementi, E. J . Chem. Phys. 1978, 69, 4885. (3) Cossi, M.; Persico, M.Theor. Chim. Acta 1991, 81, 157. (4) Curtiss, L. A.; Halley, J. W.; Hautman, J.; Rahman, A. J . Chem. Phvs. 2319. . ..,. . 1987. -.- , 86. .., . .. ( 5 ) Bounds, P. G. Mol. Phys. 1985, 54, 1335. (6) Cordeiro, M.N. D. S.; Gomes, J. A. N. F.; GonAlez-Lafont, A.; Lluch, J. M; B e r t h , J. Chem. Phys. 1990, 141, 379. (7) Rode, B. M.; Islam, S. M.Z . Naturforsch., Teil A 1991, 46, 357.
Pappalardo and Sanchez Marcos (8) Curtiss, L. A.; Halley, J. W.; Hautman, J. Chem. Phys. 1989, 133, 89. Curtiss, L. A.; Jurgens, R. J . Phys. Chem. 1990, 94, 5509. (9) Magini, M.; Licheri, G.; Paschina, G.; Piccaluga, G.; Pinna, G. X-Ray of Ions in Aqueous Solutions: Hydration and Complex Formation; CRC Press: Boca Raton, 1988; pp 142-148. (10) Sanchez Marcos, E.; Pappalardo, R. R.; Barthelat, J. C.; Gadea, F. X . J. Phys. Chem. 1992, 96, 516. (1 1) (a) Yongyai, Y. P.; Kokpol, S.;Rode, B. M.Chem. Phys. 1991,156, 403. (b) Yongyai, Y.; Kokpol, S.;Rode, B. M.J. Chem. SOC.,Faraday Trans. 1992, 88, 15f7. (12) Rinaldi, D.; Rivail, J. L.; Rguini, N. J . Comput. Chem. 1992, 191, 287. (13) Matsuoka, 0.;Clementi, E.; Yoshimine, M. J . Chem. Phys. 1976, 64, 1351. (14) Curnoyer, M. E.; Jorgensen, W. L. J . Am. Chem. SOC.1978, 106, 4942. (15) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J . Chem. Phys. 1953, 21, 1087. (16) Romano, S.;Clementi, E. Gazz. Chim. Nal. 1978,108, 319. (17) (a) Lluch, J. M.; Galera. S. Privatecommunication. (b) Galera, S. Ph.D. Thesis, Universidad Aut6noma de Barcelona, 1991. (18) Clementi, E.; Corongiu, G.; Jdnsson, B.; Romano, S.J . Chem. Phys. 1980, 72, 260. (19) Ortega-Blake, I.; Hernandez, J.; Novaro, 0. J . Chem. Phys. 1984, 81, 1894. (20) Gonzalez-Lafont, A.; Lluch, J. M.; Oliva, A.; Bertran, J. Chem. Phys. 1987, 111, 241. (21) Sanchez Marcos, E.; Pappalardo, R. R.; Rinaldi, D. J . Phys. Chem. .. 1991, 95, 8928. (22) Marcus, Y. Ion Solvation; John Wiley: Chichester, 1985. (23) Pullman. A.: Demoulin. D. Int. J . Ouantum Chem. 1979. XVI. 641. (24) Torrie, G. hi.;Patey, G.N. E l e c k h i m i c a Acta 1991, 36, 1677. (25) Pearson, R. G.; Ellen, P. C. In Physical Chemistry. An Advanced Treatise; Eyring, H., Ed.; Academic Press: New York, 1975; Vol. VII, pp 228-230.
