J. Phys. Chem. 1993,97,5547-5552
5547
Proton Solvation in Liquid Water. An ab Initio Study Using the Continuum Model Iiiaki T&6n and Estanislao Silla' Department de Qulmica-Rsica, Universitat de Valbncia, 461 00 Burjassot (Valbncia), Spain
Juan B e r t h Departament de Qulmica. Universitat Autbnoma de Barcelona, 08193 Bellaterra (Barcelona), Spain Received: October 13, I992
The proton solvation process is studied by means of a continuum model using a quantum description of the solute charge distribution. The hydronium ion has been selected as the starting point for all our calculations. The ion and the ion plus a first solvation shell have been used as solutes in the continuum calculations, and the results have beem compared with a discrete model with two explicit solvation shells. According to various experimental studies, we have found a well-defined first solvation shell with three water molecules, while the preferred location of a fourth water molecule is in the second shell, both in the gas phase and in solution. Although the combined discrete-continuum model gives an energy value in good agreement with experiment, some deficiencies, related to the noninclusion of charge transfer beyond the cavity, have been found.
Introduction Proton behavior in solution plays an important role in many chemical and biological systems, and it has been one of the more interesting problems since the early days of physical chemistry. Proton solvation in water can be seen as a particular case of a proton-transfer reaction, one of the most important elementary reactions in chemistry. The large negative free energy associated with thisstepdetermines theacid-baseequilibria in water solution and thus is of particular interest for elucidating many catalytic and enzymatic reaction mechanisms. This is why many experimental and theoretical works have dealt with both the energetics and the structure of proton solvation. The precision with which the proton solvation free energy can be determined experimentally is limited because the standard hydrogen potential cannot be obtained by measurement alone. Several values appear in the literature, one of the most recent measurements being that of Reiss and Heller.' They give an estimate of the absolute potential of the hydrogen electrode equal to-4.43 V. Thisvalueisidenticalwith that calculated by Gurevich and Pleskov2 and very similar to those of Lohmann3 (-4.48 V) and Trasatti4 (-4.44 V). Gomer and Trysod estimated a value of -4.73 V. These five independent measurements lead to a proton solvation free energy ranging from -254 to -261 kcal/mol with a mean value of -259.5 kcal/mol.6 This is very similar to values given before such as that of Noyes7 (-260.5 kcal/mol). Besides this uncertainty, it seems clear, since its existence was postulated in 1907,a that the largest part of this energy comes from the H30+formation. In fact, in the formation of the hydronium ion, an enthalpy of -166.3 kcal/mol and a free energy of -1 58.9 kcal/ mol9 have been measured. This represents more than the 60% of the total solvation free energy. The rest of the energy comes from interaction with other water molecules forming the first, second, and other solvation shells. The structure and coordination number of the hydronium ion have also been the object of some controversy. Accurate quantum mechanical calculationsI0have confirmed the nonplanar nature of this ion with a HOH angle close to 112O, although basis sets without polarization functions give a planar geometry. Regarding the coordination number, it has been shown by means of cluster studies," molar volume measurements,'* and rate of proton dissociation determinations13 that the structure of theequilibrated hydrated protonI4could be H904+(i.e., with three water molecules in the first solvation shell) which retains a tetrahedral configuration. However, some diffraction's and X-ray studies16 show 0022-3654/93/2097-5547$04.00/0
the possibility of four nearest-neighbor H20 molecules for each H30+, the orientation of the fourth water molecule having two possibilities: proton donor or charge-dipole oriented. Much attention has been devoted to the theoretical studies of the protonated water clusters H+(H20),. Small clusters, with n = 2 or 3, have been extensively studied by means of quantum mechanic^.'^-^^ Quantum mechanical studies have also been carried out for n = 4-528-32 (Le., completing at least the first minimum solvation shell). These studies generally agree well with experimental energy values in the gas phase, but the number of water molecules used in the calculations is limited by the capability of the present computers. In these studies, it was found that the preferred location of the fourth water molecule is in the second solvation shell. In one of this papers,2aNewton considered, in a limited search, the possibility of four water molecules in the first shell. Although a minimum was not found, he concluded that if space-filling requirements demand a fourth ligand, it must be charge-dipole oriented. By using a more complete basis set for the central hydronium ion, Karlstrom3*found a minimum for the fourth water molecule in the first solvation shell, with the oxygen atom toward the central ion and with a stabilizationenergy of about 4 kcal/mol. Others studies, such as nonquantum mechanical methods, deal with a much larger number of water molecule^^^-^^ (for example, the case of n = 21, which presents an unusual stability, has been examined by means of these techniques). Molecular dynamics (MD) has been used to deal with the ionic equilibrium of water.35 In this work, both the coordination numbers 3 and 4 were considered for the H30+ion. Monte Carlo (MC) techniques can alsodeal with larger H+(H20),~ystems.~~J~ In a workconceming clusters with n = 1-9,36 Kochanski found the possibility of three or four water molecules in the first solvation shell, depending on the number of water molecules included in the simulation. This was also found by Fornili et al.37in a MC study of a hydronium ion and 215 water molecules. However, in none of these studies were the three-body forces explicitly considered, and it has been shown that they play an important role, especially for the first The inclusion of the three-body potential solvation she11.3a~39 (polarization and exchange repulsion) usually reduces the coordination number. Thus, the two main problems concerning proton solvation (total free energy and coordination number) can be studied both by means of a quantum treatment of the correspondingsupermolecule and nonquantum methods (MD and MC). The problem of a 0 1993 American Chemical Society
5548 The Journal of Physical Chemistry, Vol. 97, No. 21, 1993
TuAdn et al.
quantum study arises from the number of solvent molecules that can be explicitly considered, while in MD and MC methods the potential function used is the critical point. One strategy for solving the problems of both the inclusion of the bulk solvent effects (Le., the consideration of the effect of a large number of solvent molecules) and the consideration of a complete interaction potential (Le., with three and higher body forces) is to make use of cavity models combined with a quantum mechanicaltreatment of the solute. In this work, we have used a combined discretecontinuum model for the proton solvation. We have used the hydronium ion as the starting point in the study of proton solvation. First, we have investigated the structure of the first solvation shell and the number of solvent molecules that can be placed in this shell. Second, we have modeled the rest of the solvent with a continuum model by using both the hydronium ion and the hydronium ion plus the first solvation shell as solutes. The results are compared with a full discrete model with two solvation shells of the solvent molecules.
(-0.764) -0.755
(-0.924) 4.869 H2/01
\
0.434 (0.462)
HI
I
f
H
84
Computational Details The ab initio molecular orbital calculations in the gas phase have been carried out at HartreeFock and second-order of Maller-Plesset perturbation theory levels with the 6-3 1G*@basis set. Geometries have been optimized at these levels with Berny's algorithm4' by means of the GAUSSIAN8842 package of programs. All hydronium-water interaction energies were corrected for basis set superposition error (BSSE) using the counterpoise correction proposed by Boys and Bernardi.43 In order to calculate free energies, zero-point vibrational energies, thermal corrections, and entropy changes were considered. Harmonic vibrationalfrequencies and entropy calculationswithin the ideal gas-rigid rotor-harmonic oscillatormodel were computed analytically at the Hartree-Fock level and numerically at the MP2 level at the optimized geometries using standard mechanical form~lae.~% For~ 5the H+ species, we have used the value of 26.012 eu for the entropy in the gas phase46calculated for the translational and electronic entropy. The bulk solvent effects have been taken into account by means of the cavity model developed by Rivail and co-workers.4749 In this approach, the solvent is assimilated to a continuum characterized by the dielectricpermittivity which surrounds the cavity were the solute is placed. The electrostatic interaction for the free energy of solvation is given by
.
m
+I
where ( M p )is a component of the multipole moment of order 1 and ( R p ) is the corresponding component of the reaction field. These components can be analytically calculated for spherical and ellipsoidal cavities. In practice, a good convergence is obtained extending the sum upto I = 6. This electrostatic interaction term can be included into a SCF method by introducing the corresponding operator into the Hamiltonian. Moreover,the analytical energy derivatives of the electrostatic term have been recently developed,so leading to an efficient geometry optimization procedure.51 In our case, we have used ellipsoidal cavities surrounded by a continuum of dielectric constant equal to that of water, i.e., 78.4. The volume of the initial ellipsoidal cavity is constrained to be equal to the average molecular volume in the liquid. The molecular volume is obtained by means of an empirical equation relating thevan der Waals and the molecular volume.5*However, as has been shown p r e v i o ~ s l ythis , ~ ~approach ~ ~ ~ is no longer valid when one is dealing with large solute-solvent clusters. In this case, the ellipsoidalcavity used is the minimal ellipsoid that ensures that the van der Waals surface remains inside. In both cases, the axes of the cavity are related to the axes of inertia of a solid of
H 0.487
1.896
Figure 1. Structures studied in this work and their atomic Mulliken charges, in the gas phase and in solution (in parentheses).
uniform density limited by the van der Waals surface.51 These definitions allow the use of deformable cavities during the geometry optimization procedure. In the case of the water molecule, we have tested that both cavity definitions lead to nearly the same final geometry in solution. Two more terms are needed in order to complete the energetic description of the solvation free energy: the dispersion and the cavitation terms. Thedispersion term has been approximated by means of an atom-atom potential in the framework of the continuum model of the ~ o l v e n t . Finally, ~ ~ * ~ ~the cavitation term has been computed with the Pierotti formulae derived from the scaled particle theory.56 All solvent effect calculations were carried out at 298.15 K. Re5ultS Figure 1 shows the structures studied in this work. Structures 111, IV, and V give the three possibilities for the first solvation shell studied here. Structure VI has been included for the purposes of comparison with IV and V. Finally, structure VI1 has been studied as a model of solute plus two solvation shells. From structure I to V, the study has been carried out both in the gas phase and in a continuum, while structures VI and VI1 (where a second solvation shell is considered) have been studied only in the gas phase. In this section, we present first the energy results, because they give us an idea about the quality of the methodology used. In the second place, the nuclear geometries and electronic structures are discussed. Energetics. In Table I, we present the potential and free energiesforthereactionsH++ n(H2O)-H+(H20),,(n = 1-10). Values have been obtained at the HF/6-3 lG* and MP2/6-3 lG* levels, and the hydronium-water interaction energies were corrected for basis set superpositionerror. BSSE representsabout
The Journal of Physical Chemistry, Vol. 97, No. 21, 1993 5549
Proton Solvation in Liquid Water
-
TABLE I: Potential and Free Energies for the Berrctions H+ nH 0 H(H*O),,+Calculated at the HF/631G* and MP2f631G* Levels and ExDerimental Values (in kcal/moW
+
~~
~
HF/6-31G* structure I1 I11 IV V VI VI1
AE AG -174.82 -160.27 -250.75 -206.84 -254.92 -204.21 no minimum found -263.74 -211.78 -314.77 -218.58
TABLE II: Continnum Contributions (Electrostatic plus Induction, Dispersion, and Cavitation) to the Solvation Free Energy for Several Structures (in kcaVmol)'
~~
~
MP2/6-31G*
AE AG -174.68 -159.54 -257.57 -209.73 -263.43 no minimum found -271.75
structure
I I1
AG(~X~)~ -158.9 -205.8
I11 IV V
-211.4
Values are referred to 1 atm and 298.15 K.
5% of the interaction energy at the HF level and 7% at the MP2 level in agreement with previous work.3* Experimental results are taken from ref 9. Other experimentalvalues have been offered, but they are very similar to those used here." Frequency calculations at the MP2 level have only been carried out for the smaller clusters (n = 1 and 3) because of the computational cost. For H+(H20)lo (structure VII), it has not been possible to make MP2 calculations because of its size and the zero-point energy, and entropy calculations at the HF level were made by using frequencies obtained for H+(H2O)5 (structure VI), H+(H20)6, and H+(H20),. It can be seen that the HF free energies are very close to the experimental values. This results from a cancellation of errors between the calculated interaction energy and the calculated entropy. The interaction energies and the loss of entropy are both overestimated. These two trends nearly cancel each other out, giving a final free energy very close to the experimentalvalues. MP2 interaction energies are always greater than the HF values because of the inclusion of the dispersion forces. Special attention must be paid to structures IV, V, and VI (i.e., with five solvent molecules). The absolute minimum, both in free and in potential energy, is given by structure VI. However, structure IV is also stationary in the potential energy hypersurface(PES) but not in the freeenergy hypersurface (FES). Finally, structure V does not present any minima either in the FES or in the PES. If we consider free energies, both structures IV and V dissociate into structure 111 and a water molecule. Considering the possibility of four water molecules in the first solvation shell, we can compare the energetic differencesbetween structures IV and VI. Structure VI is 7.6 kcal/mol more stable than IV when considering free energies and 8.8 kcal/mol when considering potential energies. This example shows the importance of taking into account the zero-point energy and entropy effects on the formation of clusters. Thus, it seems that the hydronium ion prefers the first shell with three water molecules and the next water molecules placed in the second shell. These results are in agreement with the previous work of Newtonz8 about the coordination of a central H30+ ion. The estimated free energy for the structure VI1 also shows the importance of the entropic and zero-point terms. The large negative interaction potential energy is, in great part, compensated for the loss of entropy and zero-point energy. The electrostatic plus induction (AG,l) and nonelectrostatic (AG,,, and AGdls) contributions to the solvation energies have been calculated in the framework of the continuum model as discussed in the previous section. Geometries have been reoptimized in the cavity (Le., nuclear and electronic relaxationsunder the effect of the reaction field have been considered). Variations in the zero-point vibrational energies, thermal corrections, and entropieswith respect to the in vacuo values53have been included intotheelectrostatic term. Calculationswiththecontinuummodel have been carried out at the HF/6-31G* level. Table I1 shows the electrostatic, dispersion, and cavitation components of the solvation energies of structures I-V. As was expected, the electrostatic contribution grows quickly when going from H2O to H30+ and diminishes when we include water molecules explicitly (structure 111) because we are increasing the volume
AGcav -3.34 4.32 -3.59 5.42 -9.79 12.39 no minimum found no minimum found &is
-5.76 -78.07 -54.84
Values are referred to the 1M
-
G o t
4.78 -76.23 -52.25
1M process.
TABLE IIk Total Solvation Free Energies Obtained as the Sum of Discrete and Continuum Contributions' structure I1 111 VI1
&iserste
-160.27 -206.84 -218.58
AGcontb -74.34 -50.36
nAGvapb 2.89 11.56 28.90
G o t
-231.72 -245.64 -1 89.68
- -
The free energy of vaporization of n water molecules is included to close the thermodynamic cycle. Final values refer to the 1 atm 1 M 1 atm. process. This includes the expansion work from 1 M @
of the cavity maintaining the same total charge. It must also be pointed out that the possibility of four water molecules in the first solvation shell disappears completely. No minimum has been found for structures IV and V when geometry optimization is carried out mside a cavity. Geometry optimization of structure IV surrounded by a continuum results in dissociationinto structure I11 and a water molecule. This is not an unexpected result since insolution wearedealing withfreeenergies. Byusingacontinuum model that includes a repulsive interaction term with the cavity, Karlstrom32 found a minimum for structure IV. To ensure that this dissociation of structure IV is not an artifact due to the cavity we have defined, we have also carried out an intermolecular distance optimization with the PCM where the cavity is defined by a set of interlocking spheres.59.60 With this model, we have also found that structure IV is no longer a minimum. Moreover, this fact is consistent with the features of structure IV invacuo and can be explainedby looking at the Mulliien charge# over the atoms appearing on Figure 1. It can be seen that there exists a zero charge transfer between the fourth water molecule and the H(HzO)3+ ion. Thus, we merely have an electrostatic interaction between a dipole and an ion that is weaker when a continuum of high dielectric constant surrounds them. As structure IV was not a strong minimum, the presence of a high polarizable medium makes this minimum disappear. In fact, we are substituting a discrete electrostatic interaction of a dipole with structure I11 by an electrostatic interaction of a continuum with structure 111. We have also used the PCM model to compare the results obtained for structure 111. Within this model, we have obtained an electrostaticcontribution to the solvationfreeenergy of -55.27 kcal/mol, very close to that given by the ellipsoidal cavity model that we are using. Moreover, the conclusions obtained about the changes in the geometry and electronic structure when passing from the gas phase to solution are very similar in both models. It seems reasonable that in species with a high degree of symmetry, as structure I11 is, the choice of cavity shape is not as decisive as it could be in other species. In Table 111the total free energies of solvation, calculated as the sum of continuum and discrete contributions, are shown. To close the thermodynamic cycle, the vaporization free energy of n water molecules is included. To be consistent, all values used are calculated at the HF/6-31G* level. Structure I1 represents a pure continuum model of the solvation of the H30+ion and structure I11 a continuum-discrete model, while structure VI1 represents a pure discrete model with the two first solvation shells. The pure continuum model gives a good value (Le., 90% of the total energy), if we consider that the volume of the cavity used does not take into account that the positive charge will produce
SSSQ The Journal of Physical Chemistry, Vol. 97, No. 21, 1993
Tuil6n et al.
TABLE IV Main Distances (in A) for the Studied Structures in the Cas Phase (e = 1) and in Solution (e = 78.4) at the HF/631C9 Level’ structure e R(OlH1) R(02H) R( 0102) R(0105) R(0205) R(0208) I
1 78.4 1 78.4 1 78.4 1 78.4 1 1
I1 I11 IV
VI VI1 0
0.947 (0.969) 0.948 0.969 (0.991) 0.968 0.986 (1.016) 0.987 0.984 (1.015)
0.950 (0.971) 0.953 0.950 (0.972)
2.619 (2.593) 2.636 2.628 (2.606)
1.002 (1.041) 0.994
0.963 (0.990) 0.956
2.546 (2.513) 2.574
3.270 (3.220) 2.778 (2.720) 2.860
2.873
MP2/6-31Gt distances in the gas phase are given in parentheses.
a contraction in the cavity because of the large attraction for the solvent molecules, leading to a greater solvation energy. This point will be discussed further. The mixed discrete-continuum model, with the first complete solvation shell, gives an energy veryclosttotheexperimentalone (95%). Finally, thethirdmodel, in which the second shell is well-defined but the bulkelectrostatic contributions are not considered, gives the worst approximation to the experimental value. This was an expected result since the interaction of a point charge with a dielectric is proportional to 1/ r ,so a large number of water molecules are needed to saturate the ionic part of the solvation energy. Geometry. In Table IV, the main bond lengths as defined in Figure 1 are shown both in the gas phase (at the H F and MP2 levels) and in solution (at the H F level) for those structures for which a minimum was found. In the gas phase, the MP2 and H F distances do not differ very much. The OH distances are slightly larger at the MP2 level, and the hydronium-water or waterwater distances are slightly shorter as could be expected from the inclusion of the dispersion energies. In any case, they follow the same tendencies and as said are very similar. This also supports the decision to use the H F level in the continuum calculations. For the water molecule, the continuum produces a lengthening of the OH bond as in the discrete representation of the water dimer.5’ Although the continuum model does not take into account explicitly hydrogen bond formation, inasfar as the main component in this bond arises from the electrostatic interaction, it seems that the continuum simulates the formation of hydrogen bonds of the water molecules with the rest of the solvent molecules. In contrast, for the hydronium cation, the presence of a continuum medium produces a shortening of the OH bond, while, as can be seen in structure 111, the formation of hydrogen bonds lengthens this bond. This fact can be explained assuming that the continuum “sees” a global 1charge whose solvation is favored as the cavity diminishes. As will be discussed below, the final effect depends on the cavity size used. In structure 111, the continuum produces a lengthening of all molecular distances. The case of the 0 2 H bond is the same as in the water molecule (structure I) but increased because of the polarization induced by the positive charge. This is the first illustration of the cooperative effect. The O l H l bond suffers a little lengthening. It has been previously pointed out that the OH bonds are weakened by the electronic polarization solvent effect.5’ This fact suggests that the hydrogen bonds are somewhat stronger in the liquid. The 0 1 0 2 distance is lengthened because both thedipole and the ion are better solvated if they remain apart, increasing the solute-solvent interface. All the aspects discussed until now are confirmed by the force constants given in Table V. Returning to Table IV, we can see that for structure IV, the 0 1 0 2 distance is smaller than the 0 1 0 5 distance by about 0.5 A. This aspect also confirms the existence of a well-defined first shell with only three water molecules. The fourth water molecule is bonded only by electrostatic interactions without any chargetransfer process, as has been discussed before. As has been pointed out in the previous section, geometry optimization of structure IV in solution leads to structure I11 and a water molecule.
+
Force C ~ ~ t r r n(in t s au) for the OH Bonds of Some studied structures structure e OlHl 02H I 1 0.613 78.4 0.608 I1 1 0.532 78.4 0.534 I11 1 0.414 0.605 78.4 0.404 0.585 TABLE V
TABLE VI: Mulliken Charges (in au) for Different Croups of the Studied Structures at the HF/6-31C+ Level’ structure e I w., ws WC WR ~
I11 IV VI VI1
1 78.4 1 78.4 1 1
~~
~
0.834 0.816 0.838
0.055 0.061 0.054
0.055 0.061 0.054
CO.001
0.825 0.796
0.038 0.018
0.052 0.018
0.034 0.025
~~
0.025
I stands for the central H30+ion and W, for the water molecules, where i indicates the number of the oxygen as given in Figure 1.
The geometric changes produced in structure VI are the same that those of structure VI1 where the second shell is completed. So, we will discuss them both referring to structure VII. However, it is interesting to notice that the 0 1 0 2 distance diminishes when one water molecule is added to the second shell (attached to the water molecule that contains the 0 2 atom), but it increases when water molecules of the second shell are added to other water molecules of the fist shell. Inclusion of an explicit second solvation shell produces a lengthening of the O l H l and 0 2 H bonds, as the continuum produces on the O l H l bond of structure 111, but the effect is now substantially larger. It seems that dielectric calculations underestimate the effect of solvation in water but give the right trends. However, the 0 1 0 2 distance behaves completely differently, because in this case the formation of hydrogen bonds with the second shell produces a shortening of the distance, while in the continuum model the opposite happens. The distance obtained in the discrete model with two solvation shells (2.57 A) compares better with the experimental distance (2.52 f 0.02 A7)than that obtained in the mixed continuumdiscrete model (2.64 A). This different behavior between the continuum and the discrete models will be discussed below. Electronic Structure, The electronic distribution of the structures studied here appears in Figure 1 where the atomic Mulliken charges, both in the gas phase ( a = 1) and in solution (c = 78.4) are given. In Table VI, the Mulliken charges are given over the H30+ion and the different water molecules. It can be seen, from the atomic Mulliken charges in the gas phase and in solution, that thewater molecule (structure I) suffers an important increment of the atomic charges under the effect of a polarizable continuum. This fact, together with the nuclear relaxation,s’ produces an increment of the dipolar moment of the water molecule up to 2.4 D. The hydronium ion also suffersan electronic polarization, as can be concluded comparing atomic charges in the gas phase and in solution. However, this effect is significantly
Ptoton Solvation in Liquid Water lower than in the case of the water molecule which is also related to the size of the cavity used as happened in the geometry optimization. In structure 111, the presence of a continuum produces a greater delocalization of the positive charge over the surrounding water molecules as reflected in Table VI. The new atomicchargedistributionleads to a strengtheningof the hydrogen bond, coherent with the lengthening of the O l H l bond. The peripheric water molecules show a similar net polarization process as that of the isolated water molecule (structure I). For structure IV, as mentioned before, no significant charge-transfer process is observed toward the last water molecule (W,in Table VI), although from the atomic charges of Figure 1 it can be deduced that this water molecule is polarized by the presence of the H+(H20)4 ion. As in the case of geometry optimizations, the electronic distribution of structure VI presents the same main new features that were produced in structure VI1 in which the second shell is completed. It is interesting to compare the electronic distribution of structure I11 in a continuum with structure VII. The positive charge on the central H30+ ion diminishes in both cases with respect to structure I11 in the gas phase, but the effect is more pronounced in structure VII. The delocalization of the central charge is 3 times more intense in the second shell than in the first shell. Evidently, the continuum model does not permit a charge transfer beyond the cavity interface, and this effect does not appear. Thus, it seems that with the first solvation and a continuum,we can describe correctly the electronic distribution of the central ion but not that of the first shell. This fact is closely related to the opposite trends in the 0 1 0 2 distances between the mixed discrete-continuum and the pure discrete model. As can be seen in Figure 1, the atomic charge on the oxygen atom of the first-shell water molecules of structure VI1 is larger than that in structure I11 both in the gas phaseand in solution, while the H1 atomiccharge hardly changes. This fact leads to a stronger hydrogen bond in structure VI1 reflected in the lengtheningof the O l H l bond and the shortening of the 0 1 0 2 distance. In structure 111, although the presence of a continuum increases the atomic charge of the oxygen atom of the water molecules which produces a strengthening of the hydrogen bond (the O l H l bond is lengthened), this effect is not strong enough to compensate the trend of the dielectricto separate an ion and a dipole to solvate them better.
Discussion In this work, we have used three different models to study the solvationof the hydronium ion: a pure continuum model, a mixed discrete-continuum model, where the first solvationshell has been explicited, and a pure discrete model with two solvation shells. Each of them presents some advantages in its use. Thus, in the continuummodel, the specific interactions are not well-represented but the long-range electrostatic interactions are considered, while in a discrete model with a reduced number of solvent molecules, the last are not included. The main advantage of the pure continuum model (i.e., a H3O+ ion inside a cavity) is its simplicity. This model can inform us quickly about some features such as electronic polarization and geometricchanges induced by the solvent and about the magnitude of the solvation free energy. However, this model depends radically on the cavity size used. In this work, we have used for the hydronium ion a cavity size derived from a relationship established for neutral species, while a contraction of the cavity is expected for charged ones.48 In this way, we have obtained a good energy estimation (about go%), but the geometric changes (shortening of the OH bonds) are just the opposite from those expected (lengtheningof the OH bonds simulating the formation of hydrogen bonds with the solvent). In fact, by scaling down the cavity size by a factor of 0.1, a lengthening of 0.002 A is obtained for the OH bonds, and the solvation energy is closer to the experimentalone. However, it is difficult to find an objective
The Journal of Physical Chemistry, Vol. 97, No. 21, 1993 5551 way to define the cavity size for charged species although some possibilities have been recently proposed based on Mulliken charges62 or in ab initio studies of solute-solvent complexes.63 The mixed discrete-continuummodel, where the first solvation shell is explicitly included in the solute definition, seems to be a more plausible way to study the solvation of the hydronium ion and other charged species.53 The inclusion of the first shellmakes this model not so dependent on the cavity size selection. This model gives a better energy estimation and reproduces many of the geometric and electronic changes induced by the solvent, but when compared with structure VII, it seems that it underestimates these changes. Moreover, it has been shown how this model predicts a lengthening of the 0 1 0 2 distance, while when we use a discrete representation of the second solvation shell, this 0 1 0 2 distance is shorter. This fact has been shown to be related to the charge-transfer process to the second solvation shell that is not represented in our discrete-continuummodel. It must be pointed that the use of the smaller basis set could lead to erroneous conclusionsabout the goodness of this model for proton solvation. Thus, by using a 3-21G basis set, a 0 1 0 2 distance of 2.50 A is obtained for structure I11 in the gas phase. While this distance is shorter than the experimental distance determined in solution (2.52 0.02 A), one could conclude that the lengthening of this distance obtained in the continuum calculations (+0.015 A) is correctly representing the solvent effects. We have shown, comparing with an ab initio calculation including the second shell, that this could not be so. We have also used the pure discrete model with two solvation shells that has been included for comparative purposes. Evidently, this model is not able to give a good estimation of the solvation energy of a charged species, because the interaction of a charge with a dielectric only decreases very slowly with distance. However, this model correctly predicts the experimental intermolecular distances in solution and is very useful to compare the predictions of the continuum or discrete-continuum models. Besides the discussion of the validity of the different models used, what is evidenced in our calculations is the importance of the nonadditive effects on both the geometric and the electronic changes in the structures studied. This fact is very important in the development of new interaction functions for Monte Carlo or molecular dynamics simulationsof the hydronium ion or other ions in liquid water. Regarding the geometric aspects, it has been shown how the cooperative effects determine the response of the OH bond of the water molecule to the solvent and the intermolecular ion-water distances. Regarding the electronic structure, our results show the importance of the charge transfer in describing the hydrogen bonds. The nonadditivity of this term, together with that arising from polarization, should be taken into account in the development of new theoretical potential functions. Then, although the continuum approach could be used to derive new potentials because the solute polarization is included, caution must be exercised in those cases where charge transfer could be important. Thus, in a pure continuum representation of the hydronium in solution,it has a +1 charge, but when we introduce the first solvation shell, this charge diminishes to 0.830. The consideration of a second solvation shell reduces the charge in 0.153 au because of the charge transfer to this new shell.
*
Conclusions In this work, we have investigated the structure and the energetics of the solvated proton by using a combined quantum continuum approach where a solute is inscribed in a cavity surrounded by a continuum. Various definitions of solute have been used. Regarding the structure, our calculations,both in the gas phase and in solution, show that the hydronium ion in solution has a well-definedfirst solvationshell with three water molecules,while the preferred location of the fourth water moleculeis in the second
5552 The Journal of Physical Chemistry, Vol. 97, No. 21, 1993
Tufi6n et al.
(19) Meyer, W.; Jakubetz, W.; Schuster, P. Chem. Phys. Lett. 1973,21, shell. Entropic and electrostatic considerations support this conclusion,whichisinagreement withseveralpreviousworks.11-*4 97.(20) Delpuech, J. J.; Serratrice, G.; Strich, A.; Veillard, A. Mol. Phys. A good description of the intermolecular ion-water distance is 1975, 29, 849. (21) Kochanski, E. Nouv. J . Chim. 1984,8, 605. only obtained when two solvation shells are explicitlyconsidered. (22) Yamabe, S.; Minato, T.; Hirao. K. J. Chem. Phys. 1984,80, 1576. The geometric changes produced in the presence of a continuum (23) Del Bene, J. E.; Frisch, M. J.; Pople, J. A. J . Phys. Chem. 1985,89, are in agreement with the two solvation shells model except in 3669. the ion-water distance. The noninclusion of the charge transfer (24) Frisch, M. J.; Del Bene, J. E.; Binkley, J. S.; Schaefer, H. F.,111. J . Chem. Phys. 1986, 84, 2279. beyond the cavity seems to be the reason for this discrepancy. In (25) Del Bene, J. E. J . Comput. Chem. 1987,8, 810. spite of this, the mixed discrete-continuum model gives the best (26) Del Bene, J. E. J . Phys. Chem. 1988, 92, 2874. estimation of the solvation free energy, with a value that represents (27) Scheiner, S.; Yu, W. 0. Int. J. Quantum Chem. 1991, 18, 37. (28) Newton, M. D. J . Chem. Phys. 1977,67, 5535. around 95% of the experimental energy. Moreover, this value (29) Deakyne, C. A.; Meot-Ner, M.; Campbell, C. L.; Hughes, M. G . ; is not so dependent on the cavity size as in the pure continuum Murphy, S. P. J. Chem. Phys. 1986, 84,4958. model, making this result more objective. (30) Lee, E. P. F.;Dyke, J. M.; Wilders,A. E.; Watts,P.Mol. Phys. 1990, Finally, the importance of the nonadditive terms in the 71, 207. (31) Peeters, D.; Leroy, G.THEOCHEM 1990, 209, 263. description of solute behavior has also been evidenced in (32) Karlstrom, G.J . Phys. Chem. 1988,92, 1318. accordance with previous work.51 Geometric changes, together (33) Holland, P. M.; Castelman, A. W. J. Chem. Phys. 1980, 72, 5984. with the polarization and the charge-transfer terms, should be (34) Searcy, J. Q.; Fenn, J. B. J . Chem. Phys. 1975, 63,4114. (35) Guissani, Y.; Guillot, B.; Bratos, S. J . Chem. Phys. 1988.88, 5850. considered in the development of theoretical pair potential for (36) Kochanski, E. J . Am. Chem. SOC.1985, 107, 7869. describing the behavior of the hydronium ion in solution and (37) Fornili, S. L.; Migliore, M.; Palazzo, M. A. Chem. Phys. Lett. 1986, perhaps of other charged species. 125, 419.
Acknowledgment. Calculations were made on a IBM 9021 of the Centre de Calcul de la Universitat de Valencia and on a RISC-6000/530 of the Departament de Quimica Fisica de la Universitat de Valencia (SEUI Project OP90-0042 and a grant of the Conselleria de Educaci6 i Cihcia de la Generalitat Valenciana). I. T. thanks the Departament de Quimica de la Universitat Aut6noma de Barcelona for their warm hospitality and to the Ministeriode Educaci6n y Ciencia (Spain) for a doctoral fellowship. This work was supportedin part by DGICYT (Project PS90-0264). References and Notes Reiss, H.; Heller, A. J . Phys. Chem. 1985, 89, 4207. Gurevich, Y. Y.; Pleskov, Y. V. Electrokhim. 1982, 18, 1477. Lohmann, F. Z . Naturforsch. A 1967, 22, 843. (4) Cited in ref 1. (5) Gomer. R.: Trvson. G.J. Chem. Phvs. 1977.66.4413. (6j Lim, C:; Bashfbrd, 0.;Karplus, M. j . Phys. Chem. 1991,95,5610. (7) Noyes, R. M. J. Am. Chem. SOC.1962, 84, 513. (8) Goldschmidt, H.; Udby, 0. Z . Phys. Chem. 1907, 60, 728. (9) Meot-Ner, M.; Speller, C. V. J . Phys. Chem. 1986,90,6616. MeotNer, M. J . Am. Chem. SOC.1986, 108, 6189 and references therein. (10) (a) Kollman, P. A.; Bender, C. F. Chem. Phys. Lett. 1973, 21,271. (b) Diercksen, G. H. F.; Kraemer, W. P.; Roos, B. 0. Theoret. Chim. Acta 1975, 36, 249.
(1 1) (a) Lau, Y. K.; Ikuta, S.; Kebarle, P. J . Am. Chem. SOC.1982,104, 1462. (b) Lee, N.; Keesee, R. G.; Castleman, A. W. J . Colloid. InterfaceSci. 1980, 75, 555. (1 2) Bockris, J. OM.; Reddy, A. K. N. Modern Electrochemistry; Plenum Press: New York, 1978. (13) (a) Lee, J.; Robinson, G.W.; Webb, S. P.; Philips, L. A.; Clark, J. H.J. Am. Chem. SOC.1986,108,6538. (b) Robinson, G. W. J . Phys. Chem. 1991, 95, 10386. (14) Eigen, M.; DeMaeyer, L.In TheStructureofEIectrolyticSolutions; Hamer, W. J., Ed.; Wiley: New York, 1959; Chapter 5. (15) Triolo, R.; Narten, A. H. J . Chem. Phys. 1975, 63, 3624. (16) Almlof, J. Chem. Scr. 1973, 3, 73. (17) Newton, M. D.; Ehrenson, S. J . Am. Chem. Soc. 1971, 93, 4971. (18) Alagona, G.;Cimiraglia,R.; Lamanna, U. Theoret. Chim. Acta 1973, 29, 93.
(38) Kochanski. E. Chem. Phvs. Lett. 1987. 133. 143. (39) Dang, L. X.; Rice, J. E.; kaldwell, J.; Kollman, P. E. J. Am. Chem. SOC.1991, 113, 2481. (40) Clark, T.; Chandrasekhar, J.; Spitznagel, . - G.W.; Schleyer, P. v. R. J . Comput. Chem. 1983,4, 294. (41) Schlegel, H. B. J. Comput. Chem. 1982,3,214. (42) Frisch, M. J.; Head-Gordon, M.; Schlegel, H. B.; Raghavachari, K.;
Binkley, J. S.;GonzBlez, C.; DeFrees, D. J.; Fox, D. J.; Whiteside, A,; Seeger, R.; Melius, C. F.; Baker, J.; Martin, R. L.;Kahn, L. R.; Stewart, J. J. P.; Fluder, E. M.;Topiol, S.; Pople, J. A. Gaussiun8b; Gaussian, Inc.: Pittsburgh, PA. (43) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (44) McQuarrie, D. Statistical Mechanics; Harper & Row: New York, 1976. (45) Hehere, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab initio Molecular Orbital Theory; Wiley: New York, 1986. (46) Stull, D. R.; Prophet, H. JANAF Thermochem. Tables 1971. (47) Rivail, J. L.; Rinaldi, D. Chem. Phys. 1976, 18, 233. (48) Rinaldi, R.; Ruiz-Lbpez, M. F.; Rivail, J. L. J . Chem. Phys. 1983, 78, 834. (49) Rivail, J. L.;Terryn, B.; Rinaldi,D.; Ruiz-Lbpez, M. F. THEOCHEM 1985, 120, 387. (50) Rinaldi, D.; Rivail, J. L.; Rguini, N. J . Comput. Chem. 1992, 13, 675. (51) Bertrin, J.; Ruiz-Lbpez, M. F.; Rinaldi, D.; Rivail, J. L. Theoret. Chim. Acta 1992,84, 181. (52) Pappalardo, R. R.; Sinchez Marcos, E. J . Chem. Soc.,Faraday Trans. 1 1991.87. 1719. -(53) Sinchez Marcos, E.; Pappalardo, R. R.; Rinaldi, R. J . Phys. Chem. 1991, 95, 8928. (54) Floris, F.; Tomasi. J. J . Comout. Chem. 1990. 10. 616. (55) Floris, F.; Tomasi, J.; Pascuai-Ahuir, J. L. J . Comput. Chem. 1991, 12, 784. (56) Pierotti, R. A. Chem. Reu. 1976, 76, 715. (57) Miertus, S.; Scrocco E.; Tomasi, J. Chem. Phys. 1981, 55, 117. (58) Bonaccorsi, R.; Cimiraglia, R.; Tomasi, J. J . Comput. Chem. 1983, 4, 567. (59) Pascual-Ahuir, J. L.; Silla, E.; Tomasi, J.; Bonaccorsi, R. J . Comput. Chem. 1987,8, 778. (60) (a) Pascual-Ahuir, J. L.;Silla, E. J . Comput. Chem. 1990,lI , 1047. (b) Silla, E.; Tufibn, I.; Pascual-Ahuir, J. L. J. Comput. Chem. 1991, 12, 1077. (61) Mulliken, R. S. J. Chem. Phys. 1955.23, 1833. (62) Aguilar, M. A.; Olivaresdel Valle, F. J. Chem. Phys. 1989,129,439. (63) Tuiibn, I.; Silla, E.; Pascual-Ahuir, J. L. J . Am. Chem. Soc., in press.
