- Email: [email protected]
The determinant role of water in the ionic dissociation of HO2
THEO CHEM ELSEVIER
Journal of Molecular
Structure (Theochem)
371 (1996) 143-152
The determinant role of water in the ionic dissociation of HOZ’ Carlos P&ez de1 Valle”, Carmela Valdemoro”,
Juan J. Novoah-*
~Institutode Matemcitica Aplicada y Fish Fundamental, CSIC Cl. Serrano 123, 28~-~udr~d, Spain hDept. Quimica Fisica, Fat. Quimica, Univ. de Barcelona Av. Diagonal 647, 08028-Barcelona. Spain Received
Abstract The energetics of the ionic dissociation of the hydroperoxyl radical (HOZ) into the H’ and 0; ions in solution is studied using ab initio methods and water clusters to mimic the solvent effects. The clusters selected are HO,..(H,O),, H’..(H@), and 0;..+(H20),, using various combinations of n, 1 and m designed to saturate the first solvation shell of the involved systems, the largest of which contains 10 water molecules (n = 10) plus the electrostatic field created by the outer solvation shell to these 10 water molecules. We have optimized the geometries of these clusters at the Hartree-Pock level using the 6-31 + + G(d,p) basis set. The impact of the electronic correlation, using the same basis set, was included using the second-order Moller-Plesset method and the Becke-Lee-Yang-Parr density functional, obtaining similar results from both methods. Our results indicate that in order to obtain reasonable values of the ionic dissociation energy, one has to use large clusters which saturate the first solvation shell of the HO2 molecule and that of the ionic products in an even way. The structures of the clusters are given and rationalized. Keywords:
Solvation; Ionic dissociation;
Hydroperoxyl
radical; Ab initio calculation
1. Introduction The ionic dissociation of the hydroperoxyl radical (H.02) into the H+ and 02 ions according to the following reaction HO:!tiH+
+O;
0)
is thought to be the key process in the decomposition of ozone in the atmosphere in the presence of cloud droplets [I]. The experimental data available on this reaction indicate that it is an equilibrium reaction, with equilibrium constant at 298K of 3.5 x lo-’ 121.
From this constant one obtains a AG between the dissociated and neutral forms in solution of 6.07 kcai mol-’ 121. At the same time, experimental measures concluded that the speed of the reactions involved in the equilibrium is slow. Why is this reaction important to the atmospheric chemistry of ozone? The answer comes from the fact that reaction (I) provides a route to generate 0; anions in the cloud droplets [2-41 and these anions once generated can undertake the attack of the ozone molecules present in the water droplets, according to the following reaction: 0,
* Corresponding author. ’ Dedicated to Professor .I. Bertrrin on his 65th birthday 0166-128O/Y6/$1.5.00
0 1996 Elsevier
PII SO1 66-l 280(96)04633-7
Science
+0~+~~0~2O~~OHiOH-
(2)
Reaction (2) is very fast: its rate constant amounts to 1.5 x lo9 mol 1-r s-r at 298K [Z], so it is a very
B.V. All rights reserved
144
C. Pirez de1 Valle et aLlJournal of Molecular Structure (Theochem) 371 (1996) 143-152
effective pathway for the decomposition of ozone. As a consequence, the equilibrium (1) is displaced to the product side and the concentration of HOa is reduced. This explains why the concentration of HO2 within the clouds is substantially lower than the value expected from its solubility equilibrium and that gaseous HO1 in cloud-free areas is larger than in cloudy ones. In summary, the existence of the combined reactions (1) + (2) in the clouds creates a mechanism for the decomposition of ozone in cloudy zones of the atmosphere, a mechanism which is capable of generating a pronounced depletion in the ozone levels compared to cloud-free regions [l]. From the previous paragraph, the importance of understanding the factors governing reaction (1) in aqueous media is clear. Preliminary computations carried out by us indicated that this reaction is highly endothermic in the gas phase, not surprisingly due to the ionic nature of the products. At the same time, the available experimental data indicate the presence of an equilibrium in which the ionic products lie just about 6 kcal mol-’ above the reactant. Therefore, the solvent plays an important role in stabilizing these products and making possible this reaction. A way of studying the role played by the solvent in reaction (1) is by carrying out ab initio methods on water clusters of the type H02.**(HZO), with n large enough to mimic the solvent effects. Within this approach reaction (1) is modeled as: H02...(H20),
R H’ .*(H20)1 + 0; . ..(H.O),
(3)
where 1 + m = n. Therefore, as a first step one has to know the optimum geometries and associated total energies of the H02.**(Hz0),, H30+...(H20)1-t and O;..(H,O),,, clusters. With these values in hand, one can compute the values of U for reaction (3) and compare it with the experimental value for AC. This will help to determine the size of the cluster needed to study the kinetics of the process, which should be our last step. In this work we will focus our attention in studying the variation of AE with n for various values of n, 1 and m that include the first solvation shell for the three systems. We will carry out our study using the Hartree-Fock and MP2 methods and, at the same time, with the BLYP functional to test the applicability of this density functional to the study of this
type of cluster. Finally, using a simplified model for the outer sphere of coordination of water molecules, we will make an extrapolation of our results to the bulk.
2. Computational details
The cluster geometries considered here were fully optimized at the Hartree-Fock (HF) level using the 631++G(d,p) basis set and the GAUSSIAN-@DF-~ program [5]. Knowing the total energies associated to each optimum geometry, we can obtain the values of the energy difference between products and reactants (hereafter identified by the symbol De). To have an estimate of the importance of the electron correlation on the values of D, we have carried out second order Moller-Plesset (MP2) computations [6] on the optimum Hartree-Fock geometries. As the MP2 computations become very time consuming for large clusters, we have also carried out density functional computations using the gradient corrected exchange functional of Becke [7] in conjunction with the gradient corrected correlation functional of LeeYang-Parr (BLYP) [8], as implemented in the GAUSSIAN-~~/DFT program [5]. It has been previously shown that the gradient corrected density functionals are capable of reproducing the known properties of hydrogen-bonded complexes [9-141. In particular, the BLYP functionals have been shown before to give results in good agreement with the MP2 ones for the structure and energetic properties of water clusters [13] and for the study of the ionization process of water molecules, acids and bases [l&16]. All the MP2 and BLYP computations carried out in this work were done on the ground state of these clusters (singlet for the H+*..(HzO)l cluster and doublet for the H02..(H20),, and O;.*.(H,O),,, clusters) using the 6-31 + + G(d,p) basis set.
3. Results and discussion We have started our study locating the minimum energy structures for the various X...(H20), clusters (n 5 4) at the Hartree-Fock level. We located the most likely starting structures using the computed molecular electrostatic potential (MEP) maps [17] of
C. Pirez de1 Valie et aLlJournal of Molecular Structure (Theochem) 371 (1996) 143-152
145
Fig. 1. Molecular electrostatic potential for the HO* (upper row, left), O-r ( u pp er row, center), HrO’ (upper row, right) and Hz0 (lower row) - 30 (deep dark), - 10 (shaded dark) and + 10 (light) molecules. For the HO2 molecule, we have represented three zones of potential: kcal mol.‘. For the 0, molecule, the zones of potential are: - 190 (deep dark) and - 110 (shaded dark) kcal mol-’ and the positive region (not drawn for simplicity) lies close to the oxygen atoms. For the HrO’ molecule, they are: + 120 kcal mol-t (light color), as no negative potential regions are found. Finally, for the Hz0 molecule the potentials represented are: - 55 (deep dark), - 15 (shaded dark) and 15 (light) kcal mol-‘.
HOz, O;, H30+ and HZ0 (see Fig. 1). These MEP maps indicate the zones in which the electrostatic interaction of a proton with the molecule is attractive (negative values of the energy, indicated in dark in our figures) and the zones in which the interaction is repulsive (positive values of the energy, in light in our figures). In deep dark we have marked in our figures the most stable region of the attractive zone, that is, the most stable region where electrostatically one can place the proton. The overlap of negative and positive regions of two molecules creates a favorable condition for the stability of its dimer, while the overlap of areas of the same character gives rise to unstable interactions from the electrostatic point of view. Within the X**.(H20), clusters one just finds O-H.--O bonds and this type of hydrogen bond has an interaction energy dominated by the electrostatic component. Therefore the MEP overlap analysis can
give a very good idea of the areas where the dimers, trimers and so on are expected to be electrostatically more stable. The analysis of the MEP map of the water molecule (Fig. 1) indicates that the molecule can act as a base (proton acceptor) in its oxygen and as an acid (proton donor) in the two H, in good agreement with the previous reported amphoterous behavior of this molecule. The HO2 map (Fig. 1) shows that this molecule can act as a base (proton acceptor) in its terminal oxygen and as an acid (proton donor) in its only H. At the same time, the central oxygen atom is weakly basic (proton acceptor). Therefore, the substitution of a hydrogen of the water molecule by a oxygen to produce the HO2 molecule modifies the MEP map around the central oxygen and the region of highest electrostatic potential, with a value of -35 kcal mol-‘, is located on the terminal oxygen. A
146
C. Pkrez de1 Valle et al.IJournal of Molecular Structure (Theochem) 371 (1996) 143-152
secondary local minima of -23 kcal mol-’ is also located on the same oxygen, similar in value to the local minima sitting on the central oxygen. The hydrogen is still in a repulsive region, as is the water molecule. These values indicate that the HO2 molecule can behave as an acid or as a base. The MEP map of the 0; anion (see Fig. 1) just shows positive potential in the regions far away from the nuclei, which are the important ones for our analysis. Therefore, this anion behaves as a base, with four strong proton acceptor points located on the position of the lone pairs for a sp* oxygen, that is, around 120” from the O-O axis. Also interesting is the relationship between the HO2 and 0; MEP maps. Qualitatively, the HO2 molecule is formed by adding a proton to one of the lone pairs of the 0; anion, thus eliminating one of the previously indicating proton attachment sites. The MEP map of the HO2 molecule should present three attachment sites located on the three positions of the 0; anion not involved in the protonation which generated the new O-H bond. This is the qualitative shape found for the MEP map of the HO2 molecule (see Fig. 1). The small changes in strength of the local minima just reflect the modifications induced in the wave function of the 0; anion by the addition of the proton (polarization, charge redistribution,. . .). The remaining MEP is that for the H30’ cation. It is nearly spherical and positive everywhere. In consequence, this compound is expected to behave only as a proton donor until all protons are attached to a base, in good agreement with previous computational and experimental studies. Using the information of the MEP maps and the overlap criteria described above, the addition of a water molecule to the HO2 molecule can take place with the water molecule acting as an acid or as a base. In the first case, the water molecule will try to create a O,-H.e.0 bond against the HO2 molecule and the most likely place where this bond is going to be formed is on the lateral oxygen atom, where the most negative area is located (see Fig. 1). Altematively, the O,-H.s.0 bond could also be formed with the central oxygen, but it is clear from the MEP map that this is not a situation as stable as the previous one. When the water molecule acts as a base, the only choice for the 00-H..*O,,, bond is the hydrogen end of the HOZ molecule. One can also carry out a
similar procedure to determine the best entry points for a second, third, or other number of water molecules. However, one must take into account that the addition of a water molecule can modify the shape of the MEP map of the isolated molecule due to the polarization and charge transfer induced in the molecule. One can expect to attach up to four molecules, one acting as base attached on the only hydrogen and three on the local minima located on the oxygens. However, giving that not all the attachment points have the same strength, one can expect a strong preference for the water molecules to go to the strongest regions and then to attach to that solvated molecule instead of attaching to the central oxygen. This means that one can find water clusters where the water molecules are mainly located in the extremes, creating the so-called “attached clusters” to distinguish them from the solvated clusters where the attached water molecules uniformly surround the solvated molecule. A similar analysis on the 0; anion indicates that one can expect up to four water molecules in the first solvation shell of this anion, with all O-H...0 bonds siting in the same plane. The HsO+ cation is expected to behave as a proton donor with the capability of forming up to three 0-H.e.0 bonds with the nearby water molecules. The previous overlap analysis, combin;ng the MEP maps and the overlap of regions, is just one step forward with respect to the so-called “Legon and Millen rules” [18]. These rules establish that one should expect the O-H*.*0 bonds to be formed in the direction of the lone pairs of the acceptor molecule. There is a qualitative agreement for our molecules between using the MEP maps and following Legon and Millen rules: for the 0; anion, where the two sp2 lone pairs sitting on each oxygen atom are located in exactly the same location where the minima of the MEP map are found and a similar behavior is found for the HO2 molecule. However, direct application of these rules does not justify why the H30’ cation is not acting as an acceptor on the lone pair located on the central oxygen, thus coordinating a fourth water molecule in its first coordination sphere. On the other hand, Legon and Millen rules, due to their qualitative nature, are not capable of predicting the relative stability of the different attachment sites or possible variations in the strength of each site due to the structural modifications.
C. Ptrez de1 Valle et aLlJournal of Molecular Structure (Theochem) 371 (1996) 143-152
Fig. 2. Optimum
energy structures
of the clusters HO,...(H,O),
Taking as starting geometries those obtained from the previous analysis, we have fully optimized, at the HF level and using the 6-31++G(d,p) basis set, the geometries of the H+..(H20)1 (I = l-4) and 02..(H20), (m = l-4 cluster, together with some of the corresponding n = 1 + m combinations of the H02.*.(H20)n cluster. The optimum geometries of these clusters are represented in Fig. 2, Fig. 3 and Fig. 4, while the corresponding total energies are collected in Table 1. The final optimum geometries for each cluster are similar to the starting ones obtained from the previous MEP analysis. As water molecules are added to the solvated molecule to build up the cluster, they progressively saturate the attachment points present in the MEP maps. As a general rule, each water molecule tries to make as many as possible O-H..*0 bonds with the solvated molecule, using as points of
Fig. 3. Optimum energy structures
141
(n = 1, 2, 3, 4).
attachment the ones predicted in the MEP analysis described above. This implies, in some cases, relaxing the optimum geometry of the O-H**-0 bonds from linearity. When enough water molecules are present, some 0-H. *.O contacts are against other water molecules, a fact that allows some relaxation in the strain imposed in each O-H...0 bond for a small number of water molecules (compare, for instance, the structures of 0;...(H20)2 and O;...(H20)3 in Fig. 3). Now, let us focus our attention on the energies collected in Table 1. We should start by saying that the hydration energies collected in it are similar to the previously reported theoretical [19] and experimental [20] data. The small changes between our values and the previous theoretical data just come from small variations in the geometries and basis sets employed. We have not corrected our values by the basis set superposition error effect, which previous studies
of the clusters H30’...(Hz0),
(n = 1, 2, 3).
148
C. Pbrez de1 Valle et al/Journal of Molecular Structure (Theochem) 371 (1996) 143-152
Fig. 4. Optimum energy structures
of the clusters 02...(HaO).
[19] have shown to be about 3 kcal mol-’ for the H+***(H20),, clusters, given the exploratory nature of this study. Here the emphasis is put on the model and the errors induced by the use of inadequate models are much larger than 3 kcal mol-‘. It is also interesting to note in Table 1 the similarity between the HP and MP2 values independently of the size of the cluster, indicting that the solute-solvent interaction is mainly electrostatic and correlation does not play an important role. The BLYP values are very close to the
(n - 1, 2, 3, 4).
MP2 results and, in consequence, follow the same trend. Using the values of total energy for each complex, and the associated hydration energy included in Table 1, one can deduce the most stable attachment site for the addition of a water molecule. Not surprisingly, one finds that the addition to the neutral molecule is less stable than the addition to the ionic fragments, independently of the number of water molecules. Between the two ions, Table 1 results indicate that it is always more stable to attach water
Table 1 Total energy (in au) and, in parenthesis, the associated hydration energy (in kcal mole’) computed H’..(HrO)., O;..(HrO), and HOr..(H20), clusters, as a function of n Cluster
n
c-w”
0
-76.031309
1
-152.070627
(5.02)
-152.475113
(6.29)
-152.837191
(5.45)
0. 1 2 3 4
0. -76.311196 -152.392565 -228.461143 -304.523659
(175.51) (31.39) (23.26) (19.58)
0. -76.507614 -152.796707 -229.069179 -305.335239
(172.61) (35.46) (25.04) (21.02)
0. -76.684954 -153.157301 -229.610086 -306.055788
(169.76) (36.43) (24.16) (19.72)
0 1 2 3 4
-149.594084 -225.655205 -301.712524 -377.765524 -453.816111
(18.69) (16.31) (13.60) (12.09)
-149.956396 -226.223409 -302.486053 -378.743185 -454.999688
(21.62) (18.88) (15.42) (15.03)
-150.339910 -226.787929 -303.229797 -379.668940 -456.104894
(21.17) (17.32) (15.61) (13.61)
0 1 2 3 4
-150.181981 -226.225231 -302.270015 -378.306226 -454.345946
(7.49) (8.45) (3.07) (5.27)
-150.519156 -226.767250 -303.018776 -379.258477 -455.500915
(9.76) (11.90) (4.49) (6.21)
-150.902118 -227.330293 -303.762304 -380.182330 -456.603997
(8.74) (11.14) (3.62) (4.65)
H+..(H20),
0;...(H20),
H0r..(H20),
HF
for the optimum energy structures of the
BLYP
MP2
-76.414250
-76.232540
C. Ptrez
del V&e et &Journal
ofMolecularStructure
Table 2 Dissociation energy (in kcal mole’) computed for the reaction o H’..$H20), + 0;...(H20)m as a function of H02...(Hz0), the number of molecules of water (n) n
1
In
0 1
0 1 0 2 1 0 3 2 1 0 4 3 2 1 0
0 0 1 0 1 2 0 1 2 3 0 1 2 3 4
2
3
4
10 22
D,(HF) 368.66 200.64 357.46 177.70 190.39 349.60 157.40 162.08 177.16 339.07 143.10 143.98 151.04 168.83 332.25 91.53 34.42
D,(MP2)
D,(BLYP)
352.90 190.16 341.04 166.60 180.45 334.07 146.05 149.47 166.06 323.14 131.24 130.64 136.80 156.85 314.31
352.56 191.53 340.12 166.24 181.49 333.93 145.70 148.68 167.80 321.94 130.62 129.17 136.02 156.84 312.98 82.59 23.23
moleculess to the H’ cation than to the 0; anion, being remarkable high the value of the first hydration of the proton to form the H30’ cation, similar to a typical covalent bond, a fact that has induced some authors to say that such a cation has three covalent bonds. Using Table 1 values, one can also compute the dissociation energy (De) for the ionic dissociation indicated in reaction (3) as a difference between the total energies of the fragments and the neutral starting cluster. We began by computing these values at the Hartree-Fock level for increasing values of n, the number of water molecules in the non-dissociated H02...(H20),, cluster. For a given value of n there are many possible combinations of the ionic H’...(HZO)~ and 0;***(H20)m clusters with the restriction that n = 1 + m. Table 2 collects the values of D, computed at the Hartree-Fock level for n I 4. Small values of D, indicate a higher stability of the dissociated forms, that is, a stronger solvent effect. According to Table 2 values, the most stable (I,m) option at the Hartree-Fock level is to attach all the water molecules to the H30’ cation. The (2,m) combinations in which the water molecules are equally distributed between the two ions, when available, do not lie energetically much higher. Finally, the least favorable (l,m)combination is that in which all
(Theo&em)
371 (1996)
143-152
149
water molecules are attached to the 0; anion. Even selecting the most stable combination we find positive values for D,, thus indicating that the ionic dissociation in these water clusters is always an endothermic process. Also interesting is the decrease in the values of D, as the number of water molecules increases: with four water molecules we have decreased the dissociation energy to less than half the value in the absence of water, an indication of the important cooperative effects present in this type of cluster. However, the value of the dissociation energy for the most stable of the largest clusters is still far from the 6 kcal mol-’ value estimated experimentally (see above). In order to fill this gap, we decided to investigate the effects of the inclusion of the electron correlation and the use of larger clusters. To understand the importance of the electron correlation on the values of the dissociation energy for reaction (3) we decided to compute the total energy at the MP2 level for each cluster at its optimum HartreeFock geometry. At the same time, as the computational cost becomes more and more prohibitive as the size of the cluster is increased, we tested the ability of the BLYP density functional to carry out the same task. The MP2 and BLYP values of D, computed at the optimum Hartree-Fock geometry are collected in the second and third columns of Table 2. These values show that the correlation energy decreases the values of D, in a near constant value and the difference (D ,(HF) - D e(MP2)) is always in the range between 9 and 17 kcal mol-‘. This is in good agreement with the mainly electrostatic character of the solute-water interactions in these clusters. These values also indicate that the inclusion of the electron correlation at the MP2 level is not going to fill the energy gap between the experimental and computed values. The same is valid for the BLYP values, because the MP2 and BLYP values do not differ by more than 2 kcal mol-’ in any case, giving support to the use of the BLYP functional on larger water clusters when we want to include the correlation energy. Finally, the relation of the relative stability of the various (Z,m) combinations indicate that the inclusion of the correlation does not modify the main trends reported at the Hartree-Fock level and the only change is in the relative order of the (I = 4, M = 0) and (I = 3, m = 1) dissociated clusters, at both the MP2 and BLYP levels.
150
C. Pirez de1 Valle et aLlJournal of Molecular Structure (Theochem) 371 (19%) 143-152
Given the small effect that the electron correlation has in decreasing the gap we decided to explore the effect of increasing the number of water molecules in these clusters. Although 12= 4 is already a large cluster for ab initio computations on the neutral cluster, it is not large enough to allow the formation of ionic clusters where the solvent effects are described with a similar degree of precision as in the neutral case. The reason is that while with 12= 4 one can complete the first solvation shell of the HO* molecule, we have not enough water molecules to saturate the first solvation shell of the two ionic fragments at the same time. This problem can be solved by increasing the value of n to at least the sum of the first shell of both ionic fragments. The first shell of the H’ ion is four (equivalent to three in H30’) [21], although there has been some controversy in the literature concerning the possibility of attaching one more molecule. For the 0, anion, the first shell can be taken as four in view of the above results. However, this anion still has nonattached lone pairs on the a electrons where more molecules could be attached. Therefore, we tested the possibility of having a water molecule attached in these positions by optimizing the structure of the O;--.(H20)6 cluster starting from the optimum geometry of the 0;-.(HZ0)4 cluster and adding two more water molecules, one above and another below the sp2 plane. Our computations showed the presence of a minimum at this position, whose HF energy is -605.912754 au, with the structure shown in Fig. 5. In this figure one finds that each water molecule is attached directly to one oxygen atoms of the 0; anion and to one of the water molecules of the first solvation shell, so these two new water molecules can be seen as part of the first or of the second solvation shell. For a better representation of the solvent effects on the 0; anion we decided to include these two molecules in our new studies with larger clusters designed to saturate the first solvation shells. We explored the effect of increasing the cluster size using a cluster of 10 water molecules, that is, a starting neutral HOz--*(HzO)to cluster dissociating into the H+-..(H20), and O;...(HZO)~ ionic clusters. To lower the computational cost of the test, the structure of the HOP* *(HzO) 10cluster used was the one obtained after an optimization at the semiempirical PM3 level [22] and shown in Fig. 6. The structure of the H+*+.(HzO)d cluster appears in Fig. 3 and that of
Fig. 5. Optimum energy structure of the cluster 02...(H20)6.
0;.*(H20)6 is in Fig. 5. The dissociation energy computed for the ionic dissociation of the H02**.(H20)t0 cluster is 91.53 kcal mol-’ at the HF level, and 82.59 kcal mol-’ at the BLYP level, a value much lower than the best obtained with the previous clusters and closer to the experimental result. Given the remarkable improvement obtained, we decided to estimate the impact of adding the water molecules of the outer shells of the HOZ-=-(H;?O) to, H’**(Hz0)4 and 0;*.(H20)6 clusters. We estimated this effect in an approximate form, representing the external shell water molecules just by point charges located in the positions where the H and 0 atoms of each water molecule would be sitting. This is the crudest approximation as it just includes part of the electrostatic effect. The positions of these atoms were
Fig. 6. Optimum energy structure of the cluster H0~...(H20)lo.
C. Pkrez de1 Valle et al.IJournal of Molecular Structure (Theochem) 371 (1996) 143-152
defined in such a way that a water molecule would be attached to all the non-bonded lone pairs of the clusters and each water molecule would be located in a tetrahedral environment as the tetrahedral coordination is the most important one in water solutions at all temperatures [23]. The O-H...0 distance is taken to be the one in the ice crystal. No further optimizations were carried out due to the approximate nature of our computation. The point charges used were + 0.423 on H and -0.846 on 0, obtained after an electrostatic potential fit of the water HF/6-31 + + G(d,p) wave function. The dissociation energy for reaction (3) after the inclusion of these point charges, simulating the effect of the second solvation shell, becomes 34 kcal mol-’ at the HF level and 23.23 kcal mol-’ at the BLYP level. The BLYP value is reasonably close to the experimental value, if one considers all approximations made in its obtainment and the fact that our D, values do not include the zero-point corrections due to the vibrational motions, which is known to make the computed D, value smaller for this type of clusters. Also the entropic factor could play an important role in decreasing the gap. One can expect that by including all these factors for the IZ= 10 computations on the fully optimized geometries one can obtain a result for the computed AG value reasonably close to the experimental value. Otherwise, one has to include the polarizing effect of the second outer shell water molecules, probably using a better model than the point charge one (such as a continuous model [24]). However, our study has shown that in order to properly describe the energetics of the ionic dissociation in water solution, reaction (3) one has to use large and properly selected clusters designed to saturate the first solvation shell of the two ions and the neutral molecule, and that the agreement is improved if one adds an estimate of the field created by the second solvation shell. The results obtained from this study will help to find better model clusters to describe the energetics of reaction (3) and/or model clusters designed for the study of the kinetics of this reaction in solution.
Acknowledgements This work is dedicated to Professor J. Bertran on his 65th birthday. It reflects our appreciation for his
151
advice and support in the difficult times, as well as his sense of humor in the good times. This research was partially supported by the Human Capital and Mobility Program, Access to Large Installations, under contract CHGE-(X92-0009 “Access to supercomputing facilities for European researchers” established between The European Community and CESCAKEPBA. This work was also supported by the DGICYT and CIRIT under Projects PB-920655-CO2-02 and GRQ94-1077. One of us (C.P.V.) also thanks the Eusko Jaurlaritza for his doctoral grant (grant number BF91.140).
References [l] J.E. Johnson and ISA. Isaksen, J. Atmos. Chem., 16 (1993) 99. (21 J. Lelieveld and P.J. Crutzen, Nature, 343 (1990) 227. [3] SE. Schwartz, Nature, 343 (1990) 209. [4] P. Warneck, Ber. Bunsenges. Phys. Chem., 96 (1992) 454. [5] GAusstAN-92, REV. ~2, M.J. Frisch, G.W. Trucks, M. HeadGordon, P.M.W. Gill, M.W. Wong, J.B. Foresman, B.G. Johnson, H.B. Schlegel, M.A. Robb, E.S. Replogle, R. Gomperts, J.L. Andres, K. Raghavachari, J.S. Binkley, C. Gonzalez, R.L. Martin, D.J. Fox, D.J. Defrees, J. Baker, J.J.P. Stewart and J.A. Pople, Gaussian, Inc., Pittsburgh, PA, 1992. [6] C. Moller and M.S. Plesset, Phys. Rev., 46 (1934) 618; R. Krishnan, M.J. Frish and J.A. Pople, J. Chem. Phys., 72 (1980) 4244. [7] A.D. Becke, Phys. Rev. A, 38 (1988) 3098; A.D. Becke, J. Chem. Phys., 96 (1992) 2155. [8] C. Lee, W. Yang and R.G. Parr, Phys. Rev. B, 37 (1988) 785. [9] C. Lee, H. Chen and G. Fitzgerald, J. Chem. Phys., 102 (1995) 1266; C. Lee, H. Chen and G. Fitzgerald, J. Chem. Phys.. 101 (1994) 4472. [lo] F. Sim, A. St-Amant, I. Papai and D.R. Salahub, J. Am. Chem. Sot., 114 (1992) 4391. [ll] K. Laasonen, M. Parrinello, R. Car, C. Lee and D. Vanderbilt, Chem. Phys. Lett., 207 (1993) 208. [12] K. Kim and K.D. Jordan, J. Phys. Chem., 98 (1994) 10089. [13] S.S. Xantheas, J. Chem. Phys., 102 (1995) 4505. [14] J.J. Novoa and C. Sosa, J. Phys. Chem., 99 (1995) 15837. [15] C. Lee, C. Sosa and J.J. Novoa, J. Chem. Phys., 103 (1995) 4360. [16] C. Lee, C.Sosa, M. Planas and J.J. Novoa, J. Chem. Phys., submitted for publication. [17] E. Scrocco and J. Tomasi, Adv. Quant. Chem., 11 (1978) 115; P. Politzer and J.S. Murray, in K.B. Lipkowitz and D.B. Boyd (Eds.), Reviews in Computational Chemistry. Vo1.2. VCH, New York, 1991, Chapter 7. [18] A.C. Legon and D.J. Millen, Chem. Rev., 6 (1986) 635.
152
C. Pirez de1 Valle et aLlJournal of Molecular Structure (Theochem) 371 (19%) 143-152
[19] D. Peeters and G. Leroy, J. Mol. Struct. (Theochem), 209 (1990) 263; E.P.F. Lee and J.M. Dyke, Mol. Phys., 73 (1991) 375. [20] Y.K. Lau, S. Ikuta and P. Kebarle, J. Am. Chem. Sot., 104 (1982) 1462; M. MeorNer and C.V. Speller, J. Phys. Chem., 90 (1986) 6616. [21] J.O. Lundgren and J.M. Williams, J. Chem. Phys., 58 (1973) 788; R. Trio10 and A.H. Narten, J. Chem. Phys., 63 (1975) 3624; M.D. Newton, J. Chem. Phys., 67 (1978) 5535; E. Kochanski, J. Am. Chem. Sot., 107 (1985) 7869;
A. Rahmouni, E. Kochanski, R. Wiest, P.E.S. Wormer and J. Langlet, J. Chem. Phys., 93 (1990) 6648. [22] J.J. Stewart, J. Comput. Chem. 10 (1989) 221; M. Zerner, in K.B. Lipkowitz and D.B. Boyd (Eds.), Reviews in Computational Chemistry, VoI.2, VCH, New York, 1991, Chapter 8. [23] E. Clementi and G. Corongiu, in Simulated Water Structure, ACS Symposium Series 568, The American Chemical Society, Washington, 1994. [24] J. Tomasi and M. Persico, Chem. Rev., 94 (1994) 2027, and references cited therein.
Journal of Molecular
Structure (Theochem)
371 (1996) 143-152
The determinant role of water in the ionic dissociation of HOZ’ Carlos P&ez de1 Valle”, Carmela Valdemoro”,
Juan J. Novoah-*
~Institutode Matemcitica Aplicada y Fish Fundamental, CSIC Cl. Serrano 123, 28~-~udr~d, Spain hDept. Quimica Fisica, Fat. Quimica, Univ. de Barcelona Av. Diagonal 647, 08028-Barcelona. Spain Received
Abstract The energetics of the ionic dissociation of the hydroperoxyl radical (HOZ) into the H’ and 0; ions in solution is studied using ab initio methods and water clusters to mimic the solvent effects. The clusters selected are HO,..(H,O),, H’..(H@), and 0;..+(H20),, using various combinations of n, 1 and m designed to saturate the first solvation shell of the involved systems, the largest of which contains 10 water molecules (n = 10) plus the electrostatic field created by the outer solvation shell to these 10 water molecules. We have optimized the geometries of these clusters at the Hartree-Pock level using the 6-31 + + G(d,p) basis set. The impact of the electronic correlation, using the same basis set, was included using the second-order Moller-Plesset method and the Becke-Lee-Yang-Parr density functional, obtaining similar results from both methods. Our results indicate that in order to obtain reasonable values of the ionic dissociation energy, one has to use large clusters which saturate the first solvation shell of the HO2 molecule and that of the ionic products in an even way. The structures of the clusters are given and rationalized. Keywords:
Solvation; Ionic dissociation;
Hydroperoxyl
radical; Ab initio calculation
1. Introduction The ionic dissociation of the hydroperoxyl radical (H.02) into the H+ and 02 ions according to the following reaction HO:!tiH+
+O;
0)
is thought to be the key process in the decomposition of ozone in the atmosphere in the presence of cloud droplets [I]. The experimental data available on this reaction indicate that it is an equilibrium reaction, with equilibrium constant at 298K of 3.5 x lo-’ 121.
From this constant one obtains a AG between the dissociated and neutral forms in solution of 6.07 kcai mol-’ 121. At the same time, experimental measures concluded that the speed of the reactions involved in the equilibrium is slow. Why is this reaction important to the atmospheric chemistry of ozone? The answer comes from the fact that reaction (I) provides a route to generate 0; anions in the cloud droplets [2-41 and these anions once generated can undertake the attack of the ozone molecules present in the water droplets, according to the following reaction: 0,
* Corresponding author. ’ Dedicated to Professor .I. Bertrrin on his 65th birthday 0166-128O/Y6/$1.5.00
0 1996 Elsevier
PII SO1 66-l 280(96)04633-7
Science
+0~+~~0~2O~~OHiOH-
(2)
Reaction (2) is very fast: its rate constant amounts to 1.5 x lo9 mol 1-r s-r at 298K [Z], so it is a very
B.V. All rights reserved
144
C. Pirez de1 Valle et aLlJournal of Molecular Structure (Theochem) 371 (1996) 143-152
effective pathway for the decomposition of ozone. As a consequence, the equilibrium (1) is displaced to the product side and the concentration of HOa is reduced. This explains why the concentration of HO2 within the clouds is substantially lower than the value expected from its solubility equilibrium and that gaseous HO1 in cloud-free areas is larger than in cloudy ones. In summary, the existence of the combined reactions (1) + (2) in the clouds creates a mechanism for the decomposition of ozone in cloudy zones of the atmosphere, a mechanism which is capable of generating a pronounced depletion in the ozone levels compared to cloud-free regions [l]. From the previous paragraph, the importance of understanding the factors governing reaction (1) in aqueous media is clear. Preliminary computations carried out by us indicated that this reaction is highly endothermic in the gas phase, not surprisingly due to the ionic nature of the products. At the same time, the available experimental data indicate the presence of an equilibrium in which the ionic products lie just about 6 kcal mol-’ above the reactant. Therefore, the solvent plays an important role in stabilizing these products and making possible this reaction. A way of studying the role played by the solvent in reaction (1) is by carrying out ab initio methods on water clusters of the type H02.**(HZO), with n large enough to mimic the solvent effects. Within this approach reaction (1) is modeled as: H02...(H20),
R H’ .*(H20)1 + 0; . ..(H.O),
(3)
where 1 + m = n. Therefore, as a first step one has to know the optimum geometries and associated total energies of the H02.**(Hz0),, H30+...(H20)1-t and O;..(H,O),,, clusters. With these values in hand, one can compute the values of U for reaction (3) and compare it with the experimental value for AC. This will help to determine the size of the cluster needed to study the kinetics of the process, which should be our last step. In this work we will focus our attention in studying the variation of AE with n for various values of n, 1 and m that include the first solvation shell for the three systems. We will carry out our study using the Hartree-Fock and MP2 methods and, at the same time, with the BLYP functional to test the applicability of this density functional to the study of this
type of cluster. Finally, using a simplified model for the outer sphere of coordination of water molecules, we will make an extrapolation of our results to the bulk.
2. Computational details
The cluster geometries considered here were fully optimized at the Hartree-Fock (HF) level using the 631++G(d,p) basis set and the GAUSSIAN-@DF-~ program [5]. Knowing the total energies associated to each optimum geometry, we can obtain the values of the energy difference between products and reactants (hereafter identified by the symbol De). To have an estimate of the importance of the electron correlation on the values of D, we have carried out second order Moller-Plesset (MP2) computations [6] on the optimum Hartree-Fock geometries. As the MP2 computations become very time consuming for large clusters, we have also carried out density functional computations using the gradient corrected exchange functional of Becke [7] in conjunction with the gradient corrected correlation functional of LeeYang-Parr (BLYP) [8], as implemented in the GAUSSIAN-~~/DFT program [5]. It has been previously shown that the gradient corrected density functionals are capable of reproducing the known properties of hydrogen-bonded complexes [9-141. In particular, the BLYP functionals have been shown before to give results in good agreement with the MP2 ones for the structure and energetic properties of water clusters [13] and for the study of the ionization process of water molecules, acids and bases [l&16]. All the MP2 and BLYP computations carried out in this work were done on the ground state of these clusters (singlet for the H+*..(HzO)l cluster and doublet for the H02..(H20),, and O;.*.(H,O),,, clusters) using the 6-31 + + G(d,p) basis set.
3. Results and discussion We have started our study locating the minimum energy structures for the various X...(H20), clusters (n 5 4) at the Hartree-Fock level. We located the most likely starting structures using the computed molecular electrostatic potential (MEP) maps [17] of
C. Pirez de1 Valie et aLlJournal of Molecular Structure (Theochem) 371 (1996) 143-152
145
Fig. 1. Molecular electrostatic potential for the HO* (upper row, left), O-r ( u pp er row, center), HrO’ (upper row, right) and Hz0 (lower row) - 30 (deep dark), - 10 (shaded dark) and + 10 (light) molecules. For the HO2 molecule, we have represented three zones of potential: kcal mol.‘. For the 0, molecule, the zones of potential are: - 190 (deep dark) and - 110 (shaded dark) kcal mol-’ and the positive region (not drawn for simplicity) lies close to the oxygen atoms. For the HrO’ molecule, they are: + 120 kcal mol-t (light color), as no negative potential regions are found. Finally, for the Hz0 molecule the potentials represented are: - 55 (deep dark), - 15 (shaded dark) and 15 (light) kcal mol-‘.
HOz, O;, H30+ and HZ0 (see Fig. 1). These MEP maps indicate the zones in which the electrostatic interaction of a proton with the molecule is attractive (negative values of the energy, indicated in dark in our figures) and the zones in which the interaction is repulsive (positive values of the energy, in light in our figures). In deep dark we have marked in our figures the most stable region of the attractive zone, that is, the most stable region where electrostatically one can place the proton. The overlap of negative and positive regions of two molecules creates a favorable condition for the stability of its dimer, while the overlap of areas of the same character gives rise to unstable interactions from the electrostatic point of view. Within the X**.(H20), clusters one just finds O-H.--O bonds and this type of hydrogen bond has an interaction energy dominated by the electrostatic component. Therefore the MEP overlap analysis can
give a very good idea of the areas where the dimers, trimers and so on are expected to be electrostatically more stable. The analysis of the MEP map of the water molecule (Fig. 1) indicates that the molecule can act as a base (proton acceptor) in its oxygen and as an acid (proton donor) in the two H, in good agreement with the previous reported amphoterous behavior of this molecule. The HO2 map (Fig. 1) shows that this molecule can act as a base (proton acceptor) in its terminal oxygen and as an acid (proton donor) in its only H. At the same time, the central oxygen atom is weakly basic (proton acceptor). Therefore, the substitution of a hydrogen of the water molecule by a oxygen to produce the HO2 molecule modifies the MEP map around the central oxygen and the region of highest electrostatic potential, with a value of -35 kcal mol-‘, is located on the terminal oxygen. A
146
C. Pkrez de1 Valle et al.IJournal of Molecular Structure (Theochem) 371 (1996) 143-152
secondary local minima of -23 kcal mol-’ is also located on the same oxygen, similar in value to the local minima sitting on the central oxygen. The hydrogen is still in a repulsive region, as is the water molecule. These values indicate that the HO2 molecule can behave as an acid or as a base. The MEP map of the 0; anion (see Fig. 1) just shows positive potential in the regions far away from the nuclei, which are the important ones for our analysis. Therefore, this anion behaves as a base, with four strong proton acceptor points located on the position of the lone pairs for a sp* oxygen, that is, around 120” from the O-O axis. Also interesting is the relationship between the HO2 and 0; MEP maps. Qualitatively, the HO2 molecule is formed by adding a proton to one of the lone pairs of the 0; anion, thus eliminating one of the previously indicating proton attachment sites. The MEP map of the HO2 molecule should present three attachment sites located on the three positions of the 0; anion not involved in the protonation which generated the new O-H bond. This is the qualitative shape found for the MEP map of the HO2 molecule (see Fig. 1). The small changes in strength of the local minima just reflect the modifications induced in the wave function of the 0; anion by the addition of the proton (polarization, charge redistribution,. . .). The remaining MEP is that for the H30’ cation. It is nearly spherical and positive everywhere. In consequence, this compound is expected to behave only as a proton donor until all protons are attached to a base, in good agreement with previous computational and experimental studies. Using the information of the MEP maps and the overlap criteria described above, the addition of a water molecule to the HO2 molecule can take place with the water molecule acting as an acid or as a base. In the first case, the water molecule will try to create a O,-H.e.0 bond against the HO2 molecule and the most likely place where this bond is going to be formed is on the lateral oxygen atom, where the most negative area is located (see Fig. 1). Altematively, the O,-H.s.0 bond could also be formed with the central oxygen, but it is clear from the MEP map that this is not a situation as stable as the previous one. When the water molecule acts as a base, the only choice for the 00-H..*O,,, bond is the hydrogen end of the HOZ molecule. One can also carry out a
similar procedure to determine the best entry points for a second, third, or other number of water molecules. However, one must take into account that the addition of a water molecule can modify the shape of the MEP map of the isolated molecule due to the polarization and charge transfer induced in the molecule. One can expect to attach up to four molecules, one acting as base attached on the only hydrogen and three on the local minima located on the oxygens. However, giving that not all the attachment points have the same strength, one can expect a strong preference for the water molecules to go to the strongest regions and then to attach to that solvated molecule instead of attaching to the central oxygen. This means that one can find water clusters where the water molecules are mainly located in the extremes, creating the so-called “attached clusters” to distinguish them from the solvated clusters where the attached water molecules uniformly surround the solvated molecule. A similar analysis on the 0; anion indicates that one can expect up to four water molecules in the first solvation shell of this anion, with all O-H...0 bonds siting in the same plane. The HsO+ cation is expected to behave as a proton donor with the capability of forming up to three 0-H.e.0 bonds with the nearby water molecules. The previous overlap analysis, combin;ng the MEP maps and the overlap of regions, is just one step forward with respect to the so-called “Legon and Millen rules” [18]. These rules establish that one should expect the O-H*.*0 bonds to be formed in the direction of the lone pairs of the acceptor molecule. There is a qualitative agreement for our molecules between using the MEP maps and following Legon and Millen rules: for the 0; anion, where the two sp2 lone pairs sitting on each oxygen atom are located in exactly the same location where the minima of the MEP map are found and a similar behavior is found for the HO2 molecule. However, direct application of these rules does not justify why the H30’ cation is not acting as an acceptor on the lone pair located on the central oxygen, thus coordinating a fourth water molecule in its first coordination sphere. On the other hand, Legon and Millen rules, due to their qualitative nature, are not capable of predicting the relative stability of the different attachment sites or possible variations in the strength of each site due to the structural modifications.
C. Ptrez de1 Valle et aLlJournal of Molecular Structure (Theochem) 371 (1996) 143-152
Fig. 2. Optimum
energy structures
of the clusters HO,...(H,O),
Taking as starting geometries those obtained from the previous analysis, we have fully optimized, at the HF level and using the 6-31++G(d,p) basis set, the geometries of the H+..(H20)1 (I = l-4) and 02..(H20), (m = l-4 cluster, together with some of the corresponding n = 1 + m combinations of the H02.*.(H20)n cluster. The optimum geometries of these clusters are represented in Fig. 2, Fig. 3 and Fig. 4, while the corresponding total energies are collected in Table 1. The final optimum geometries for each cluster are similar to the starting ones obtained from the previous MEP analysis. As water molecules are added to the solvated molecule to build up the cluster, they progressively saturate the attachment points present in the MEP maps. As a general rule, each water molecule tries to make as many as possible O-H..*0 bonds with the solvated molecule, using as points of
Fig. 3. Optimum energy structures
141
(n = 1, 2, 3, 4).
attachment the ones predicted in the MEP analysis described above. This implies, in some cases, relaxing the optimum geometry of the O-H**-0 bonds from linearity. When enough water molecules are present, some 0-H. *.O contacts are against other water molecules, a fact that allows some relaxation in the strain imposed in each O-H...0 bond for a small number of water molecules (compare, for instance, the structures of 0;...(H20)2 and O;...(H20)3 in Fig. 3). Now, let us focus our attention on the energies collected in Table 1. We should start by saying that the hydration energies collected in it are similar to the previously reported theoretical [19] and experimental [20] data. The small changes between our values and the previous theoretical data just come from small variations in the geometries and basis sets employed. We have not corrected our values by the basis set superposition error effect, which previous studies
of the clusters H30’...(Hz0),
(n = 1, 2, 3).
148
C. Pbrez de1 Valle et al/Journal of Molecular Structure (Theochem) 371 (1996) 143-152
Fig. 4. Optimum energy structures
of the clusters 02...(HaO).
[19] have shown to be about 3 kcal mol-’ for the H+***(H20),, clusters, given the exploratory nature of this study. Here the emphasis is put on the model and the errors induced by the use of inadequate models are much larger than 3 kcal mol-‘. It is also interesting to note in Table 1 the similarity between the HP and MP2 values independently of the size of the cluster, indicting that the solute-solvent interaction is mainly electrostatic and correlation does not play an important role. The BLYP values are very close to the
(n - 1, 2, 3, 4).
MP2 results and, in consequence, follow the same trend. Using the values of total energy for each complex, and the associated hydration energy included in Table 1, one can deduce the most stable attachment site for the addition of a water molecule. Not surprisingly, one finds that the addition to the neutral molecule is less stable than the addition to the ionic fragments, independently of the number of water molecules. Between the two ions, Table 1 results indicate that it is always more stable to attach water
Table 1 Total energy (in au) and, in parenthesis, the associated hydration energy (in kcal mole’) computed H’..(HrO)., O;..(HrO), and HOr..(H20), clusters, as a function of n Cluster
n
c-w”
0
-76.031309
1
-152.070627
(5.02)
-152.475113
(6.29)
-152.837191
(5.45)
0. 1 2 3 4
0. -76.311196 -152.392565 -228.461143 -304.523659
(175.51) (31.39) (23.26) (19.58)
0. -76.507614 -152.796707 -229.069179 -305.335239
(172.61) (35.46) (25.04) (21.02)
0. -76.684954 -153.157301 -229.610086 -306.055788
(169.76) (36.43) (24.16) (19.72)
0 1 2 3 4
-149.594084 -225.655205 -301.712524 -377.765524 -453.816111
(18.69) (16.31) (13.60) (12.09)
-149.956396 -226.223409 -302.486053 -378.743185 -454.999688
(21.62) (18.88) (15.42) (15.03)
-150.339910 -226.787929 -303.229797 -379.668940 -456.104894
(21.17) (17.32) (15.61) (13.61)
0 1 2 3 4
-150.181981 -226.225231 -302.270015 -378.306226 -454.345946
(7.49) (8.45) (3.07) (5.27)
-150.519156 -226.767250 -303.018776 -379.258477 -455.500915
(9.76) (11.90) (4.49) (6.21)
-150.902118 -227.330293 -303.762304 -380.182330 -456.603997
(8.74) (11.14) (3.62) (4.65)
H+..(H20),
0;...(H20),
H0r..(H20),
HF
for the optimum energy structures of the
BLYP
MP2
-76.414250
-76.232540
C. Ptrez
del V&e et &Journal
ofMolecularStructure
Table 2 Dissociation energy (in kcal mole’) computed for the reaction o H’..$H20), + 0;...(H20)m as a function of H02...(Hz0), the number of molecules of water (n) n
1
In
0 1
0 1 0 2 1 0 3 2 1 0 4 3 2 1 0
0 0 1 0 1 2 0 1 2 3 0 1 2 3 4
2
3
4
10 22
D,(HF) 368.66 200.64 357.46 177.70 190.39 349.60 157.40 162.08 177.16 339.07 143.10 143.98 151.04 168.83 332.25 91.53 34.42
D,(MP2)
D,(BLYP)
352.90 190.16 341.04 166.60 180.45 334.07 146.05 149.47 166.06 323.14 131.24 130.64 136.80 156.85 314.31
352.56 191.53 340.12 166.24 181.49 333.93 145.70 148.68 167.80 321.94 130.62 129.17 136.02 156.84 312.98 82.59 23.23
moleculess to the H’ cation than to the 0; anion, being remarkable high the value of the first hydration of the proton to form the H30’ cation, similar to a typical covalent bond, a fact that has induced some authors to say that such a cation has three covalent bonds. Using Table 1 values, one can also compute the dissociation energy (De) for the ionic dissociation indicated in reaction (3) as a difference between the total energies of the fragments and the neutral starting cluster. We began by computing these values at the Hartree-Fock level for increasing values of n, the number of water molecules in the non-dissociated H02...(H20),, cluster. For a given value of n there are many possible combinations of the ionic H’...(HZO)~ and 0;***(H20)m clusters with the restriction that n = 1 + m. Table 2 collects the values of D, computed at the Hartree-Fock level for n I 4. Small values of D, indicate a higher stability of the dissociated forms, that is, a stronger solvent effect. According to Table 2 values, the most stable (I,m) option at the Hartree-Fock level is to attach all the water molecules to the H30’ cation. The (2,m) combinations in which the water molecules are equally distributed between the two ions, when available, do not lie energetically much higher. Finally, the least favorable (l,m)combination is that in which all
(Theo&em)
371 (1996)
143-152
149
water molecules are attached to the 0; anion. Even selecting the most stable combination we find positive values for D,, thus indicating that the ionic dissociation in these water clusters is always an endothermic process. Also interesting is the decrease in the values of D, as the number of water molecules increases: with four water molecules we have decreased the dissociation energy to less than half the value in the absence of water, an indication of the important cooperative effects present in this type of cluster. However, the value of the dissociation energy for the most stable of the largest clusters is still far from the 6 kcal mol-’ value estimated experimentally (see above). In order to fill this gap, we decided to investigate the effects of the inclusion of the electron correlation and the use of larger clusters. To understand the importance of the electron correlation on the values of the dissociation energy for reaction (3) we decided to compute the total energy at the MP2 level for each cluster at its optimum HartreeFock geometry. At the same time, as the computational cost becomes more and more prohibitive as the size of the cluster is increased, we tested the ability of the BLYP density functional to carry out the same task. The MP2 and BLYP values of D, computed at the optimum Hartree-Fock geometry are collected in the second and third columns of Table 2. These values show that the correlation energy decreases the values of D, in a near constant value and the difference (D ,(HF) - D e(MP2)) is always in the range between 9 and 17 kcal mol-‘. This is in good agreement with the mainly electrostatic character of the solute-water interactions in these clusters. These values also indicate that the inclusion of the electron correlation at the MP2 level is not going to fill the energy gap between the experimental and computed values. The same is valid for the BLYP values, because the MP2 and BLYP values do not differ by more than 2 kcal mol-’ in any case, giving support to the use of the BLYP functional on larger water clusters when we want to include the correlation energy. Finally, the relation of the relative stability of the various (Z,m) combinations indicate that the inclusion of the correlation does not modify the main trends reported at the Hartree-Fock level and the only change is in the relative order of the (I = 4, M = 0) and (I = 3, m = 1) dissociated clusters, at both the MP2 and BLYP levels.
150
C. Pirez de1 Valle et aLlJournal of Molecular Structure (Theochem) 371 (19%) 143-152
Given the small effect that the electron correlation has in decreasing the gap we decided to explore the effect of increasing the number of water molecules in these clusters. Although 12= 4 is already a large cluster for ab initio computations on the neutral cluster, it is not large enough to allow the formation of ionic clusters where the solvent effects are described with a similar degree of precision as in the neutral case. The reason is that while with 12= 4 one can complete the first solvation shell of the HO* molecule, we have not enough water molecules to saturate the first solvation shell of the two ionic fragments at the same time. This problem can be solved by increasing the value of n to at least the sum of the first shell of both ionic fragments. The first shell of the H’ ion is four (equivalent to three in H30’) [21], although there has been some controversy in the literature concerning the possibility of attaching one more molecule. For the 0, anion, the first shell can be taken as four in view of the above results. However, this anion still has nonattached lone pairs on the a electrons where more molecules could be attached. Therefore, we tested the possibility of having a water molecule attached in these positions by optimizing the structure of the O;--.(H20)6 cluster starting from the optimum geometry of the 0;-.(HZ0)4 cluster and adding two more water molecules, one above and another below the sp2 plane. Our computations showed the presence of a minimum at this position, whose HF energy is -605.912754 au, with the structure shown in Fig. 5. In this figure one finds that each water molecule is attached directly to one oxygen atoms of the 0; anion and to one of the water molecules of the first solvation shell, so these two new water molecules can be seen as part of the first or of the second solvation shell. For a better representation of the solvent effects on the 0; anion we decided to include these two molecules in our new studies with larger clusters designed to saturate the first solvation shells. We explored the effect of increasing the cluster size using a cluster of 10 water molecules, that is, a starting neutral HOz--*(HzO)to cluster dissociating into the H+-..(H20), and O;...(HZO)~ ionic clusters. To lower the computational cost of the test, the structure of the HOP* *(HzO) 10cluster used was the one obtained after an optimization at the semiempirical PM3 level [22] and shown in Fig. 6. The structure of the H+*+.(HzO)d cluster appears in Fig. 3 and that of
Fig. 5. Optimum energy structure of the cluster 02...(H20)6.
0;.*(H20)6 is in Fig. 5. The dissociation energy computed for the ionic dissociation of the H02**.(H20)t0 cluster is 91.53 kcal mol-’ at the HF level, and 82.59 kcal mol-’ at the BLYP level, a value much lower than the best obtained with the previous clusters and closer to the experimental result. Given the remarkable improvement obtained, we decided to estimate the impact of adding the water molecules of the outer shells of the HOZ-=-(H;?O) to, H’**(Hz0)4 and 0;*.(H20)6 clusters. We estimated this effect in an approximate form, representing the external shell water molecules just by point charges located in the positions where the H and 0 atoms of each water molecule would be sitting. This is the crudest approximation as it just includes part of the electrostatic effect. The positions of these atoms were
Fig. 6. Optimum energy structure of the cluster H0~...(H20)lo.
C. Pkrez de1 Valle et al.IJournal of Molecular Structure (Theochem) 371 (1996) 143-152
defined in such a way that a water molecule would be attached to all the non-bonded lone pairs of the clusters and each water molecule would be located in a tetrahedral environment as the tetrahedral coordination is the most important one in water solutions at all temperatures [23]. The O-H...0 distance is taken to be the one in the ice crystal. No further optimizations were carried out due to the approximate nature of our computation. The point charges used were + 0.423 on H and -0.846 on 0, obtained after an electrostatic potential fit of the water HF/6-31 + + G(d,p) wave function. The dissociation energy for reaction (3) after the inclusion of these point charges, simulating the effect of the second solvation shell, becomes 34 kcal mol-’ at the HF level and 23.23 kcal mol-’ at the BLYP level. The BLYP value is reasonably close to the experimental value, if one considers all approximations made in its obtainment and the fact that our D, values do not include the zero-point corrections due to the vibrational motions, which is known to make the computed D, value smaller for this type of clusters. Also the entropic factor could play an important role in decreasing the gap. One can expect that by including all these factors for the IZ= 10 computations on the fully optimized geometries one can obtain a result for the computed AG value reasonably close to the experimental value. Otherwise, one has to include the polarizing effect of the second outer shell water molecules, probably using a better model than the point charge one (such as a continuous model [24]). However, our study has shown that in order to properly describe the energetics of the ionic dissociation in water solution, reaction (3) one has to use large and properly selected clusters designed to saturate the first solvation shell of the two ions and the neutral molecule, and that the agreement is improved if one adds an estimate of the field created by the second solvation shell. The results obtained from this study will help to find better model clusters to describe the energetics of reaction (3) and/or model clusters designed for the study of the kinetics of this reaction in solution.
Acknowledgements This work is dedicated to Professor J. Bertran on his 65th birthday. It reflects our appreciation for his
151
advice and support in the difficult times, as well as his sense of humor in the good times. This research was partially supported by the Human Capital and Mobility Program, Access to Large Installations, under contract CHGE-(X92-0009 “Access to supercomputing facilities for European researchers” established between The European Community and CESCAKEPBA. This work was also supported by the DGICYT and CIRIT under Projects PB-920655-CO2-02 and GRQ94-1077. One of us (C.P.V.) also thanks the Eusko Jaurlaritza for his doctoral grant (grant number BF91.140).
References [l] J.E. Johnson and ISA. Isaksen, J. Atmos. Chem., 16 (1993) 99. (21 J. Lelieveld and P.J. Crutzen, Nature, 343 (1990) 227. [3] SE. Schwartz, Nature, 343 (1990) 209. [4] P. Warneck, Ber. Bunsenges. Phys. Chem., 96 (1992) 454. [5] GAusstAN-92, REV. ~2, M.J. Frisch, G.W. Trucks, M. HeadGordon, P.M.W. Gill, M.W. Wong, J.B. Foresman, B.G. Johnson, H.B. Schlegel, M.A. Robb, E.S. Replogle, R. Gomperts, J.L. Andres, K. Raghavachari, J.S. Binkley, C. Gonzalez, R.L. Martin, D.J. Fox, D.J. Defrees, J. Baker, J.J.P. Stewart and J.A. Pople, Gaussian, Inc., Pittsburgh, PA, 1992. [6] C. Moller and M.S. Plesset, Phys. Rev., 46 (1934) 618; R. Krishnan, M.J. Frish and J.A. Pople, J. Chem. Phys., 72 (1980) 4244. [7] A.D. Becke, Phys. Rev. A, 38 (1988) 3098; A.D. Becke, J. Chem. Phys., 96 (1992) 2155. [8] C. Lee, W. Yang and R.G. Parr, Phys. Rev. B, 37 (1988) 785. [9] C. Lee, H. Chen and G. Fitzgerald, J. Chem. Phys., 102 (1995) 1266; C. Lee, H. Chen and G. Fitzgerald, J. Chem. Phys.. 101 (1994) 4472. [lo] F. Sim, A. St-Amant, I. Papai and D.R. Salahub, J. Am. Chem. Sot., 114 (1992) 4391. [ll] K. Laasonen, M. Parrinello, R. Car, C. Lee and D. Vanderbilt, Chem. Phys. Lett., 207 (1993) 208. [12] K. Kim and K.D. Jordan, J. Phys. Chem., 98 (1994) 10089. [13] S.S. Xantheas, J. Chem. Phys., 102 (1995) 4505. [14] J.J. Novoa and C. Sosa, J. Phys. Chem., 99 (1995) 15837. [15] C. Lee, C. Sosa and J.J. Novoa, J. Chem. Phys., 103 (1995) 4360. [16] C. Lee, C.Sosa, M. Planas and J.J. Novoa, J. Chem. Phys., submitted for publication. [17] E. Scrocco and J. Tomasi, Adv. Quant. Chem., 11 (1978) 115; P. Politzer and J.S. Murray, in K.B. Lipkowitz and D.B. Boyd (Eds.), Reviews in Computational Chemistry. Vo1.2. VCH, New York, 1991, Chapter 7. [18] A.C. Legon and D.J. Millen, Chem. Rev., 6 (1986) 635.
152
C. Pirez de1 Valle et aLlJournal of Molecular Structure (Theochem) 371 (19%) 143-152
[19] D. Peeters and G. Leroy, J. Mol. Struct. (Theochem), 209 (1990) 263; E.P.F. Lee and J.M. Dyke, Mol. Phys., 73 (1991) 375. [20] Y.K. Lau, S. Ikuta and P. Kebarle, J. Am. Chem. Sot., 104 (1982) 1462; M. MeorNer and C.V. Speller, J. Phys. Chem., 90 (1986) 6616. [21] J.O. Lundgren and J.M. Williams, J. Chem. Phys., 58 (1973) 788; R. Trio10 and A.H. Narten, J. Chem. Phys., 63 (1975) 3624; M.D. Newton, J. Chem. Phys., 67 (1978) 5535; E. Kochanski, J. Am. Chem. Sot., 107 (1985) 7869;
A. Rahmouni, E. Kochanski, R. Wiest, P.E.S. Wormer and J. Langlet, J. Chem. Phys., 93 (1990) 6648. [22] J.J. Stewart, J. Comput. Chem. 10 (1989) 221; M. Zerner, in K.B. Lipkowitz and D.B. Boyd (Eds.), Reviews in Computational Chemistry, VoI.2, VCH, New York, 1991, Chapter 8. [23] E. Clementi and G. Corongiu, in Simulated Water Structure, ACS Symposium Series 568, The American Chemical Society, Washington, 1994. [24] J. Tomasi and M. Persico, Chem. Rev., 94 (1994) 2027, and references cited therein.