16136
J. Phys. Chem. 1995, 99, 16136-16143
Continuum Electrostatics Model for an sN2 Reaction in Supercritical Water Gerald E. Bennett? Peter J. Rossky,**’and Keith P. Johnston*$? Departments of Chemical Engineering and Chemistry and Biochemistry, University of Texas, Austin, Texas 78712 Received: June 13, 1995; In Final Form: August 15, 1995@
The electrostatic free energy of solvation AA for the sN2 reaction of the chloride ion with methyl chloride is calculated with a continuum model. The predictive capacity of the continuum electrostatic model is evaluated by comparing the results with simulation. This comparison provides insight into the role of clustering (electrostriction) of solvent about the reacting species. The results for AA are in reasonable agreement with simulation at high and low densities, but the model underpredicts AA significantly in the mid-density region. This behavior is consistent with earlier spectroscopic and simulation studies which demonstrated that clustering (electrostriction) is the most prevalent in the mid-density region.
of the effects of clustering on reactions, the redox reaction of 1. Introduction I2 to I- was investigated electrochemically and modeled with a Supercritical water (SCW) is of interest in power cycle modified Born equation which included electrostriction.20.2’ technology, in hydrothermal breeding of crystals, in extraction However, no spectroscopic data or theoretical models were of organics from soil,’ as a reaction medium in geochemistry, available to test the model. Recently, the excited state deproand in hydrothermal synthesis and waste destruction reaction^.^-^ tonation of /?-naphthol was measured from ambient water to To facilitate development of SCW technologies, a greater SCW.22 The loss of solvation by water above 200 “C strongly knowledge is needed about molecular interactions in SCW inhibits its ability to accept a proton and accelerates the solutions and how they influence chemical reactions, as has reprotonation rate. Under basic conditions where the proton recently been exploned with molecular dynamics computer can be transferred directly to base, water clustering and simulation for an sN2 reaction.6-8 As the nature of SCW reprotonation are much less important. changes from typical aqueous to more typical nonaqueous with Molecular dynamics (MD) computer simulation offers an density, the variation in solvation with temperature and pressure, opportunity to relate the free energy change of reaction to the and consequently changes in reaction rate constants, can be solvation structure at the molecular level. The structure of enormous. supercritical water about various nonionic and ionic solutes has An interesting feature of supercritical fluid solutions is that begun to be characterized by MD ~ i m u l a t i o n . ~Recently, ~~~-~~ solute-solvent interactions increase the local density above the MD simulation has been used to study a chemical reaction in bulk density, so-called solute-solvent clustering, to a much SCW, the S N substitution ~ reaction of C1CH3C1. The greater degree than in the liquid state. Three characteristic Helmholtz free energy change in solution, AA, the internal regions of solvation in supercritical fluids have been defined energy, AE, and entropy, AS, along the reaction coordinate were by Sun et al.:9.10the gas-phase region, the mid-density region, related to the molecular details of solvation of species of varying and the high-density “liquid-like” region. Clustering in the midpolarity.6-8 The conclusion was that the solvating power of density and high-density regions was first characterized with water is remarkably robust. Extreme conditions, e.g. a reduced and in the gas-phase and mid-density absorption density of 0.05 (0.016 g/cm3) at T, = 1, are required to remove regions with a fluorescence probe.I3 Many additional spectrohalf the water molecules from the first solvation shell about scopic studies have demonstrated concepts of clustering with c1-. f l u o r e s c e n ~ e ,Raman,16 ~ ~ ~ ~ ~and ’ ~ EPR1’ ~ ~ ~ measurements. More In comparison to simulation, continuum solvent models for recently, the three regions of solvation have been observed in the calculation of the electrostatic contribution to the solvation SCW on the basis of the x-n* spectral shift of benzophenone free energy require far less computer time. Using a continuum and the n-x* spectral shift of acetone.’* In the low-density values are obtained via a numerical solution to Poisson’s gaslike region, solvation is stoichiometrically d i ~ c r e t e . ~In* ’ ~ ~ model, ~~ equation, which can provide fast and accurate calculations of the high-density liquid-like region, solute-solvent interactions free energies of solvation under ambient conditions.*6 The resemble those in liquid solution. In the mid-density region, results are in good agreement with experiment for various neutral the microscopic picture of solute-solvent interactions is and ionized molecule^.^^-^^ Given the free energy of solvation complex due to the large degree of clustering. The three regions of clustering are directly analogous to the behavior observed in for each species in solution, other thermodynamic quantities adsorption of a solute at a walL83I9 A key goal of the present such as pK,, reaction free energy, and rate of reaction can be obtained. Karplus et al.,30 for example, compared pKa and study is to compare very recent molecular dynamics computer simulation with a continuum electrostatic model to solvation free energy results computed from a continuum model with integral equation results for different ionizable side chains determine how much these local densities influence the free energy change for a reaction in solution. of amino acids. Calculations of pK, for more complex solute The mechanism of solvation has been studied for relatively shapes, such as macromolecules like DNA3’ and protein^,^^,^^ also have been modeled with electrostatic continuum models. few reactions in SCW by in-situ technisues. In a previous studv Recently, continuum electrostatic model contributions to the free ’ Department of Chemical Engineering. energy of activation for the hydrolysis of anisole in supercritical Department of Chemistry and Biochemistry. Abstract published in Advance ACS Absrracfs, October I , 1995 water have been determined.34 For the hydrolysis reaction
+
@
0022-365419512099-16136$09.00/0 0 1995 American Chemical Society
J. Phys. Chem., Vol. 99, No. 43, 1995 16137
S N Reaction ~ in Supercritical Water studied, it was found that at supercritical conditions a large degree of solvent-solute clustering is present. The objective of the present study is to investigate how well the electrostatic continuum theory describes the free energy, enthalpy, and entropy of solvation for a model chemical reaction in SCW over a wide range in temperature and density. Differences between the continuum model and simulation for these quantities will be correlated with simulated local density augmentation (clustering) of water about the reacting species and transition state. Such a comparison will be used in an effort to define thermodynamic conditions where continuum models are sufficiently accurate for predicting free energies of reaction. Another goal is to explore the possibility that the range of usefulness of continuum models may be expanded by using local values of the dielectric constant. We consider the S N substitution ~ reaction of methyl chloride with a chloride ion (eq 1 ) . C1-
+ CH3Cl’
[Cl- - -CH3- - -C1’]
(1)
This reaction has been studied in detail along the reaction coordinate in vacuo and in water at ambient conditions by Monte Carlo ~imulation,~~ integral equation theory,36molecular dynami c ~ and , ~ polarizable ~ continuum calculation^.^^ The recent molecular dynamics simulations of this reaction6-8 in SCW will be used as a basis to test the cohtinuum calculations. For the continuum calculations, we use the implementation of the numerical solution to the Poisson equation contained in the program
2. Theory and Methods 2.1. Free Energy of Reaction in Solution. The free energy surface, or potential of mean force, for the sN2 reaction is given by
Here, rA is the asymmetric reaction coordinate which is obtained by taking the difference between the distances from the methyl carbon to each of the chloride atoms. AA(rA) is the free energy change in the solvent at rA, and the corresponding change in energy for the reaction in the gas phase is designated AUgas(rA). At Ira1 = 8 A, the reactants are far apart and the AA is defined as 0 with AUgasset to 0 for separated reactants. At ambient conditions the potential of mean force at rA = 8 A is less than 1 kcal/m01.~~As the chloride ion approaches the methyl chloride molecule, charge is transferred from the ion to the methyl chloride. At lrAl = 1.44 A, a minimum in the gas phase internal energy is obtained. This state is an ion dipole complex (IDC). Further approach with transfer of charge from the ion to the methyl chloride results in a symmetric transition state complex (TSC) at lrAl = 0 A. At this state, the distance from each chloride to the carbon atom and the charge on each of the chloride atoms are exactly equal, and the free energy is a maximum. 2.2. Continuum Electrostatic Model. In the continuum electrostatic model, the system is represented by a solute which is described by a discrete set of charges qi at positions 0, that are located inside a low dielectric cavity that is imbedded in a solvent. The solvent is represented by a polarizable dielectric continuum with spatially dependent dielectric response function €0. Here, the dielectric constant is taken to be eo outside the cavity, and at the solute/solvent interface, it changes to ~i inside
the cavity. The electrostatic potential solution to Poisson’s equation:
V.E(3)V$(3)
4 0 is obtained from a
+ 4np(3) = 0
(3)
Only for very restricted solute charge distributions and solute cavity shapes can the Poisson equation be solved analytically. For a point charge which lies in a spherical cavity that is imbedded in a medium with a fixed dielectric constant, the Bom equation is obtained.4O Similarly, for a solute described by a set of point charges that lie within either a spherical, elliptical, or cylindrical cavity shape, analytical solutions can be ~ b t a i n e d . ~ IFor - ~ ~other shapes, the Poisson equation may be solved by a finite difference method. A particularly effective implementation is provided by the program DelPhi.39,M.45 In the DelPhi implementation, the continuous functions 40, €0, and p 0 in eq 4 are approximated by discrete values on a three-dimensional grid inside a box. The grid is made up of N nodes in each of three dimensions. The number of nodes (N) is set to 65, which is an optimum resolution for both speed and accuracy.46 For each calculation along the reaction coordinate, the geometry and charge distribution of the [Cl- -CH3- - -C1’] complex are obtained from gas-phase ab-initio calculations performed at the 6-31G* Each atomic partial charge is then placed on the grid by assigning fractional charges to each of the eight vertices of the grid box containing the charge using a trilinear weighting f u n ~ t i o n . 4 ~ ~ ~ ~ Key parameters are the atomic radii of the solute atoms, R, defining the cavity.27 We will consider two different cases here. In the first case, called the “R’ case, the atomic radii R are taken directly from the solute model used in all atom solvent simulations. The Lennard-Jones diameters, cr, are obtained from the optimized potentials for liquid simulations (OPLS) force field, and R = (2)-5’6cr.28.47 In the second case, which we will denote the “R”’ case, the atomic radii are chosen to reproduce experimental data for the free energy of solvation in ambient water using the continuum model. A large data base of R’ values has been p ~ b l i s h e d . ~Of~ particular .~~ note, the atomic radius for the hydrogen atom R’ is adjusted to 1.0 8, .26.29 Similar adjustment to the atomic radius for the hydrogen atom has been made in Poisson calculations for solvation free energies of methanol and acetamide in ambient water to give better agreement with simulated free energy perturbation results.48 For the model-based R case, the set of atomic radii at various positions along the reaction coordinate lrAl are determined directly from the Lennard-Jones parameters used in the simulation (see Table 1). From the reactant to TSC, both chlorines becomes larger, and the R for carbon and hydrogen are kept constant. Here the cavity radius of the chloride ion at the reactant state (rA = 8 A) is 2.612 A. For the R’ case, the cavity radius of the chloride ion is 1.937 A, the value that gives the experimental free energy of hydration in the Born equation at ambient ~onditions.4~ This value is for the reactant state. To obtain the cavity radii of the two chlorine atoms along the reaction coordinate for the R’ case, we chose to fix the relative change from the reactant state to the TSC equal to the relative changes in the R case. As rA is decreased, R’ for the chlorine nucleophile is scaled by the following relationship:
R’(rA)= (1.937/2.612)R(rA)
(4)
In addition to the loci of the atomic charges with their corresponding atomic radii, a solvent probe radius is required by DelPhi. The probe radius is set to 1.48, for water to describe regions of space that are inaccessible to the solvent molecules,29 thus leading to greater accuracy in treating the dielectric
16138 J. Phys. Chem., Vol. 99, No. 43, I995
Bennett et al.
-
-
TABLE 1: Partial Charge and Cavity Radii Used To Describe the Chloride Methyl Chloride Complex [Cl- -CHJ- XI’] along the Reaction Coordinate R’ (A) partial charge le) rA
(A)
8.00 5.71 4.21 2.70 1.44 1.12 0.77 0.54 0.00
c1
C
H
C1‘
c1
C
H
C1’
c1
C
H
C1’
-0.991 -0.989 -0.986 -0.981 -0.973 -0.968 -0.949 -0.912 -0.770
-0.184 -0.185 -0.186 -0.187 -0.189 -0.190 -0.195 -0.204 -0.240
0.125 0.127 0.129 0.132 0.136 0.140 0.151 0.174 0.260
-0.200 -0.207 -0.214 -0.226 -0.247 -0.260 -0.309 -0.405 -0.770
2.612 2.614 2.616 2.620 2.627 2.632 2.648 2.683 2.855
1.946 1.946 1.946 1.946 1.946 1.946 1.946 1.946 1.946
0.927 0.927 0.927 0.927 0.927 0.927 0.927 0.927 0.927
2.624 2.630 2.635 2.645 2.660 2.669 2.698 2.747 2.855
1.937 1.939 1.940 1.943 1.948 1.952 1.964 1.990 2.117
1.946 1.946 1.946 1.946 1.946 1.946 1.946 1.946 1.946
1.000 1.000 1.000 1 .om 1.om 1.om 1.om 1.000 1.om
1.946 1.950 1.954 1.961 1.972 1.979 2.001 2.037 2.117
TABLE 2: Comparison of Dielectric Constant from Simulation with Archer’s Model at the Same Values of Either T and p or T, and R~ dielectric constant state SPC simulation Archer’s dielectric model and Archer’s dielectric model and EOS of Hill at same T, and p:’.52 at T and p25 EOS of Hill at same T and p51.52 P (g/cm3) T (K) 0.125 0.0833 0.0674 0.0579 1.oo
700 800 900 1000 298
2.23 1.65 1.48 1.33 16
boundary on the grid.45,50 The solute parameters used are summarized in Table 1. The electrostatic continuum is described by a dielectric constant inside the solute cavity and E for the solvent. The E inside the cavity was chosen as unity. This is consistent with the MD studies, sihce electronic polarizability is not considered in either case. The S N reaction ~ has been simulated in simple point charge (SPC) water at three different supercritical condition^.^.^ For these conditions of temperature and density p, the bulk dielectric constant has not been determined for the SPC model; however, it has been determined by molecular dynamics simulation at other supercritical c~nditions~~(see Table 2). Table 2 also lists values of the actual dielectric constant of water determined from the equation of state of Hill5’ and the dielectric model of Archere5*The simulation and model values of E are compared at constant T and p as well as constant Tr and pr. Here, the reduced variables are TI = T/Tc and pr = p/pc, where the subscript c denotes critical values. To best match the model and experimental conditions, one has to compensate for the error in the critical properties of SPC water (Tc= 587 K, pc = 0.290 g / ~ m compared ~ ) ~ ~ to that in Hill’s model (Tc = 647 K and pc = 0.323 g / ~ m ~ ) We , ~ Icompare, in Table 2, the predicted value of the bulk solvent dielectric constant for real water and for SPC water on the basis of either equating the actual T and p values or equating the values of Tr and pr, where Hill’s accurate values of Tc and pc are used in each case. As is clear from Table 2, the matching based on reduced parameters rather than absolute values is very successful. Hence, this approach was adopted for obtaining E for the continuum solvent Delphi calculations. Another requirement to solve the Poisson equation in finite difference form is the specification of boundary conditions for the sides of the cubic grid. Initially, Coulombic boundary conditions are used with a cubic box length approximately 5 times larger than the maximum extent of the molecule. The results of this calculation provide the boundary conditions for the next calculation on a finer grid. This method, in which successive calculations are performed and where at each time the grid resolution increases, is termed focusing. Computations of the solvation free energy for the complex at each position along the reaction coordinate were thus performed with initially
2.07 1.54 1.37 1.28 78.67
2.23 1.62 1.42 1.32 78.69
20% of the box filled using Coulombic boundary conditions. Using focusing boundary conditions, further calculations were performed with the percentage of the box filled by the molecule successively increasing from 40% to 60%, 80%, and finally 90% * From the solution of the Poisson equation by Delphi, the electrostatic potential is obtained along all the grid points inside the cubic box. The electrostatic solvation free energy is obtained by summing over all the grid points the product of the potential #i and the corresponding charge, divided by 2. The calculation is also performed in vacuum (EO = l), and subtraction leads to the result:
Note that at ambient conditions the electrostatic Gibbs free energy is essentially equal to the Helmholtz free energy, as the PV term is negligible. At supercritical conditions, where the fluid is compressible, the electrostatic work of charging a volume element is directly related to the Helmholtz free en erg^.^^.^^ Here volume changes upon solvation can be large and the Gibbs free energy can be quite different from the Helmholtz free energy. The entropy change is computed through numerical differentiation of free energy with respect to temperature at constant density:
AS = -a(AA”‘(r,))/aT The intemal energy change (AE) is then obtained simply by taking the difference of AA and -TAS. 3. Results and Discussion
The free energy of the S N reaction ~ in water at ambient conditions is compared in Figure 1 with the molecular dynamics simulations of J ~ r g e n s e n . ~For ~ the particular numerical implementation used here, the uncertainty in AA for Delphi is about 1 kcdmol at each position along the reaction coordinate, which is adequate for our purposes.29 The AA is defined as 0 for r A = 8 A, where the C1- is far from the CH3C1 dipole. The monotonic increase in AA from the reactant state to the TSC reflects the weakening in solvation with charge dispersal. With the OPLS-based values of R , the change in free energy is
sN2 Reaction in Supercritical Water
J. Phys. Chem., Vol. 99, No. 43, 1995 16139
T=25"C p=l &m3
at T ~ 1 . 0 ;~ ~ 4 . ~5 =; 9 . 7 3
Simulation 0 Delphi R size parameter Delphi R' size parameter
20 n I .
n m
aH
0
3
l5
w
10
Y
a
I.\
a
5 0
0
0
2
4
6
25
I
I
8
0
Figure 1. Comparison of the continuum electrostatic model and simulation for the free energy of reaction in ambient water along the reaction coordinate. (Simulation data are from Jorgensen et al.35)
underpredicted relative to the simulations for values of r A < 2. This is a reasonable expectation for a continuum model, since the electrostriction of water about the charged species is neglected. Electrostriction, due to the high electric field near the solute, causes water molecules to become more closely packed than the The electrostriction is greater about the more charged reactant state than the TSC, as shown by simulation.8 Therefore, a larger loss of solvation takes place than expected for a continuum fluid. Consequently, the continuum model underpredicts AA. The difference of 2 kcal/ mol at the TSC will, in fact, have a large influence on the rate constant, as will be discussed just before the Conclusion section. An improved result is observed for the case with the empirical R' cavity size parameters. Because R' for Cl- was fit to give the correct solvation free energy of C1-, and R' was scaled along the reaction coordinate, the results are in much better agreement with simulation along the reaction coordinate, although a significant difference persists at the TSC. The empirically derived R' values compensate for the otherwise neglected effects of a molecular solvent to give an improved solvation free energy. Consequently, the continuum results based on R' appear considerably more accurate than those based on the OPLS values of R. Figure 2 compares the results from Delphi with very recent results from simulation6-8 for the reaction in SCW at three different thermodynamic states. For the case at the highest density (pr = 1.5; Figure 2a) and the critical temperature (T, = l), the simulation results are not very different than in ambient water. The simulated value of AA at the TSC is about 2 kcaV mol smaller than in ambient water. A similar decrease in AA at this condition relative to ambient water is observed for the continuum model for both the R and R' cases. The differences in behavior between the continuum model and the simulations are very similar to those in ambient water, with comparable errors. The situation changes with a decrease in pr to 0.5 as shown in Figure 2b. The Delphi predictions for the R' case are no longer in close agreement with the simulation data. From the reactant state to the ion-dipole complex, the change in solvation is small; only about 2% of the charge originally located on the chloride ion is dispersed to the methyl chloride. Hence, the simulated AA is only 7 kcaVmo1. The distribution functions of solvent about the reacting species change only modestly from
'
at T ~ 1 . 0 ;~~'0.5;~ = 2 . 6 6 n
DelPhiwithR sizeparameten
I
aH
3
Y
W
a 25-
,
i
,
'
I
'
,
at Tp1.3; p,=0.5; -2.24
s 0 aEl 3 8
20 -
0
DelF'hiwithR sizeparsmetera
15
(c)
a I
0 0
I
2
I
,
4 rA
-
I
-
6
8
Figure 2. Comparison of the continuum electrostatic model and simulation for the free energy of reaction in supercritical water. (Simulation data are from Balbuena et al.6-8)
the reactant ~ t a t e . ~Therefore, ,~ the degree to which AA is underpredicted by the continuum model is small. At the transition state, where the charge is most dispersed, the decrease in solvation relative to the reactant state produces a AA of 19 kcaumol. Here the Delphi calculations of AA are notably well below the simulations, by about 4 kcal/mol for the R' case and 7 kcal/mol for the R case. The substantial loss in electrostriction in this lower density case is likely the main cause of the low predictions of AA. Presumably, the R' case takes into account some of the electrostriction, on the basis of the good agreement in ambient water. However, the difference between the R and R' cases is about the same as in ambient water, consistent with the fact that the electrostriction grows in SCW relative to ambient water, as is known from ~ i m u l a t i o n . ~ ~ ~
Bennett et al.
16140 J. Phys. Chem., Vol. 99, No. 43, 1995
TABLE 3: Coordination Number and Local Reduced Density of Water about C1- (Note: V , in Eq 7 Is 419 bi3 at the Reactant State and 299 bi3 for the IDC and TSC) estimated local reduced density thermodynamic conditions modeled coordination number of water around C17.* T, PI Q,P) reactants IDC TSC reactants IDC TSC ~~
1.o 1.o 1.3
1.5
0.5 0.5
9.73 2.66 2.24
7.50 7.20 5.30
8.93 8.16 1.59
Table 3 shows the local coordination number of water molecules around the chloride ion for the reactant, IDC, and TSC at the three thermodynamic state points. The coordination numbers, which were obtained from MD simulation,’.8 represent the number of nearest neighbors in the first solvation shell. The first solvation shell extends to r = 4.75 A, which is the radial distance from the center of the ion to the first minimum after the first peak in the chloride-water (oxygen in water) radial distribution function, gclo(r). The local density ( p ~ is~ ) therefore simply computed by dividing the coordination number by the corresponding estimate of the volume occupied by water in the first solvation shell V,. The computation of p~~ is performed by the following equation:
5.00
2.0
4.90 3.05
1.8 1.7
2.3 2.2 1.6
1.5 1.5 0.9
at Tr=l.O; ~ ~ ‘ 1 . 5 ;~ ~ 9 . 7 3 n
0 Delphi with R size parametas 0 Delphi with R’ size parametas
LII
aP
15
3
6
5l
10
5
0 (7) where p is the bulk solvent density. The disparities between the continuum results (as is best seen for the R’ case) and the simulations are expected to grow monotonically as a function of the degree of electrostriction. This electrostriction may be characterized by the ratio of the estimated local p (Table 3) to the bulk p . At the highest density, the local reduced densities at each point along the reaction coordinate are closest to the bulk reduced densities. Therefore, the effective dielectric response acting on the solute is closest to the bulk dielectric constant. As a result, the continuum results are closest to the simulations at this condition (Figure 1 and Figure 2a). At Tr = 1.0 and pr = 0.5, the difference in the bulk and local p is largest, as is the degree to which AA is underpredicted. The local reduced density at the reactant state is approximately 4 times larger than the bulk density. Increasing the temperature from Tr = 1.0 to Tr = 1.3 at pr = 0.5 causes a large decrease in the coordination numbers, and local density effects become less important. With this temperature increase, the effective dielectric response at the TSC is much closer to the bulk value, such that the solvent continuum model predictions improve relative to the simulation results. In summary, the degree of underprediction of the continuum solvent model is directly related to the electrostriction, as reflected by the difference between the simulated local p and the bulk p. Similar observations are obtained for the solvent internal energy of reaction in solution (AE), as shown in Figure 3a,b,c. Errors in AE are greater than those for AA,due to error resulting from differentiation of AA with temperature at constant density. The Delphi calculations with the R case underpredict the changes in AJT,since electrostriction is neglected. For the R’ case, the Delphi predictions are within statistical error of the simulations, which is quite satisfying for a derivative property. The largest underprediction at the TSC, which occurs for the condition where E = 2.66, is consistent with the result discussed above for AA. In SCW, both the continuum model and simulations give positive values of AS (Figure 4) going from the reactant state to TSC, reflecting an increase in water disorder as charge is dispersed. Because the magnitude of TAS is much smaller than that of AE or AA, and because of large uncertainties resulting from temperature differentiation, the differences between the
at Tr=l.O; p,=0.5; ~ = 2 . 6 6 n
0 Delphi wih R size parametas
L1
aP
0
Delphi wih R’ size parameters
15
3
6
5l
h m
a8
3
6
5l 2
0
4 rA
6
8
(A>
Figure 3. Comparison of the continuum electrostatic model and simulation for the internal energy of reaction in supercritical water. (Simulations references are given in Figure 2.) continuum model and simulation are generally within the uncertainties of the models. Thus, the two cannot be compared to characterize electrostriction. Whereas the simulations were performed at only a few state points, the Delphi calculations, which are much more rapid, have been done for a large number of state points, as shown in Figure 5 for the IDC and TSC. According to the continuum ’
J. Phys. Chem., Vol. 99, No. 43, 1995 16141
S NReaction ~ in Supercritical Water
-----Delphi with R size parameter Delphi with 13’size parameter
at Tr=l.O; p,=1.5; 139.73 . Delphi wih R size parametem
-6
I
.
,
,
,
.
0.0
0.2
0.4
0.6
0.8
1.0
1 1 -E Figure 5. Effect of dielectric constant on the change in free energy of
at Tr=l.O; p,=OS; ~ 2 . 6 6
reaction in solution at the transition state complex and ion dipole complex.
Simulation
o Delphi wih R size parametas Delphi wih R’ size perameters
- 4 -6
y
-4
0
2
4 rA
6
8
Figure 4. Comparison of the continuum electrostatic model and simulation for the entropy of reaction in supercritical water. (Simulations references are given in Figure 2.)
model, the only temperature dependence of AA results from the temperature dependence of E . Therefore, Figure 5 contains all temperatures and pressures for continuum water. The solvent dielectric constant was varied from 2 to 80. The results are presented in terms of (1 - Uc), which serves as a macroscopic measure of the medium polarity from the Born equation.
The results obtained from simulation at the TSC and IDC are also shown. As the dielectric constant approaches 1, the solvent becomes irrelevant and AA goes to 0. For both the R and R‘ cases the curvatures in the predictions of the continuum model
are small. In all cases the continuum model in the R case underpredicts AA from simulation, as already noted. We will interpret the results of Figure 5 in terms of three characteristic regions of solvation in supercritical fluids defined by Sun et al.:9*10 the gas-phase region, the mid-density region, and the high-density region with “liquid-like” solvation. The liquid-like region at high density may be defined for values of e greater than approximately 9.7, or (1 - 1k) > 0.9. Here, pr > 1.5. This region encompasses ambient water, water along the saturation curve, and the high-density part of the supercritical region. Here, the continuum predictions are closer to the simulation data than in the mid-density region, where 0.17 < (1 - lk) < 0.9 or 1.2 < e < 9.7. In this mid-density region, the local densities are further from the bulk density, due to the larger free volume and compressibility of the solvent. This electrostriction is evident in the radial distribution functions from simulation* and in the deviation between the continuum model and simulation for the two lowest simulation points at the TSC in Figure 5. Although the simulations did not explore AA at lower densities, they did examine the radial distribution functions.* On the basis of the distribution functions and the above cited examples of solvation in the three characteristic regions, it may be expected that the continuum model results will again be closer to simulation as the density is lowered into the gasphase region, as solvent clustering decreases. This study provides another validation of a method for examining clustering of solvent about a solute by comparing a spectroscopic (or simulation) result with predictions of a continuum This present comparison is, however, more rigorous than in earlier studies, where the continuum treated the solute as a simple sphere. To obtain the same results as simulation from Delphi for the R case, a higher dielectric constant is required. For example, at pr = 0.5 and T, = 1.0 where E = 2.7, the required dielectric constant is about 5.7. This type of enhancement in the local dielectric constant is consistent with the electrostriction observed in simulation and spectroscopy. It was also considered in the anisole hydrolysis study in supercritical water.34
16142 J. Phys. Chem., Vol. 99, No. 43, 1995
Bennett et al.
TABLE 4: Comparison of Continuum Model and Simulation for the Rate of Reaction Relative to That in an Ideal Gas Wk'g thermodynamic conditions modeled T,
PI
0.46 1.o
3.1 1.5 0.5 0.5
1.o
1.3 a
AA* (continuum model) - AA* (simulation)
with Delphi
€(T,P)
80.0 9.73 2.66 2.24
with Delphi
simulation
R cavity size
R' cavity size
R cavity size
R' cavity size
8.4E-10" 1.3E-4 7.8E-5 2.2E-3
4.3E-9 3.2E-4 4.OE-3 2.4E-2
1.1E-10 6.7E-5 1.3E-3 l.lE-2
-0.96 -1.2 -5.1 -4.0
1.2 0.85 -3.6 -2.7
Read as 8.4 x 1O-Io.
P=350 bar
differences between the model and simulation for the R' cavity set are quite reasonable given that Wk'g varies over so many orders of magnitude. In Figure 6, it is evident that ln(Wk'g) versus 1/T increases linearly at constant density. The same is not true at constant pressure. As the critical temperature is approached, and in the supercritical region, the dielectric constant rapidly decreases with 1/T at constant pressure, leading to highly nonlinear behavior. The simpler behavior of isochores compared to isobars is observed for many other properties of nonelectrolytes and electrolytes in supercritical fluids, for example the equilibrium constant for the reaction of 2-naphthol and OH- in supercritical water.60
+p=0.4 g/cm
4. Conclusion -14.0 -
1.25
1.75
2.25
2.75
lOOO/T (K-') Figure 6. Rate constant relative to the ideal gas value calculated from the continuum model for various isochores and isobars as a function of UT.
Using transition state theory, we can relate the rate constant to the potential of mean force. The rate constant is related to the transmission coefficient K , a preexponential factor A, and the activation barrier for the reaction AW by the following expression: k = KA exp(-AW*lRT)
(9)
where AW* = AW(TSC) - AW(1DC). Sophisticated calculations for K have been performed by molecular dynamics t r e a h ~ e n tand ~ ~ quatum-chemical .~~ free energy surface calculat i o n ~at~ambient ~ conditions. For the transition state theory preexponential factor A , ab-initio calculation^^^ at room temperature have also been made. To show the effects of changes in the PMF with temperature and pressure on the rate of reaction, we have plotted ln(Wkig) using the continuum model versus 1/T (Figure 6) along different isobars and isochores. Therefore, knowledge about K , A, and U as a function of temperature is not required. As a result, eq 9 reduces to ln(k/kig)= -AA*/RT
(10)
where AA* = AA(TSC) - AA(IDC). Because of the exponential relationship in eq 10, small changes in AA* will have a huge effect on the rate constant. For example, at T,= 1.0 and pr = 0.5 the difference between the continuum model and simulation in AA* is 5 kcdmol ( R cavity set); consequently, Wk'g is underpredicted by a factor of 51 (Table 4). At other conditions where electrostriction is less prevalent, the differences in rate constants are smaller. The
The predictions of the simulation data by the continuum electrostatic model are quite reasonable, except in the middensity region. The comparison of the continuum model and simulation provides insight into clustering (electrostriction of solvent about the reacting species). The choice of the cavity radii is extremely important. For a given set of size parameters, the continuum model underpredicts AA from simulation since it does not include electrostriction. The error grows from the high-density to mid-density region as the difference between the local and bulk densities increases. In the R' case, where the size parameters are fit to data in ambient water, the continuum model is more accurate at high densities and at low densities, but still underpredicts the results in the mid-density region. The greater deviation in the mid-density region is consistent with the greater degree of electrostriction in this region, as seen earlier in spectroscopic and simulation studies. In this region, an uncertainty in the free energy barrier of 3-4 kcaymol translates to an uncertainty in the rate constant of a factor of 16. The electrostatic continuum models may be improved by using local dielectric constants, which may be obtained from simulation. Acknowledgments are made to the U S . Army Research Office for University Research Initiative Grant No. DAAL 0392-6-0174 and DAAH 04-93-6-0363 and to the Separations Research Program at the University of Texas. Partial support of this work by a grant from the R. A. Welch Foundation (F0761) to P.J.R. is also acknowledged. We thank the University of Texas System Center for High Performance Computing and Cray Research Inc. for computational support. We thank Barry H. Honig and Anthony Nichols for providing the Delphi program and useful comments, and S . Tucker for providing a stimulating manuscript prior to publication. References and Notes (1) Hawthorne, S . B.; Yang, Y.; Miller, D. J. Anal. Chem. 1994, 66, 29 12-2920. (2) Shaw, R. W.; Brill, T. B.; Clifford, A. A,; Eckert, C. A,; Franck, E. U. Chem. Eng. News 1991.69, 26-39. (3) Tester, J. W.; Holgate, H. R.; Armelli, F. J.; Webley, P. A.; Killilea, W. R.; Hong, G. T.: Bamer, H. E. In ACS Symposium Series; Tedder, D.
S N Reaction ~ in Supercritical Water W., Pohland, F. G., Eds.; American Chemical Society: Washington, DC, 1993; VOI. 518, pp 35-76. (4) Gloyna, E. F.; Li, L. Waste Manage. 1993, 13, 379-394. ( 5 ) Caruana, C. M. Chem. Eng. Progr. 1995, 941, 10-18. (6) Balbuena, P. B.; Johnston, K. P.; Rossky, P. J. J. Am. Chem. SOC. 1994, 116, 2689-2690. (7) Balbuena, P. B.; Johnston, K. P.; Rossky, P. J. J. Phys. Chem. $994, 99, 1554-1565. (8) Flanagin, L. W.; Balbuena, P. B.; Johnston, K. P.; Rossky, P. J. J. Phys. Chem. 1994, 99, 5196-5205. (9) Sun, Y.-P.; Bennettt, G.; Johnston, K. P.; Fox, M. A. J. Phys. Chem. 1992, 96. 10001-10007. (10) Sun, Y.-P.; Fox, M. A.; Johnston, K. P. J. Am. Chem. SOC. 1992, 114, 1187-1194. (11) Kim, S.; Johnston, K. P. AlChE J. 1987, 33, 1603-1611. (12) Yonker, C. R.; Frye, S. L.; Kalkwarf, D. R.: Smith, R. D. J. Phys. Chem. 1986, 90, 3022-3026. (13) Kajimoto, 0.;Futakami, M.; Kobayashi, T.; Yamasaki, K. J. Phys. Chem. 1988, 92, 1347-1352. (14) Brennecke, J. F.; Tomasko, D. L.: Peshkin, J.; Eckert, C. A. lnd. Eng. Chem. Res. 1990, 29, 1682-1690. (15) Knutson, B. L.; Tomasko, D. L.; Eckert, C. A.; Debenedetti, P. G.; Chialvo, A. A. In ACS Symposium Series; Bright, F. V., McNally, M. E. P., Eds.; American Chemical Society: Washington, DC, 1992; Vol. 488, pp 60-72. (16) Ben-Amotz, D.; LaPlant, F.; Shea, D.; Gardecki, J.; List, D. In ACS Symposium Series; Bright, F. V., McNally, M. E., Eds.; American Chemical Society: Washington, DC, 1992; Vol. 488, pp 18-30. (17) Carlier, C.; Randolph, T. W. AlChE J. 1993, 39, 876-884. (18) Bennett, G. E.; Johnston, K. P. J. Phys. Chem. 1994, 98, 441447. (19) Subramanian, R.; Pyada, H.; Lira, C. T. Submitted to lnd. Eng. Chem. Res. (20) Flarsheim, W. M.; Bard, A. J.; Johnston, K. P. J. Phys. Chem. 1989, 93, 4234-4242. (21) Wood, R. H.; Quint, J. R.; Grolier, J.-P. E. J. Phys. Chem. 1981, 85, 3944-3949. (22) Green, S.: Xiang, T.: Johnston, K. P.; Fox, M. A. J. Phys. Chem., in press. (23) Cochran, H. D.; Cummings, P. T.; Karaborni, S. Fluid Phase Equilib. 1992, 71, 1-16. (24) Gao, J. J. Phys. Chem. 1994, 98, 6049-6053. (25) Cui, S. T.; Harris, J. G. Chem. Eng. Sci. 1994, 49, 2749-2763. (26) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 19781988. (27) Mohan, V.; Davis, M. E.; McCammon, J. A,; Pettitt, B. M. J. Phys. Chem. 1992, 96, 6428-6431. (28) Jean-Charles, A,; Nicholls, A,; Sharp, K.; Honig, B.; Tempczyk, A.; Hendrickson, T. F.; Still, W. C. J. Am. Chem. Soc. 1991, 113, 14541455. (29) Simonson, T.; Briinger, A. T. J. Phys. Chem. 1994, 98, 46834694. (30) Lim, C.; Bashford, D.; Karplus, M. J. Phys. Chem. 1991,95, 56105620.
J. Phys. Chem., Vol. 99, No. 43, 1995 16143 (31) Warwicker, J.: Watson, H. C. J. Mol. Biol. 1982, 157, 671-679. (32) Gilson, M. K.; Honig, B. H. Proteins 1988, 4, 7-18. (33) Gilson, M. K.; Honig, B. H. Proteins 1988, 3, 32-52. (34) Tucker, S. C.; Gibbons, E. M. In ACS Symposium Series; Truhlar, D. G., Cramer, C. J., Eds.; American Chemical Society: Washington, DC, 1994; Vol. 568, pp 196-211. (35) Chandrasekhar, J.; Scott, F. S.; William, L. J. J. Am. Chem. SOC. 1985, 107, 154-163. (36) Huston, S. E.; Rossky, P. J.; Zichi, D. A. 1.Am. Chem. SOC. 1989, 111, 5680-5687. (37) Hwang, J.-K.; King, G.; Creighton, S.; Warshel, A. J. Am. Chem. SOC. 1988, 110, 5297-5311. (38) Basilevsky, M. V.; Chudinov, G. E.; Napolov, D. V. J. Phys. Chem. 1993, 97, 3270-3277. (39) Sharp, K. A. Delphi; 3.0 ed.; Sharp, K. A,, Ed.; Columbia University, 1988. (40) Born, M. Z. Phys. 1920, I , 45-48. (41) Beveridge, D. L.; Schnuelle, G. W. J. Phys. Chem. 1975, 79, 25622566. (42) Harrison, S. W.; Nolte, H.-J.; Beveridge, D. L. J. Phys. Chem. 1976, 80, 2580-2585. (43) Felder, C. E. J. Chem. Phys. 1981, 75, 4679-4682. (44) Honig, B.; Sharp, K.; Yang, A.4. J. Phys. Chem. 1993,97, 11011109. (45) Gilson, M. K.; Sharp, K. A. J. Comput. Chem. 1987,9, 327-335. (46) Klapper, I.; Hagstrom, R.; Fine, R.; Sharp, K.; Honig, B. Proteins 1986, 47, 47-59. (47) Sharp, K.; Jean-Charles, A,; Honig, B. J . Phys. Chem. 1992, 96, 3822-3828. (48) Still, W. C.; Temczyk, A.: Hawley, R. C.; Hendrickson, T. J . Am. Chem. SOC. 1990, 112, 6127-6129. (49) Rashin, A. A.; Honig, B. J. Phys. Chem. 1985, 89, 5588-5593. (50) Lee, B.; Richards, F. M. J. Mol. Biol. 1971, 55, 379-400. (51) Hill, P. G. J. Phys. Chem. Ret Data 1990, 19, 1233-1274. (52) Archer, D. G.; Wang, P. J. Phys. Chem. Re$ Data 1990, 19, 371411. (53) dePablo, J. J.; Prausnitz, J. M.; Strauch, H. J.; Cummings, P. T. J. Chem. Phys. 1990, 93, 7355-7359. (54) Frank, H. S. J. Chem. Phys. 1955, 23, 2023-2032. ( 5 5 ) Wood, R. H.; Smith-Magawan, D.; Pitzer, K. S.; Rogers, P. S. 2. J. Phys. Chem. 1983, 87, 3297-3300. (56) Bockris, J. 0. M.; Reddy, A. K. N. Modern Electrochemistry, 1st ed.; Plenum Press: New York, 1970. (57) Kim, S.; Johnston, K. P. lnd. Eng. Chem. Res. 1987, 26, 12061213. (58) Bergsma, J. P.; Gertner, B. J.; Wilson, K. R.: Hynes, J. T. J. Chem. Phys. 1987, 86, 1356-1376. (59) Gertner, B. J.; Wilson, K. R.; Hynes, J. T. J. Chem. Phys. 1989, 90, 3537. (60) Xiang, T.; Johnston, K. P. J. Phys. Chem. 1994, 98, 7915-7922.
JP95 I6 194
