Ab Initio Theoretical Study of Temperature and Density Dependence

Apr 1, 2006 - The temperature and density dependence of the molecular and thermodynamic properties of water is investigated theoretically by means of ...
0 downloads 0 Views 270KB Size
Download PDF
J. Phys. Chem. B 2006, 110, 8451-8458

8451

Ab Initio Theoretical Study of Temperature and Density Dependence of Molecular and Thermodynamic Properties of Water in the Entire Fluid Region: Autoionization Processes Norio Yoshida,† Ryosuke Ishizuka,‡ Hirofumi Sato,§ and Fumio Hirata*,†,‡ Department of Theoretical Study, Institute for Molecular Science, Okazaki 444-8585, Japan, Department of Fundamental Molecular Science, The Graduate UniVersity for AdVanced Studies, Okazaki 444-8585, Japan, and Department of Molecular Engineering, Kyoto UniVersity, Kyoto 615-8510, Japan ReceiVed: NoVember 28, 2005; In Final Form: February 27, 2006

The temperature and density dependence of the molecular and thermodynamic properties of water is investigated theoretically by means of the ab initio electronic structure theory combined with the reference interaction site model method, so-called RISM-SCF. We consider the autoionization process (H2O + H2O h H3O+ + OH-) by regarding H2O, H3O+, and OH- as “solute” molecules in an aqueous solution and evaluate molecular geometry, electronic structure, solvation structure, and the ionic product of water (pKw) of these species as functions of thermodynamic conditions. In our previous paper, we calculated these properties by using essentially the same method in a wide range of density values (0.6-1.4 g/cm3). However, the calculation was limited at rather higher density (>0.6 g/cm3) due to the difficulty of convergence, which is inherent to the hypernetted-chain (HNC) closure. The problem is overcome in this study by employing the KovalenkoHirata (KH) closure which hybridizes the HNC and the mean-spherical approximation (MSA). Here, we present the results for the thermodynamic range of densities from 0.025 to 1.0 g/cm3 and for temperatures from 300 to 800 K including the supercritical point.

1. Introduction Thermodynamic properties, electronic structure, and solvation structure of liquid water show intriguing dependence on temperature and density. The phenomena have not only been a challenge for theoretical chemistry but also are of great interest for the industrial community, since water acts as a solvent for a variety of chemical processes in industrial applications. Especially, the supercritical state of water has drawn a lot of attention for its capability of dissolving organic compounds and accelerating chemical reactions. Despite the recent progress in theoretical chemistry including molecular simulations and quantum chemistry, theoretical study of water in the wide range of its phase diagram remains a challenge. The reason that the problem is so challenging is that the electronic structure of a molecule and many-body properties of liquid are highly coupled in the phenomena which require a treatment based on both the quantum chemistry and statistical mechanics. The typical and the most important example is autoionization, a dissociation reaction, of a water molecule represented as

H2O + H2O h OH- + H3O+

(1)

The ionic product, Kw ) [OH-][H3O+], and its logarithm, pKw ) -log Kw, are not only measures of equilibrium for the process (1) but also standard quantities for protonation (or deprotonation) processes of other chemical species. The ionic product shows strong dependence on temperature and density due to the interplay between the intra- and intermolecular processes. * Corresponding author. † Institute for Molecular Science. ‡ The Graduate University for Advanced Studies. § Kyoto University.

Several theoretical attempts have been proposed to determine pKw. The Car-Parrinello density functional method was employed by Trout and Parrinello to investigate the dissociation process of H2O in liquid water.1,2 In the subsequent work by Sprik, pKw was estimated and showed good agreement with experimental values.3 However, applications of the CarParrinello method are largely limited to relatively small systems due to its computational demand, and it may not be appropriate to study such problems as pKw which require evaluation of free energy in the thermodynamic limit of a Coulombic system. Recently, a combined quantum mechanical and molecular mechanical (QM/MM) method was adopted to compute the temperature dependence of pKw by Yagasaki et al.4 The temperature dependence of pKw showed good agreement with experimental value relatively, though different pressure was assumed. Since such calculations require numerical simulations, their computational cost is rather high. One of the requirements from the fields of experimental work and/or engineering is to evaluate pKw in a wide range of temperature and density (or pressure) values. In the previous study, we investigated the temperature and density dependence of the electronic structure, molecular geometry, solvation structure, and the ionic product of water.5-8 A single molecule in water was designated as “solute”, and other molecules were regarded as “solvent”; a combined ab initio electronic structure theory and the site-site integral equation method of liquids, referred to as RISM-SCF developed by Tenno et al. and Sato et al.,9-11 was applied to the solute-solvent system. By using this method, both the electronic and solvation structures of a solute water molecule were determined in a selfconsistent manner. Although we have calculated such properties in a wide range of density (0.6-1.4 g/cm3) and temperature (300-900 K) values, the low density region (less than 0.6 g/cm3) has not been covered due to numerical divergence. This adverse

