J. Phys. Chem. 1996, 100, 11165-11174
11165
A Continuum Solvation Model Including Electrostriction: Application to the Anisole Hydrolysis Reaction in Supercritical Water Huaqiang Luo and Susan C. Tucker* Department of Chemistry, UniVersity of California, DaVis, California 95616 ReceiVed: August 15, 1995; In Final Form: October 25, 1995X
We present a simple and efficient method for calculating the effect of electrostriction on the solvation of molecular species. This new method introduces solvent compressibility into the standard continuum dielectric description of the solvent. It is based on a numerical grid-based algorithm for solving Poisson’s equation and is therefore applicable to arbitrary charge distributions in cavities of arbitrary shape. Application to the anisole hydrolysis reaction in supercritical water shows that the effects of solvent compression can be large; for example, compression is observed to lower the free energy barrier to reaction by as much as 14 kcal/mol. Additionally, compression effects are found to qualitatively alter the pressure dependence of the reactivity. Finally, the inclusion of compression effects is shown to bring the calculated pressure dependence of the activation barrier into significantly better agreement with experiment.
I. Introduction Supercritical fluids (SCF) are of technological interest because they have been shown to promote novel solute reactivity.1 Additionally, the observed reactivity shows a strong dependence on the thermodynamic conditions of the solvent,2-7 giving the promise of controlling reactivity by adjusting the thermodynamic variables. Supercritical water (SCW) is of particular interest, because it poses no waste disposal problems and because a number of applications, including efficient oxidative destruction of hazardous wastes, have recently been developed.8-12 The dominant reaction mechanisms in SCW have been shown to be extremely sensitive to variations in pressure for temperatures near (but above) the critical temperature, Tc. In particular, it has been found that for low pressures homolytic free radical mechanisms dominate, whereas for high pressures more chemically selective heterolytic, ionic mechanisms dominate.3,13 This variation in reaction process with pressure is known to be tied to the large variation in SCW’s density, and thus its dielectric properties, with pressure in this regime.2-4,14-16 For example, at 380 °C, the dielectric constant of water varies from ) 2 at P ) 21 MPa to ) 14 at P ) 47 MPa.17 This extreme variation in density-dependent solvent properties with moderate changes in pressure is typical of SCF’s in the region of phase diagram where their isothermal compressibility is large. Because this region, which we denote as the “compressible regime”, offers the most adjustable control of solvent properties, it is the region in which most studies of reactivity in SCF’s are performed. Note that while this regime includes the so-called “critical region”, it extends over a much large region of the phase diagram above the critical point (T > Tc). In order to take advantage of the promise of solvent control of chemical reactivity in SCF’s, we must be able to predict the temperature and pressure dependence of the solvent effects. Toward this end, we propose a new method for calculating the pressure dependence of solute rate constants in SCF’s.18 We test our method by application to a charge separation SN2 reaction, the hydrolysis of anisole, which has been studied experimentally in SCW.13,19 The primary manner in which solvents affect solute reaction rates is by altering the relative stabilities of the reacting solute X
Abstract published in AdVance ACS Abstracts, May 1, 1996.
S0022-3654(95)02359-8 CCC: $12.00
at various positions along its reaction path, thus changing the activation barrier to reaction. As such equilibrium solvation effects alter the rate constant exponentially, in contrast to nonequilibrium solvent effects which alter the prefactor, we consider here only the equilibrium solvent effects.20,21 The key quantity to evaluate, once the solute’s gas phase reaction profile is known, is the free energy of solvation as a function of the solute reaction path. In general, free energies of solvation may be calculated by either molecular dynamics or Monte Carlo simulation techniques or by dielectric continuum models. While simulation techniques provide the most reliable results, they are computationally intensive and, in addition, require the design of fitted interaction potentials for each solute-solvent system studied.22 It requires a significant effort to determine the critical parameters, Tc and Pc, for each solvent model (which are often quite different from the experimental values).23-25 These parameters must be known precisely in order to run simulations in the compressible regime, due to its vicinity to the critical point. They are also required for making comparisons between calculated results and experiment, as it is the reduced thermodynamic variables which are relevant (e.g., Tr ) T/Tc, Pr ) P/Pc). In addition, it is difficult to create solvent potentials which accurately predict the pressure and temperature dependence of the important solvent properties, such as the compressibility and the dielectric constant, over large variations in density. Finally, simulation techniques become prohibitively expensive near the critical point due to convergence difficulties as the time and length scales for collective solvent motions become very long.26-28 In contrast, continuum solvation calculations provide a simple and efficient approach for evaluating the electrostatic contribution to the solvation free energies.29-31 These simple methods, in which the solute is modeled as a charge distribution in a dielectric continuum, eliminate the need for fitted potential functions. Specifically, the solute charge distribution may be taken directly from quantum chemistry calculations, and the only necessary solvent parameter is the dielectric constant, for which the experimental value may be used directly. Under ambient conditions, such models have proven to be effective for systems ranging from simple ions to complex biological molecules.29,30,32 Recently, dielectric continuum methods have been shown to fail under supercritical conditions.14,15 Both an initial study by © 1996 American Chemical Society
11166 J. Phys. Chem., Vol. 100, No. 26, 1996 Tucker and Gibbons on the anisole hydrolysis reaction in SCW14 and a subsequent study by Bennett, Johnston, and Rossky on the Cl- plus CH3Cl reaction in SCW15 postulate that standard continuum models fail in this regime because they assume that the solventsi.e., the dielectric continuumsis incompressible, while SCF’s of interest are characterized by large isothermal compressibilities. A large isothermal compressibility is the macroscopic manifestation of large density fluctuation in the pure solvent which are self-correlated over a long range. Note that as the critical point is approached form above, the correlation length, and thus the compressibility, diverge. A large solvent compressibility may affect the solvation of a solute in two ways. First, as evidenced by the large density fluctuations in the pure solvent, there is very little free energy cost for solvent compression. Thus, upon compression of the solvent around the solute, the energetic gain associated with attractive solvent-solute interactions will far outweigh the free energy cost to the solvent of this compression. As a result, enhanced solvent density in the vicinity of the soluteswhere “vicinity” is defined by the range of the solvent-solute potential interactionssis expected. We refer to this general phenomenon as clustering, as will be discussed further below. A secondary effect may also arise, as a result of the long correlation length of solvent-solvent density correlations in the compressible solvent. When this correlation length becomes long relative to the range of the solvent-solvent potential interactions, the enhanced solvent density in the vicinity of the solute (which arises as a direct result of solvent-solute potential interactions, as discussed above) may induce solvent density enhancements out to much further distances from the solute. The range of these secondary density enhancements is determined by the correlation length and, as such, are expected to become important only very near to the critical point. Additionally, the density enhancements resulting from this effect would be primarily outside of the range of the solute’s potential and thus are not expected to have much effect on the solute’s free energy of solvation. These effects are therefore neglected herein. The existence of enhanced solvent density around solutes33 in SCF’s is experimentally well established, both from partial molar volume measurements34,35 and from spectroscopic studies.4,36-40 Such density enhancements have been characterized by computer simulation studies and by integral equation theory.41-43 The formation of such solvent density enhancements can significantly effect solute reactivity, as has been demonstrated experimentally. Specifically, the decomposition of R-chlorobenzyl methyl ether was studied via spectroscopic shift measurements,4 and the reduction of I2 to I- was studied electrochemically.44 Both studies used an analysis based on dielectric continuum theory to show that some form of solvent density enhancement (range and mechanism unspecified) plays an important role in determining reactivity. The continuum model studies mentioned above also confirm that free energies of activation profiles calculated without accounting for these solvent density effects are severely in error as compared to experiment14 and computer simulation.15 Here we address the failure of the traditional, incompressible continuum theories to describe solvation in SCF’s by introducing solvent compressibility into the model. Our new model allows for density increases in the solvent which arise directly from the solvent-solute potential interactions, i.e., clustering. However, continuum models neglect short-range, molecular potential interactions and incorporate only the electrostatic interactions. Hence, our new model will include only density enhancements
Luo and Tucker resulting from solvent-solute electrostatic interactions, i.e., electrostriction. Additionally, use of a continuum description of the solvent means that detailed molecular structure for the solvent near the solute, i.e., within one or two solvation shells, is neglected. Similarly, electrostriction effects over this short range are only crudely approximated. There exists one previous continuum solvation model due to Wood, Quint, and Grolier,45 which is based on a thermodynamic exposition of Frank.46 This model allows for electrostriction, but it is restricted to the case of a single point charge in spherical cavity and hence cannot be used to study reactivity.45 These authors showed that the solvent density enhancement around ions in an SCF solvent can be quite large, although they did not evaluate the effect of this compression on the solvation free energies in this case.47 Additionally, this model was used by Flarsheim et al. to describe the clustering effects on I- in their study of the I2/I- reduction and found to be adequate.44 We have extended the ideas of Wood et al.45 to the case of a general charge distribution in a cavity of arbitrary shape by taking advantage of the finite-difference grid-based Poisson-Boltzmann equation used in the incompressible continuum method of Honig and co-workers.48-50 We test the method by application to the anisole hydrolysis reaction in SCW, i.e.
PhOCH3 + OH2 f PhO- + CH3OH2+
(1)
This reaction provides an excellent test case because it does not occur under ambient conditions (in the absence of strong acid), yet at the experimental conditions of 380 °C (Tr ) 1.01) and high pressures this mechanism is dominant. At low pressures this mechanism is subsumed by competing pyrolysis pathways. It follows that the rate constants for these mechanisms are strongly pressure dependent. The rate constant for the hydrolysis reaction has been measured experimentally and found to increase from 4.6 × 106 to 62 × 106 s-1 with a pressure increase of only ∼25 MPa.13 In addition, this reaction has been determined experimentally to be a charge separation SN2 reaction.19 Such an asymmetric ionic process in which the charge distribution changes dramatically along the reaction path will provide a stringent test of the calculation of the electrostatic solvation effects. The new method is described in section II. Relevant details and previous results for the anisole hydrolysis reaction are summarized in section III. Results of the application of the new method to this reaction are presented and discussed in section IV. Conclusions follow. II. Compressible Continuum Method A. Calculation of the Electrostatic Potential. To develop a continuum solvation model which incorporates electrostriction, we work within the framework of the finite-difference gridbased method used by Honig and co-workers to solve the standard incompressible dielectric continuum solvation problem.48,49 In this standard model, the solute is described by a charge distribution, Fe(r), which resides in a cavity of fixed dielectric constant, cav. The solute is embedded in the solvent which is represented by a fixed dielectric constant, s. The finite-difference grid algorithm is then used to iteratively find a numerical solution to Poisson’s equation,48
∇‚[(r) ∇φ(r)] ) -4πFe(r)
(2)
where φ(r) is the electrostatic potential. The position dependence of the dielectric constant, (r), is simply a three-
A Continuum Solvation Model Including Electrostriction
J. Phys. Chem., Vol. 100, No. 26, 1996 11167
dimensional step function at the cavity boundary; i.e., it has the constant value s outside the cavity and the constant value cav inside the cavity. The dielectric constant, s, characterizes the solvent’s polarization response to the electric field created by the solute. In the continuum model framework, the properties of the solvent appear only through the bulk dielectric constant, or permittivity, s. Thus, additional solvent response may be incorporated by allowing the permittivity to depend on the solute-induced field. Here, we generalize the model to allow for electrostriction, in which the local solute field induces a local solvent compression, or density enhancement, which in turn causes an increase in the local dielectric constant due to the increased dipole density. In this case, the local permittivity becomes a function of the (field-dependent) local density, i.e.
s(r;T) ) s(F[E(r)];T)
(3)
This function can be built up from the heirarchy of functions s(F,T), F(E), and E(r). The first two functions, s(F,T) and F(E), are properties of the pure solvent. They are therefore independent of the solution to Poisson’s equation, and they may be determined from experimental data. In contrast, the third function, E(r), must be determined through the self-consistent solution of Poisson’s equation (eq 2) and eq 3. In the present work, the density dependence of the dielectric constant for SCW, s(F), is evaluated at fixed temperature from the empirical function of Archer and Wang.51 A general expression for the dependence of the density on the local electric field, F(E), at constant temperature and chemical potential, was derived by Frank,46 i.e.
( )( )
∂F 1 dF ) 0F 2 ∂P
E,T
∂s ∂F
E,T
dE2
(4)
in SI units, with 0 the vacuum permittivity. Note that the change in density with electric field is directly proportional to the isothermal compressibility of the solvent, κT ) -(1/V)(∂V/ ∂P)T, and is thus to be important for highly compressible SCF’s. The density dependence of the necessary derivatives in eq 4 is evaluated numerically from the equation of state17 and s(F).51 Both derivatives are approximated by their zero E-field values, as discussed in the Appendix. We note that although the E dependence of the isothermal compressibility can be derived from thermodynamic relations,52 its effect on the present calculations is found to be negligible. The local field-dependent density function, F(E) - F(E)0), is found by numerical integration of eq 4 and tabulated for 4000 values of the upper limit, E2 (see Figure 1). To combine the functions s(F,T), F(E), and E(r) to create the position-dependent dielectric function, eq 3, requires the selfconsistent solution of Poisson’s equation, eq 2, with eq 3. This nonlocal self-consistency is easily handled by a finite-difference grid algorithm. This type of algorithm has been used in the standard incompressible model to handle the nonlocal selfconsistency of the electrostatic potential, φ, required by Poisson’s equation.48 In the incompressible continuum calculation, the solventsolute system is placed on a cubic grid of mesh size h. Poisson’s equation is then solved by making an initial guess for the potential, φi, at each of the grid points. An updated estimate of the local value of the potential at each grid point is then found according to48
Figure 1. Equilibrium density in the presence of an E-field for SCW with E ) 0 approximation in eq 4, solid lines, without (see Appendix), dashed lines; states 1-3.
φinew )
∑j jφjold + 4πqi/h ∑j j
(5)
where qi is the free charge (from Fe(r)) at grid point i, and the summations run over its six neighboring grid points, {j}. The local potential, φold j , is the previous value of φ at the grid point j, while j is the value of the dielectric constant (s or cav) at the midpoint of the line connecting grid point i to j. To include electrostriction into this model, we use eq 3 to update the value of the local dielectric constant at each grid point for each new set of solutions {φi}. Thus, in the initial iteration of eq 5, each j will be either the expected bulk value of the solvent dielectric at zero field, s0, or the value of the dielectric in the solute cavity, cav. On subsequent iterations, each j in the solvent region will be updated. Specifically, the updated local dielectric at the ij midpoint is found as j ) old s(F[Eold j ],T), where the local field at the ij midpoint, Ej , is evaluated from the difference of the (old) potential values at the neighboring grid points. Furthermore, s(F[Eold j ],T) is found from the functions F(E) and then s(F,T) described above. The functions {φi} and {j} are then updated alternately until the electrostatic potential achieves a self-consistent solution. In the present calculations a grid size of 65 × 65 × 65 points is used. The initial calculation is performed with a large grid spacing so that the box is about 5 times larger than the solute molecule, and Coulombic boundary conditions (φ(r) ) ∑iqi/4π0s0(r - ri) for r at the grid boundary) are used. Focusing49 is then used to repeat the calculation with successively smaller grid spacings, down to a spacing of 0.2 Å.53 The boundary conditions in each case are determined from the previous calculation. The solute cavity boundary is the standard solvent accessible surface defined by rolling a probe sphere over the solute surface and assigning the solvent inaccessible regions to the solute cavity.49,54 Finally, the solute surface is defined by van der Waals atomic radii, and the cavity dielectric constant, cav, is always taken to be 1. B. Free Energies of Solvation. The free energy of solvation, ∆Gs, is found by subtracting the electrostatic work of charging the solute in a vacuum, Gve , from the electrostatic work of charging the solute in the presence of the fluid, Ge. For a grid-based incompressible fluid model (I), the electrostatic work of charging the solute is
11168 J. Phys. Chem., Vol. 100, No. 26, 1996
Luo and Tucker
GIe(s) ) 1/2∫0sE2(r) dr
(6)
) 1/2∑qiφis
(7)
where the sum runs over all grid points i for which the solute charge qi is nonzero, i.e., only grid points in the molecular cavity (assuming no added salt). The superscript s on the local electrostatic potential φi is to indicate that {φi s} is the solution to Poisson’s equation (eq 2) with the constant solvent dielectric s. Thus, the vacuum result Gve is simply GIe(s ) 1). In the compressible fluid case, the vacuum work term Gve is of course unchanged from in the incompressible case. However, the work term in the presence of the compressible fluid, Ge, is more complex than GIe(s), as a result of the nonlinear relation between the electric displacement vector, D, and the field, E, induced by the compression, i.e., D ) s(E)E. As derived by Frank, the local electrostatic work outside of the solute cavity is given by46
Ge ) ∫[1/20E2(r) + we(r)] dr
(8)
where the first term on the right-hand side is the work required to create the final electrostatic field E(r) in a vacuum and the second term gives the work done on a unit volume of the compressible fluid during the charging process. The integrand factor we(E) is found from the relation46
we(E) ) ∫0 0E′ d[(s - 1)E′] E
(9)
The work done on the fluid per unit mass, we(E), is a property of the fluid and as such can be evaluated for values of E from 0 to Emax and stored for later use in the calculation of Ge. Note that, in the evaluation of eq 9, s and F are evaluated at each value of E′ in the integrand, according to the functions s[F(E)] and F(E) discussed in section II.A. To evaluate the contribution to Ge from outside of the solute cavity (eq 8), E(r) at the grid midpoints, {Ej}, is evaluated from the converged grid-based solution {φi}. The local densities and we values at these midpoints are evaluated from the functions F(E) and we(E), respectively, with E ) Ej. The integration runs over all grid points outside of the cavity. For charged solutes, a long-range Born55 correction is added to Ge to account for the contribution of the field from beyond the solute grid. The work from inside the cavity is calculated from E(r) by using the expression for GIe, eq 6. Note that the simplified expression, eq 7, requires a zero-field boundary condition, and thus it is not valid in this case, because the integral in eq 6 extends only to the cavity boundary where the field is nonzero. III. Hydrolysis of Anisole Klein and co-workers have measured the rate constant for the anisole hydrolysis reaction, eq 1, in SCW at three different state points; see Table 1.13 The three state points considered were all at T ) 380 °C (Tr ) 1.01), but included a range of pressures from 23.2 to 49 MPa. The measured rate constant for hydrolysis, k1, is shown in Table 2 along with kp, the rate constant for the competing pyrolysis pathways. As the pressure is increased from state 1 to state 3, the rate of pyrolysis, kp, decreases by 2 orders of magnitude, while that for hydrolysis increases by more than 1 order of magnitude, making hydrolysis the dominant pathway at high pressure. The measured rate constant for hydrolysis, k1, is a pseudo-first-order rate constant because the solvent provides one of the reactants. The intrinsic bimolecular rate constant, k2 ) k1/[H2O], is also shown in Table
TABLE 1: Properties of SCW at 380 °C state
P (MPa)
F
F/Fc
P/F(∂F/∂P)T
1 2 3
23.2 27.3 49
0.25 0.50 0.63
0.8 1.6 2.0
3.8 10.1 14.4
15 0.81 0.23
TABLE 2: Experimental Rate Constants13 state
kp (×106 s-1)
k1 (×106 s-1)
k2 (×105 L/(mol s))
1 2 3
3.8 3.9