10.1021/jp0568834 CCC: $33.50 © 2006 American Chemical Society Published on Web 04/01/2006

8452 J. Phys. Chem. B, Vol. 110, No. 16, 2006

Yoshida et al. standard ab initio molecular orbital (MO) calculations and concerns electronic structures of the isolated species in the gas phase, vac ∆Gvac ) ∆Evac reaction + ∆Gkin

Figure 1. Scheme of evaluation of the free energy change in aqueous solution.

behavior originates from the hypernetted-chain (HNC) approximation which was employed to close the RISM equation. A new closure relation which resolves the problem has been proposed recently by Kovalenko and Hirata.12 It hybridizes the HNC and the mean-spherical approximation (MSA) closures. By the improved closure relation, the liquid-vapor coexisting curves for a variety of liquids have been reproduced successfully.13,14 These results encouraged us to apply the same closure relation to study the combined quantum and many-body properties of water at the low density region by means of the RISM-SCF. In the present study, we apply the KH closure to investigate the properties of fluid water in the wide range of density (0.025-1.0 g/cm3) and temperature (300-800 K) values. In the calculation, a water molecule is considered to be a quantum mechanical solute, which is immersed into solvent consisting of classical water molecules. The dissociated state of the water molecule, H3O+ and OH-, is treated as a quantum solute, too. The organization of this paper is as follows. In section 2, we describe the procedure to calculate the properties of water by means of the RISM-SCF and provide computational details. The results of the calculations and their discussion are presented in section 3. The main stress is put on the phase behavior of the autoionization of water. We also discuss computational results for the molecular geometry and liquid structure, energetics, and electronic structure. Concluding remarks are given in section 4.

Since the method of calculation in the present study is very similar to the previous one, we just provide a brief outline of the procedure.7,8 2.1. Theoretical Method. The thermodynamic cycle of the free energy changes in aqueous solution (∆Gaq) associated with the autoionization reaction is illustrated in Figure 1. We regard the species concerning the reaction, i.e., H3O+, H2O, and OH-, as quantum solute molecules in a classical solvent at infinite dilution. Hereafter, we use symbols ∆ for changes of quantities associated with the chemical reaction and δ for changes due to solvation. Note that both quantities depend on temperature and density. Thus, ∆Gaq can be written in terms of the free energy changes associated with the reaction in vacuo (∆Gvac) and those due to the solvation of reacting species, i.e., δG(H2O), δG(H3O+), δG(OH-):

∆Gaq ) ∆Gvac + δG(H3O+) + δG(OH-) - 2δG(H2O) (2) The free energy change can be also rewritten as

∆Gaq ) ∆Gvac + ∆δEreorg + ∆δµ + ∆δGkin

vac where ∆Evac reaction and ∆Gkin are, respectively, electronic and kinetic energy changes of the reaction in the gas phase. The group ∆δµ corresponds to the difference in the solvation free energies between the reactant and product systems,

∆δµ ) δµ(H3O+) + δµ(OH-) - 2δµ(H2O)

(3)

where ∆Gvac is the free energy difference obtained by the

(5)

In the present study, we employ the KH closure; therefore, the solvation free energy δµ is expressed as

F δµ ) 4π β Rγ

[

1

∑∫dr r2 2(hRγ(r))2Θ(-hRγ(r)) - cRγ(r)

]

1 - hRγ(r)cRγ(r) (6) 2 where hRγ(r) and cRγ(r) are the total and direct correlation functions between solute site R and solvent site γ, respectively. The term Θ(x) is the Heaviside step function which puts the term h2 into effect in the regions of density depletion only, and β ) 1/kBT is inverse of the Boltzmann factor times the thermodynamic temperature. The electronic distortion or polarization upon solvation is described by the electronic reorganization energy, δEreorg, and its change along the autoionization process is given by ∆δEreorg. These quantities are calculated in the course of a self-consistent determination of the electronic structure and solvent distribution by means of the RISM-SCF procedure. ∆δGkin is the kinetic contribution (translational, rotational, and vibrational). Since our primary concern in the present study is the temperature and density dependence of pKw, we concentrate our attention on the quantities relative to that at temperature T ) 273.15 K and density F ) 1.0 g/cm3, that is,

∆pKw ) pKw(F, T) - pKw(F ) 1.0, T ) 273.15)

2. Method of Calculation

(4)

(7)

where pKw is related to ∆Gaq as

∆Gaq ) 2.303RTpKw

(8)

with R being the gas constant. 2.2. Computational Details. We implemented the RISMSCF method on the general atomic and molecular electronic structure system (GAMESS) code.15 For the electronic structure calculations, we employed a triple-ζ basis set augmented by p polarization functions on hydrogen atoms and d polarization and diffusive p (Rp ) 0.059) functions on oxygen atoms are employed in the restricted Hartree-Fock (RHF) method.16 Molecular geometries of all the related species (H2O, H3O+, and OH-) are optimized in both the gas and solution phases under a constraint of the point group symmetry (C2V for H2O and C3V for H3O+). For the purpose of fitting the effective point charge of each atomic site in H2O, H3O+, and OH-, respectively, we generate 628 (for H2O), 725 (for H3O+), and 535 (for OH-) grid points at which the electrostatic poteintial is calculated. To generate fitting points, we employed the geodesic point selection scheme included in GAMESS.17 To evaluate the electronic energy of solute species, we also perform MP2 calculation following the geometry optimization in the RHF level. To calculate the kinetic energy of solute species, Gkin,

Autoionization Processes

J. Phys. Chem. B, Vol. 110, No. 16, 2006 8453

Figure 2. Density dependence of the radial distribution functions between solute species and solvent water at 600 K: (a) O(solute)-O(solvent), (b) O(solute)-H(solvent), (c) H(solute)-O(solvent), (d) H(solute)-H(solvent). These plot of each species (H2O, H3O+, OH-) are shifted.

TABLE 1: Geometrical and Potential Parameters for Solvent Water

this paper, we employ the KH12 closure which hybridizes the HNC and MSA closures and is defined as

geometry r(OH) (Å) ∠HOH (deg) potential parameters σ (Å)  (kcal/mol) partical charge (e)

O

H

3.1660 0.1554 -0.82

1.0000 0.0540 0.41

second derivatives of the free energy of solute species are numerically evaluated using the analytical expression for the free energy gradient (details of the evaluation of an analytical derivative of free energy have been described in the Appendix.) The solvent-solvent and solute-solvent RISM calculations were performed on a grid of 2048 points with the grid spacing of 0.05 Å. The solvent model we use in the present study is summarized in Table 1. This model has been successfully used in the liquid phase calculations with the Lorentz-Berthlot combination rule. The same Lennard-Jones parameters are used for solute species. All calculations were carried out in a range of density of 0.025-1.0 g/cm3 and temperature of 300-800 K. 2.3. Effects of Closure Relation. Calculations in our previous papers6-8 were performed at rather higher density (>0.6 g/cm3), due to the difficulty of convergence of the solution to the RISMSCF equation with the HNC closure at low density. In

{

exp(dRγ(r)) for dRγ(r) e 0 1 + dRγ(r) for dRγ(r) > 0

(9)

dRγ(r) ) -βuRγ(r) + hRγ(r) - cRγ(r)

(10)

gRγ(r) )

1.0000 109.27

It combines the HNC approximation in the spatial regions of density depletion where gRγ(r) e 1 and the MSA in the enrichment regions of gRγ(r) > 1. Despite the heuristic nature of this closure, it turns out to be quite successful in explaining a variety of solution processes including vapor-liquid phase transition.14 To verify the effect of the closure relation, several results of calculations using various closures are shown in Table 2. The difference between the results obtained with the KH and HNC closures in this paper is quite small, while in the previsous paper it is a bit large, which originates from the dissimilarity of computational conditions, i.e., a method to generate the fitting points for the effective point charge, a grid of correlation functions, and details of a basis set of solute molecular wave functions. 3. Results and Discussions 3.1. Solvation Structures. The density dependence of the solute-solvent radial distribution functions (RDFs) of 600 K is shown in Figure 2. The reason that we have chosen this

8454 J. Phys. Chem. B, Vol. 110, No. 16, 2006

Yoshida et al.

TABLE 2: Effect of the Closure Relationa 293.15 K KH r(OH) of H2O (Å) ∠HOH of H2O (deg) dipole moment of H2O (D) ∆pKw c

HNC

573.15 K HNCb

KH

HNC

HNCb

0.949 0.949 0.956c 0.948 0.948 0.949c 105.8 105.8 106.1c 105.9 105.9 106.4c 2.81 2.83 2.91d 2.70 2.70 2.67d -0.25 -0.41 -0.71c -8.62 -5.13 -7.17c

a These values are at 1.0 g/cm3. b Result of previous calculation. Reference 7. d Reference 6.

temperature is that the critical point for water within the framework of RISM/KH theory is located at T ) 600 K and F ) 0.12 g/cm3.13 RDFs between the solute and solvent oxygen atoms are shown in Figure 2a. The density dependence of the height as well as the position of the first peak can be explained in terms of an interplay of these two forces acting between the solute and solvent atoms: the electrostatic interaction between the partial charges on these two atoms, e.g., a hydrogen bond, and the core repulsion due to the steric effect. In the case of O-O pairs, interactions are purely repulsive since both atoms have negative charges. Therefore, the position and the height of the first peak are determined essentially by the packing effect. RDFs of OH- and H2O show similar behavior in a sense that the first peak becomes broader and lower with decreasing density, which is the reflection of weakening of the packing effect. RDFs of H3O+ show the opposite tendency in terms of the peak height, indicating that some associative configuration between a pair of the species is favored in the low density region. The configuration is stabilized due to the increased positive charge of the solute hydrogen atom, which enhances the hydrogen bond with solvent water molecules. Both the first minimums and the second peaks in RDFs for all the solute species disappear at a density of 0.025 g/cm3. These features are typical of those in the vapor phase. Solvation structures of solvent hydrogen atoms around the solute oxygen atom are represented in Figure 2b. The extremely well-defined peak around R ) 1.8 Å in the RDFs of OH- and H2O represents the hydrogen bond between solute oxygen and solvent hydrogen. The intensity of the hydrogen-bond peak becomes stronger with decreasing density. At the same time, the first peak in the RDF of H3O+ shifts monotonically in the direction of a larger separation of sites. These tendencies are again due to the interplay between the packing effect and the electrostatic interactions, which change their strength with decreasing density. The hydrogen-bond peak entirely disappears in the case of H3O+, because the positive charge of H3O+ repels the positive partial charges of water hydrogens. The dependence of the H(solute)-O(solvent) RDFs upon the density, Figure 2c, can be explained in a manner similar to the case of O(solute)H(solvent), i.e., in terms of the interplay between two forces. Since the positive charges on H(solute) of H3O+ are enhanced, the hydrogen bond peak at 1.8 Å becomes higher than that in the pure water. It is also conspicuous that the associative configuration becomes more stable with decreasing density. On the other hand, the hydrogen-bond peak weakens in the case of OH-, because the negative charge of OH- repels the negative partial charges of water oxygen by the electrostatic force. RDFs of H2O are essentially the same as in the case of O(solute)H(solvent), and the difference is mostly due to polarization of solute species. (If the solute species have the same partial charges and geometrical properties as the solvent, RDFs of O-H and H-O of H2O shall be identical.) Finally, RDFs of H(solute)-H(solvent) are shown in Figure 2d. Although electrostatic interaction between the hydrogen sites is purely

Figure 3. Density and temperature dependence of the dipole moment of solute species: (a) H2O, (b) H3O+, (c) OH-. The units are Debye.

repulsive, the first peaks in RDFs of H-H are correlated with first peaks in RDFs of H-O, or the hydrogen bond, and increase with decreasing density in the case of H3O+ and H2O. The H-H RDFs around OH- shows a distinct behavior in density dependence: both the peak height and width enhance with increasing density. The increase in the peak height is due to strengthening of the O-H hydrogen bond as has been already seen in Figure 2b. The broadening of the peak indicates that the water molecules hydrogen bonded to OH- become more free to rotate around the solute, so that the H-H distance has broader distribution. 3.2. Solute Structures. The temperature and density dependence of the electronic structure of the solute species, represented by the dipole moment, is shown in Figure 3a-c for H2O, H3O+, and OH-, respectively. Dipole moments were estimated from the solute electron density distribution which was determined in a self-consistent manner. When polar solvents are concerned, enhancement of the dipole moment of a solute molecule contributes to stabilization and the electronic structure is polarized to increase the dipole moment. For all the species, the dipole moment is largely enhanced by solvation compared to that in the gas phase, as indicated in Table 3. Both OH- and

Autoionization Processes

J. Phys. Chem. B, Vol. 110, No. 16, 2006 8455

Figure 4. Density and temperature dependence of the intramolecular geometry of solute species: (a) O-H bond length of solute H2O, (b) H-O-H angle of solute H2O, (c) O-H bond length of solute H3O+, (d) H-O-H angle of solute H3O+, (e) O-H bond length of solute OH-. The units of length and angle are angstroms and degrees, respectively.

TABLE 3: Geometrical and Electronic Properties for H2O, H3O+, and OH- in the Gas Phase r(OH) (Å) ∠HOH (deg) dipole moment (D) partial charge on O (e) partial charge on H (e)

H2O

H3O+

0.943 106.927 2.199 -0.835 0.417

0.962 114.297 1.530 -0.624 0.541

TABLE 4: Effects of MP2 Correction at 300 K and 1.0 g/cm3 a

OHδEMP2(H2O) δEMP2(H3O+) δEMP2(OH-) ∆EMP2 vac ∆Gaq pKw

0.948 1.632 -1.307 0.307

H3O+ are polarized more than H2O. OH- is polarized especially strongly. This is because the molecular shape of the species is linear and the effect of induction of the dipole is effective. In general, polarization of the electronic structure is reduced with rising temperature or decreasing density and approaches the electronic structure in an isolated molecule. Three factors in increasing the dipole moment are considered. The most important factor for the dipole moment is the induced partial charge on each atom due to electrostatic field of the solvent. The other factors are bond length between oxygen and hydrogen atom, r(OH), and the bond angle of HOH, ∠HOH: increasing the bond length and/or decreasing the bond angle results in increasing the dipole moment. The temperature and density dependence of these intramolecular geometrical properties is shown in Figure 4. With rising temperature or decreasing density, both r(OH) and ∠HOH of H2O and H3O+ change monotonically. The bond angle seems to be determined by the subtle balance between two effects due to polarization. One is an intramolecular repulsion between hydrogen atoms. It is enhanced by the electronic polarization upon solvation and makes the angle wider. On the other hand, the widening of the bond angle gives rise to lesser stabilization since it decreases the dipole moment of the molecule. Changes of r(OH) and ∠HOH enhance the dipole moment of solute species with increasing “solvation effects”.

a

RHF

RHF-MP2

41.20 30.02

0.04 -2.01 -1.21 -7.87 30.02 21.88

The units of the energy components are kilocalories per mole.

The temperature and density dependence of r(OH) of OHshows different behavior from those of H2O and H3O+. The bond length of OH- is shortened with decreasing density and increasing temperature until the density becomes ∼0.2 g/cm3. It takes a broad minimum around F ) 0.1-0.2 g/cm3 and then elongates again in the lower density region. The reason for this phenomenon can be explained in terms of the electrostatic stabilization and the cavity formation. The shorter bond length is favorable at higher density in terms of the cavity formation, because greater work is required to make a cavity with larger size. On the other hand, the longer bond length is more stable at higher density from the point of view of the electrostatic interaction. As the density decreases, the two tendencies balance to make a minimum at some density. 3.3. Free Energy and pKw. Theoretical results of pKw at T ) 300 K and F ) 1.0 g/cm3 are shown in Table 4. The theoretical value without the MP2 correction is significantly overestimated compared with the experimental value, 14. It should be noted that the MP2 correction makes drastic improvement. It is well-known that the electron correlation of an anion is large. The electron correlation of OH- becomes dominant in the MP2 component of energy change of the autoionization reaction in vacuo, ∆EMP2 vac . This effect drives the autoionization

8456 J. Phys. Chem. B, Vol. 110, No. 16, 2006

Figure 5. Density and temperature dependence of ∆pKw, the difference between the pKw(F, T) and pKw(1.0 g/cm 3, 273.15 K), where F and T are the solvent water density and the temperature of the system. The theoretical result of this work and an experimental value19 are shown in parts a and b, respectively.

Yoshida et al. reaction, eq 1, toward the right-hand side. The MP2 correction for the change in solvation energy is also effective for the ionic species, unlike the case of H2O. Despite the fact that rather large disagreement in magnitude between the experimental and theoretical results for pKw still remains, their qualitative behaviors seem to be in harmony. Therefore, we concentrate our attention on quantities relative to those at a temperature of T ) 273.15 K and a density of F ) 1.0 g/cm3 as defined in eq 7. In Figure 5, the density and temperature dependence of ∆pKw is shown. The ∆pKw increases monotonically with decreasing temperature and density. Although the behavior of ∆pKw qualitatively shows agreement with experimental results,18,19 the extent of change in ∆pKw with changing temperature and/or density is overestimated. The computer simulation result of the temperature dependence of pKw at 20 MPa by Yagasaki et al.4 shows good agreement with the experimental result at 32 MPa. Considering that their calculation treats a cluster including several water molecules as a solute by QM, the disagreement of our results with the experiments is mainly due to the lack of the quantum interactions, i.e., charge transfer, between solute species and solvent water in our calculations. Since pKw is related to ∆Gaq by eq 8, we consider components of ∆Gaq to investigate the trend of pKw. The reason why the temperature dependences of ∆pKw and ∆Gaq become opposite is that pKw is inversely proportional to T. The ∆Gaq is decomposed into several components appearing in the left-hand side of eq 3. The results of calculations are shown in Figure 6. The behavior of ∆Gvac shows a small temperature dependence in Figure 6c, since it only depends on an intramolecular vibrational energy of solute species. The terms ∆δGkin and ∆δEreorg also show week dependence on density and temperature. Therefore, the most

Figure 6. Density and temperature dependence of the Gibbs free energy change for the autoionization process of water: (a) ∆Gaq, (d) ∆δµ. Parts b and c are energy components of ∆Gaq plotted against density and temperature, respectively: (filled squares) ∆Gaq; (filled circles) ∆Gvac; (open circles) ∆δGkin; (filled triangle) ∆δEreorg; (open squares) ∆δµ. The units of the y-axis are kilocalories per mole.

Autoionization Processes

J. Phys. Chem. B, Vol. 110, No. 16, 2006 8457 4. Summary

Figure 7. Density and temperature dependence of the excess chemical potential: (a) H2O, (b) H3O+, (c) OH-. The units of the y-axis are are kilocalories per mole.

important contribution to the temperature dependence of ∆Gaq comes from ∆δµ, as shown in Figure 6d. The value of ∆δµ increases with decreasing density and with increasing temperature. Components of ∆δµ in the left-hand side of eq 5 are shown in Figure 7. Although ∆δµ as a sum of δµ(H3O+), δµ(OH-), and -2δµ(H2O) shows a monotonic behavior, the minimum points appear in its components, δµ(H2O) and δµ(OH-) (see Figure 7a and b). The behavior indicates the existence of two competing effects which interplay. The excess chemical potential of solvation is affected by two major effects, which are electrostatic stabilization and cavity formation. The components δµ(H2O) and δµ(OH-) exhibit a similar monotonic decrease with increasing density in the low density region, which can be attributed essentially to the density dependence of the free energy of electrostatic stabilization. On the other hand, they increase with increasing density in the high density region, which can be attributed to effects of cavity formation. The term δµ(OH-) decreases monotonically with increasing density, which is due to strong electrostatic stabilization effects, as we have seen in Figure 2.

The temperature and density dependence of the molecular and thermodynamic properties of water was investigated theoretically by means of the RISM-SCF theory. We considered the autoionization process (H2O + H2O h H3O+ + OH-) by regarding H2O, H3O+, and OH- as solute molecules in an aqueous solution and evaluate the molecular geometry, electronic structure, solvation structure, and ionic product of water (pKw) of these species as functions of thermodynamical conditions. These calculations were performed in a wide range of density (0.025-1.0 g/cm3) and temperature (300-800 K) values by adopting the KH closure. Solvation structures of the three species, represented by the site-site radial distribution functions (RDFs) between the solute and solvent, have also exhibited a marked dependence on the temperature and density. The density dependence of an RDF at short separations can be explained in terms of an interplay between essentially two different “solvation” structures, the contributions of which change as the density decreases. The first of those is that characteristic to a simple fluid in which the packing effect dominates the configuration, which is favored at a higher density. The other structure is due to hydrogen bonding between the solute and solvent species, the contribution of which is greater at a lower density. The temperature and density dependence of the solute electronic structure and geometrical properties was also evaluated. In all the species, the dipole moments were largely enhanced by solvation compared to those in the gas phase. The geometrical properties of solute species were changed by solvation so that the dipole moment increases. The bond length r(OH) was elongated with increasing density and decreasing temperature, and ∠HOH became narrow with increasing density and decreasing temperature. On the other hand, the temperature and density dependence of r(OH) of OH- shows a different behavior from H2O and H3O+. The reason for this phenomenon can be explained in terms of electrostatic stabilization and cavity formation. The results for ∆pKw obtained from the theory have shown a monotonic increase with increasing density at all the temperatures investigated in good accord with the experimental observation. The behavior is determined essentially by the difference in the excess chemical potentials. An excess chemical potential of solvation is affected by two major effects, which are electrostatic stabilization and cavity formation. Although the MP2 correction improves the theoretical results of pKw, a nontrivial error still remains. This error may have been caused by the assumptions that we adopted: only two water molecules are concerned with the reaction (1) and quantum interaction or charge transfer between solute species and solvent water is neglected. Therefore, the estimation of pKw is likely to be improved by taking such effects into consideration. Acknowledgment. This work is supported by the Grant-in Aid for Scientific Research on Priority Area of “Water and Biomolecules” of the Ministry of Education in Japan, Culture, Sports, Science and Technology (MONBUKAGAKUSHO). We are also grateful for the support by the grant from the NAREGI Nanoscience Project of the same ministry. Numerical calculations were partly carried out in the Research Center for Computational Science, Institutes of Natural Science. Appendix The Lagrangian of the system for the RISM-SCF/MCSCF theory with the HNC closure11 was already defined in detail

8458 J. Phys. Chem. B, Vol. 110, No. 16, 2006

Yoshida et al.

previously. Therefore, we concisely describe the Lagrangian for the case of the KH closure system. We define the total Helmholtz free energy A as the sum of the electronic energy of the solute molecule Esolute and the excess chemical potential concerning the solute-solvent interaction δµ, given by eq 6

A ) Esolute + δµ

(11)

Esolute can be estimated by ab initio electronic structure methods such as MCSCF. For the chemical potential term, we adopt the free energy derived from the KH closure relation by Kovalenko and Hirata,12

δµ )

[

β

4πF β

∑∫

∑ R′γ′

2 2

]

[

∫dr r2 {e-βu ∑ Rγ

{

Rγ(r)+tRγ(r)

}

]

1 1 - tRγ(r) - hRγ(r)tRγ(r) + hRγ(r)2 2 4π dk k2 -cˆ Rγ(k)Fhˆ Rγ(k) + 3 (2π) β Rγ 1 cˆ Rγ(k) cˆ R′γ′(k) ω ˆ RR′(k)χˆ γγ′(k) 2R′γ′

∑∫

{

∑

}

im(Sim - δim) ∑I CI2 - 1) - ∑i ∑ m

Variations with respect to the functions yield

δCI +

Γijklgkl - ij ∑ij δφi|γijh + ∑ kl

(

∂

-

λ

4πF β

∫dr r2e-βu ∑ Rγ

Rγ(r)+tRγ(r)

)| 〉

φj (14)

References and Notes

(12)

where χγγ′ are the pure solvent site density pair correlation functions, ωRR′ are intramolecular correlation function containing the structural information of the solute, uRγ is the solute-solvent intearaction potential, and a circumflex over the quantity, represents the Fourier transform of the corresponding function. The total Helmholtz free energy A can be regarded as a functional of correlation functions hRγ, cRγ, tRγ, of the MO coefficient, and of the CI coefficient. Imposing constrains of orthonormality of configuration state functions and one particle orbitals, we define the Lagrangian of the system as

L ) A - E(

〈

}

where all the notations are the same as in ref 11. We find that the variational condition for the free energy functional leads to the KH approximation, the RISM-OZ equations, and the modified MCSCF equations involving the solvation effect. The first derivative of the free energy with respect to the nuclear coordinates of the solute molecule is derived from eq 14 employing the same strategy as described in ref 11. Using the variational conditions, eq 14, we can find a quite similar expression of the energy gradient as the previous one.

}Θ(-hRγ(r)) +

(tRγ(r) - uRγ(r))2 Θ(hRγ(r)) -

2

{

∑I ∑J HIJCJ - ECI

∑λ bλ ∂q

γij

tRγ(r) - uRγ(r) + 1 +

1

∫dr r2 ∑ Rγ

[{(e-βuRγ(r)+tRγ(r) - hRγ(r) - 1)Θ(-hRγ(r)) + (tRγ(r) - uRγ(r) - hRγ(r))Θ(hRγ(r))}δtRγ(r) + {-tRγ(r) - hRγ(r) - cRγ(r)}δhRγ(r)] 4π dk k2{Fhˆ Rγ(k) + 3 (2π) β Rγ cˆ R′γ′(k) ω ˆ RR′(k) χˆ γγ′(k)}δcRγ(r) +

∫dr r2 2(hRγ(r))2Θ(-hRγ(r)) ∑ Rγ 1 cRγ(r) - hRγ(r)cRγ(r) ) 2

-

4πF

1

4πF β

δL ) -

(13)

(1) Trout, B. L.; Parrinello, M. Chem. Phys. Lett. 1998, 288, 343. (2) Trout, B. L.; Parrinello, M. J. Phys. Chem. B 1999, 103, 7340. (3) Sprik, M. Chem. Phys. 2000, 258, 139. (4) Yagasaki, T.; Iwahashi, K.; Saito, S.; Ohmine, I. J. Chem. Phys. 2005, 112, 144504. (5) Maw, S.; Sato, H.; Ten-no, S.; Hirata, F. Chem. Phys. Lett. 1997, 276, 20. (6) Sato, H.; Hirata, F. J. Chem. Phys. 1999, 111, 8545. (7) Sato, H.; Hirata, F. J. Phys. Chem. A 1998, 102, 2603. (8) Sato, H.; Hirata, F. J. Phys. Chem. B 1999, 103, 6596. (9) Ten-no, S.; Hirata, F.; Kato, S. Chem. Phys. Lett. 1993, 214, 391. (10) Ten-no, S.; Hirata, F.; Kato, S. J. Chem. Phys. 1994, 100, 7443. (11) Sato, H.; Hirata, F.; Kato, S. J. Chem. Phys. 1996, 105, 1546. (12) Kovalenko, A.; Hirata, F. J. Chem. Phys. 1999, 110, 10095. (13) Kovalenko, A.; Hirata, F. Chem. Phys. Lett. 2001, 349, 496. (14) Hirata, F., Ed. Molecular Theory of SolVation; Kluwer: Dordrecht, 2003. (15) Schmidt, M. W.; Baldrige, K. K.; Boats, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (16) Dunning, T. H., Jr. J. Chem. Phys. 1970, 53, 2823; 1971, 55, 716. (17) Spackman, M. A. J. Comput. Chem. 1996, 17, 1. (18) Mesmer, R. E.; Marshall, W. L.; Palmer, D. A.; Simonson, J. M.; Holmes, H. F. J. Solution Chem. 1988, 17, 699. (19) Palmer, D. A.; Fernandez-Prini, F R.; Harvey, A. H. Aqueous Systems at EleVated Temperatures and Pressures; Elsevier/Academic Press: Oxford, 2004.