J. Phys. Chem. B 1997, 101, 1063-1071
1063
A Compressible Continuum Model Study of the Chloride plus Methyl Chloride Reaction in Supercritical Water Huaqiang Luo and Susan C. Tucker* Department of Chemistry, UniVersity of California, DaVis, California 95616 ReceiVed: September 18, 1996; In Final Form: NoVember 25, 1996X
The recently introduced compressible electrostatic continuum solvation model [J. Phys. Chem. 1996, 100, 1165] is found to accurately reproduce molecular simulation predictions [Flanagin et al. J. Phys. Chem. 1995, 99, 5196] of the solvation energies and nonstructural density distributions for the title reaction under supercritical conditions. The compressible continuum model is unique in its ability to isolate the importance of the solventsolute clustering in defining the energetics of this reaction at different thermodynamic states. The thermodynamic conditions required for application of the model are also discussed.
I. Introduction Supercritical fluids are promising solvents for use in industrial chemical processes, both because of the novel reactivity observed in these solvents and because of the thermodynamically controlled tunability they exhibit.1-12 In particular, at thermodynamic conditions near the supercritical fluid solvent’s critical point, small changes in pressure cause large changes in the solvent’s density and correspondingly large changes in its solvating properties. In supercritical water (SCW), for example, the dielectric constant of the fluid may be varied from nonpolar values (∼2) to polar values (∼14) with a change in pressure of less than 100 atm.13 Supercritical fluid solvent tunability is greatest in the “compressible regime”, i.e., in the part of the phase diagram beyond the critical point where the isothermal compressibility of the fluid, κT, is large. Thus, it is in this compressible regime that we would like to understand how the supercritical fluid solvent affects solute reaction dynamics and solute reactivity. Theoretically, the gross effect of a solvent on a solute reaction can be characterized by the change in the reaction’s activation barrier upon solvation.14-16 That is, the free energy of solvating the solute may vary along the solute’s reaction path. Such differential solvation along the reaction path will alter the relative stability of the solute species along this path, consequently altering the reaction’s activation barrier. Thus, the necessary quantity for understanding solvent effects is the solute’s reaction-path-dependent free energy of solvation. The most accurate theoretical tools for computing relative solvation free energies are Monte Carlo and molecular dynamics (MD) computer simulation,17 and these methods have proven to be effective in normal liquids.18-20 Unfortunately, these techniques become costly and inefficient in the compressible regime around the solvent’s critical point where long-range correlations in the solvent’s density fluctuations become important. Note that it is just these long-range correlations which are at the origin of the large isothermal compressibilities discussed above.21 Monte Carlo or MD treatment of these longrange correlations requires large simulation cells and hence large numbers of particles.22 Additionally, such long-range correlations mean slow relaxation times, and these large simulations would need to be carried out over very long times.22,23 It is thus of significant interest to explore the use of other, less * To whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, January 15, 1997.
S1089-5647(96)02860-X CCC: $14.00
computationally intensive methods for examining solvation energies in supercritical fluids. Dielectric continuum solvation models provide a promising alternative to computer simulation.24-27 These models are simple and computationally efficient, and they have been used successfully in the study of many solution-phase reactive systems under ambient conditions. These methods rest on the premise that the solvation free energy is just the difference in the electrostatic work of charging the solute in the presence of the solvent and that of charging the solute in a vacuum, although corrections are sometimes added for other solvation effects, such as cavitation or dispersion.25,27 To evaluate this electrostatic work, the solute is represented by a set of point charges in a vacuum cavity, while the solvent is represented by a continuous dielectric medium surrounding the solute cavity. The choice of solute cavity is not unique, although it is usually defined by interlocking atom-centered spheres. Sets of atomic radii have been developed for the purpose of defining these molecular cavities.25-28 Unfortunately, by using a single, fixed dielectric constant to represent the solvent, conventional continuum models assume, de facto, that the solvent retains a fixed, uniform density despite the presence of the solute. For a compressible fluid, electrostrictionsthe process by which the solvent compresses toward the solute in response to the solute-generated electrostatic field (E)sbecomes important. The fact that electrostriction is nonnegligible, even in liquids under ambient conditions, was demonstrated by Rick and Berne through a comparison of incompressible continuum calculations with simulation.29 An important feature of systems in which electrostriction plays a role is that the effective dielectric constant of the solvent becomes nonuniform as the solute is charged. A few previous continuum solvation studies have considered nonuniform dielectric constants, but their focus has been primarily on the treatment of first solvation-shell effects.30-34 One such model,34 which involves dividing the space outside of the solute cavity into a near and a far region and assigning a different value of the dielectric constant to each region, has been applied to the electrostriction problem in SCW.4 However, the required parameters of the modelsthe boundary position and inner dielectric constantsdepend upon the electric displacement (D ) (r)E) generated by the solute and had to be taken from computer simulation. Recently, we proposed a continuum model which includes an explicit, self-consistent treatment of field-induced solvent compression and the resultant nonuniform © 1997 American Chemical Society
1064 J. Phys. Chem. B, Vol. 101, No. 6, 1997 dielectric constant.36,37 This compressible continuum model, which is applicable to molecules and thus to reacting solutes, is based upon an earlier model developed by Wood et al.38 for atomic ions. Since electrostriction increases in importance with increasing solvent compressibility, this phenomenon will obviously be much more significant for highly compressible supercritical fluids than for their liquid counterparts. Tucker and Gibbons recently illustrated this point, demonstrating that for the anisole hydrolysis reaction in SCW the incompressible continuum (IC) model predicts highly inaccurate solvation energies and, as a result, does not reproduce the experimental pressure dependence of the reaction rate for this system.12 Subsequently, we applied our compressible continuum (CC) model to this reaction system and found that this model provides a much better description of the pressure dependence of the reaction rate when compared to experiment.11,37,39 However, the experimental data are extremely limited, providing only the rate constant at three different pressures.11,40 Thus, no check on the detail of the CC calculationsssolvent density distributions and absolute differentials in solvation energy along the reaction pathscould be made.37 Johnston and Rossky and co-workers have provided a more detailed demonstration of the failure of the IC model by comparing results from this model against computer simulation.41 Specifically, these authors conducted intensive MD simulations of the Cl- + CH3Cl exchange reaction in ambient water and in SCW at three different state points.42,43 From each of these simulations, they evaluated the reaction path dependence of the solvation free energy. Comparison of the differential solvation along the reaction path predicted by the MD simulations with those predicted by IC model calculations showed that, indeed, the IC model fails dramatically for this reaction under certain supercritical conditions.41 We now examine the ability of our CC model to describe the Cl- + CH3Cl system in SCW. This system provides a stringent test for the CC model, both because the magnitude of the errors in the IC model are large and because the MD simulations provide detailed information about the solvent density distribution that the CC model must reproduce. In section II we present a brief overview of the CC model. The parameters of the model system are reviewed in section III, while resultssfor both solvent density distributions and solvation energeticssare presented in section IV. Conclusions follow. II. Compressible Continuum Model The compressible continuum (CC) model provides a method for determining the electrostatic work, Fel, of charging the solute in the presence of a compressible dielectric fluid. Because the fluid being considered is compressible, processes conducted under conditions of constant volume and of constant pressure could differ substantially. It is therefore important to determine whether the electrostatic work calculated by this model corresponds to the Helmholtz or to the Gibbs free energy. This question is considered in detail in the Appendix, where we argue that the electrostatic work evaluated from the present formulation of the CC model is most appropriately equated with the Gibbs free energy, but that, under the assumed condition of infinite dilution, it may be equated with the Helmholtz free energy as well. In the CC model, the solute is represented by a set of point charges Fe(r) in a molecule-shaped cavity, and the surrounding solvent is represented by a compressible dielectric fluid. The total electrostatic energy of charging, Fel, is37,38,44,45
Luo and Tucker
Fel )
∫dr {1/20E2(r) + δc∫0E′)E(r)0E′ d[(s(E′) - 1)E′]} (1)
where the integral dr is over all space including the solute cavity, the first term in brackets is the in Vacuo contribution from each volume element of space, and the second term gives the work done on each volume element of the compressible fluid upon charging. (The relation of this second term to eq 7 is made by realizing that, for unit specific volume of the fluid, VP ) 0(s - 1)E.) The factor δc indicates that the second term is zero in the solute cavity, within which there is no fluid. To evaluate the electrostatic energy by eq 1 requires knowledge of the position-dependent E field, which is governed by Poisson’s equation for the electrostatic potential φ(r), i.e.
0∇‚[(r)∇φ(r)] ) -Fe(r)
(2)
Note that here (r) equals 1 for points inside the molecular cavity and s(r) for points outside. The potential φ(r), and thus the field E(r), are seen to depend upon the nonuniform solvent dielectric constant s(r), which itself depends upon the fluid density at r, F(r). But F(r) depends on the local field E(r), via eq 10, and hence we must find a self-consistent solution of eq 2 with s(r) ) s[F(E(r))].37,38 Given s[F(E)], which may be evaluated in advance, self-consistent solution of eq 2 with electrostriction is obtained with the iterative finite grid algorithm46 described in detail previously.37 The function s[F(E)] is comprised of two parts. The field dependence of the solvent density F(E) is found from an integration of eq 10 started at the bulk density, F(E ) 0). Note that, when integrating the right hand side of eq 10, F, κE,T, and (∂s/∂F)E,T are all functions of the field E. The field dependence of the compressibility is calculated as suggested by Wood, from an integration started at κT with E ) 0,47 i.e.
( )∫
∂2s 1 κT(E) ) κT(E)0) + 0F2 2 2 ∂F
E′)E
E,T
0
κT2(E′) dE′2 (3)
where κT(E)0) is determined from an empirical function.48 Second, the density dependence of the dielectric constant s(F), which is also required for the first and second derivatives of s relative to F in eqs 10 and 3, is evaluated from an empirical function known only in the zero-field limit.49 Using this zerofield approximation amounts to neglect of dielectric saturation, the effects of which are expected to be significantly smaller than the electrostriction effects which are included by our treatment.29 We note, however, that in ref 45 it was found that dielectric saturation effects become important for low bulk density, low bulk dielectric conditions under which soluteinduced E fields may have large values. Thus, these effects are possibly non-negligible under the low-density conditions studied herein, and we intend to address this issue in future work. III. Parameters for Cl- + CH3Cl in Supercritical Water We consider the Cl- + CH3Cl SN2 exchange reaction in SCW because simulation results are available for this system.18,42,43 Thus, the parameters used to model this system have been chosen to correspond to the simulated system. A. Continuum Model Parameters. The parameter sets required for application of a continuum model to this system are the solute geometry, charge distribution, and molecular cavity at various points along the reaction path, which for this system is just the difference in the two Cl-C distances, rA. The symmetric transition state is at rA ) 0, there is an ion-dipole
Cl- + CH3Cl Reaction in Supercritical Water complex at rA ) 1.44 Å, and rA ) 8 Å is taken to represent the reactants. We use exactly the same parameter sets as were used in the IC model calculations of this system by Benett et al.,41 but we detail them below, in order to clear up some possible points of confusion. Geometries and Charges. Specifically, the all-atom geometries and charge distributions used herein were taken directly from the HF/6-31G* calculations,50 to which the MD simulation model of this reaction was fitted.18,51 We note that in the simulation study of Flanagin et al. a united-atom methyl group was used,18,42 and as a result of compensation for this approximation, the fitted three-site charge distributions used in these simulations are not identical to those achieved by summing the ab initio values. Solute Radii. The solute cavity is taken to be the solvent accessible surface (probe radius of 1.4 Å) created by the interlocking atomic spheres whose radii are taken as follows: the reaction path-dependent Cl radii are taken directly from the Lennard-Jones values used in the simulation,42 while the C and H radii are taken from OPLS Lennard-Jones values.52,53 This is the R parameter set of ref 41, and the values are tabulated there. Bennett et al.41 also considered a second radius set, R′, which was optimized to accurately predict the experimental solvation free energy of an isolated chloride ion in water under ambient conditions. These authors compared the relative free energies of solvation along the reaction path computed from the IC model with each of the two radius sets in water under ambient conditions, i.e., s ) 80, with those computed from simulation under the same conditions.18 From simulation, the solvation of the transition state was found to be 20.8 kcal/mol less favorable than that of the reactants, whereas the IC continuum model predicts 18.6 kcal/mol with set R and 23.0 kcal/mol with set R′. We repeated these calculations using the CC model and found that electrostriction effects enhance the chloride ion solvation slightly more than that of the transition state. The resultant differential is 19.1 and 23.5 for the R and R′ parameter sets, respectively. Thus, with electrostriction explicitly included, the R parameter set result is in better agreement with simulation (within 1.7 kcal/mol) than is the R′ set result (off by 2.7 kcal/mol). Additionally, it is perhaps not appropriate to use the R′ set with explicit electrostriction as it has been suggested that these radii have been adjusted to compensate for neglected electrostriction effects.29,41 Finally, the R′ radii set is not related to the radii used in the simulations, nor was it optimized along the reaction path.54 We therefore used the R, rather than the R′, parameter set in all subsequent CC model calculations. Finite Grid Parameters. The parameters for the numerical solution of Poisson’s equation (eq 2)sgrid spacing, boundary conditions, etc.sare exactly as described in our previous study.37 The numerical error associated with using such a finite grid solution is ≈(0.5 kcal/mol. However, when comparing CC and IC model calculations for a single solute configuration, the same grid is used, and this grid-based error tends to cancel. As a result, comparisons of CC and IC results are more accurate than the standard error bars would imply. B. Properties of Supercritical Water. Relation of Simulated and CC Model Water. The simulations were performed at three supercritical state points of SPC water.55 The dielectric propertiessand presumably also other thermodynamic propertiesscomputed from this model of water agree well with those for real water when the state variables are corrected for the difference betweeen the critical variables, Tc and Fc, of SPC and real water.41,56 That is, SPC water at reduced variables Tr ) T/Tc and Fr ) F/Fc is comparable to real water at the same
J. Phys. Chem. B, Vol. 101, No. 6, 1997 1065 TABLE 1: SCW States Studieda state
F (g cm-3)
F/Fc
T/Tc
s
b κT/κIG T
1 2 3
0.485 0.165 0.162
1.50 0.51 0.50
1.0 1.03 1.3
9.73 2.66 2.24
1.28 3.28 1.34
All quantities from ref 48 except s from ref 49. b κIG T is the compressibility of an ideal gas under the same thermodynamic conditions. a
reduced variables. This correlation enables us to use experimentally determined functions for the density dependence of the dielectric constant and the compressibility in the CC model and still compare these calculations with the simulations. The bulk dielectric constant and the reduced compressibility are given in Table 1 for these three SCW states. CC Model with T ) Tc. The CC model calculations were performed at the same three state points as were the simulations, except that slightly shifted conditionssTr ) 1.03, Fr ) 0.51 instead of Tr ) 1.00, Fr ) 0.50swere used for state 2. The shift in temperature was required because, for compressible calculations with Tr ) 1.00 and a bulk density Fr < 1, regions of the solvent may be compressed to exactly critical conditions (Tr ) 1, Fr ) 1), where the zero-field compressibility, κT,E)0, diverges, potentially causing catastrophic failure of the model. This divergence is nonphysical in that it arises because we apply the compressibility of a macroscopic sample of the pure solvent at Tr ) 1, Fr ) 1 to a small volume element of solventsthe size of which determined by the grid spacing used in the numerical solution of Poisson’s equationsunder the same thermodynamic conditions. The divergence of κT in macroscopic systems is a result of the divergent correlation length of the solvent density fluctuations. If the natural correlation length exceeds the linear dimensions of the volume, the correlation length will be limited by the size of the system, and as a result, κT will not diverge.57 While this strong dependence of κT on system size occurs only very close to the critical point, it will be interesting to build this size dependence of κT into the CC model. Finally, note that the density of state 2 is also shifted slightly, so that the dielectric constant remains the same in the shifted as in the unshifted state. The CC model results for the shifted state are thus directly comparable with the IC model results for the original state. IV. Results and Discussion A. Solvent Density Distribution. 1. Density Contours. Unlike the usual IC models, the CC model allows the solvent density to respond to the solute-induced electric field, generating the nonuniform, equilibrium density distribution of solvent around the solute. It thus enables us to examine the “solventclustering” phenomena in supercritical fluids.4,9,37,58-63 A contour plot of the solvent density distribution around the iondipole and transition state complexes is given in Figure 1 for the SCW state point 3. The constant-density contoursswhich in fact also correspond to lines of constant fieldsare of very different shape around the two solute configurations. For the ion-dipole complex, which has virtually all of the excess negative charge centered on one of the chlorine atoms (qCl ) -0.97 e), the compression is highly asymmetric, being localized around the charged chlorine atom. In contrast, the compression around the transition state complex is symmetric across the CH3 plane, being centered equally around the two chlorines, each of which carry a charge of -0.77 e. In both cases the density contours are concave near the CH3 group where the charge intensity and the resultant E field are relatively weak. As the distance from the solute is increased, the contours become less
1066 J. Phys. Chem. B, Vol. 101, No. 6, 1997
Luo and Tucker
Figure 2. Radial density distribution around Cl- for states 1 and 3. CC (solid line); MD (dashed line).
Figure 1. Solvent density contours in the molecular plane for SCW state 3: (a) ion-dipole complex; (b) transition state. Solute atoms are denoted by different size solid circles. Largest to smallest represent C, Cl, and H atoms, respectively. Density is in g cm-3.
concave until they finally converge to spherical surfaces at distances which are sufficiently large that the solute acts as an effective point charge. We note that the contours around the ion-dipole complex approach spherical at much shorter distances from the solute than do the contours around the transition state. This difference is expected, despite the fact that the transition state is more compact, since the charge on the ion-dipole complex is more centralized than on the transition state. 2. Radial Density Distributions. Here we compare the density distributions predicted by the CC model with those from MD simulations.42,43 In particular, we consider the radial density distribution functions, F(r) ) Fbg(r), where Fb is the bulk density of the state point and g(r) is the usual radially averaged sitesite pair distribution function. In the simulation work, the Cl-O distribution function was evaluated, since in SPC water the oxygen carries the Lennard-Jones parameters and defines the molecular position. In order to compare these SPC-based density distributions with those from the real-water-based continuum models, the SPC bulk density is scaled to the real f water density for the same reduced variables, i.e., FSPC b SCP H2O (FSPC /F )F . For the CC model, F(r) is found simply by b c c radially averaging the full distribution, F(r), around the chloride ion. Chloride Ion; CC Vs MD. Figure 2 shows the radial density distribution around an isolated chloride ion from both the CC
Figure 3. Integrated excess number of water molecules around Clfor states 1 and 3.
model and the MD simulations for SCW state point 1 and 3. While the continuum model cannot reproduce the molecular structure appearing in the MD distribution function, it does a good job of predicting the underlying nonstructural density enhancement around the solute. In particular, for r > 4.0 Å, i.e., beyond the nearest-neighbor peak, the slow decay of the MD density distribution toward its bulk value is reproduced quite well by the continuum model, and we observe no tendency for this model to overestimate the range of the solvent density enhancement. Wood and co-workers45 have performed a similar comparison of the density distribution around a spherical ion in SCW computed from their CC model and from molecular simulation using TIP3P water.64 While these authors did observe a much longer ranged density enhancement from the CC model than from the simulation at 670 K and 0.32 g/cm3, they attributed this discrepancy to the fact that the simulation results were not corrected for the probable deviation of the TIP3P critical point away from that of real water. That our comparison with corrected simulation shows such good agreement supports this hypothesis. Figure 3 shows the integrated excess number of solvent molecules around the chloride ion, evaluated from the CC model and MD simulation density distributions. Comparing the CC model results for state 3 with those from the MD simulation, it is evident that, especially at short range, the CC model underpredicts the excess number of solvent molecules, presumably because specific molecular interactionsswhich produce
Cl- + CH3Cl Reaction in Supercritical Water increased solvent coordination numbers even in nonpolar solventsshave been neglected. At longer range the CC curve roughly parallels the MD curve, indicating that outside of the range of these specific molecular interactions the CC model does a reasonably good job of predicting the compression. While the excess numbers of solvent molecules for state 1 and state 3 are very similar, the excess number for state 1 is less than that for state 3 at short range and exceeds that for state 3 at longer range. The CC model predicts the relationship between the excess numbers for these two states quite well, accurately reflecting both the magnitude of the difference and the point of crossing (∼9 Å). This suggests that the CC model accurately predicts the underlying, nonmolecular compression effect at all ranges and that it is this underlying behavior which determines the thermodynamic state dependence of the solvent compression effects. Chloride Ion; State Dependence. It is also of interest to compare the CC model density distributions shown in Figure 2 for the two SCW states 1 and 3, because the reduced bulk compressibilities at these states are similar (see Table 1), but their bulk densities and bulk dielectric properties differ. The local density of the low bulk density state 3 approaches that of the high bulk density state 1 as r decreases toward the solute, reflecting the fact that the compression is more substantial for the low bulk density state 3 (and also explaining why, at short range, the excess number for state 3 exceeds that for state 1 in Figure 3). The difference in the compression for these two states is magnified when we consider the relative density enhancement for each state. For example, the local density 5 Å from the solute center is about twice the bulk density for the low bulk density state 3, and only 20% greater than the bulk density for the higher bulk density state 1. The greatly enhanced compression of state 3 relative to state 1 can be understood as follows. The bulk density of state 3 (Fr ) 0.5) is below the critical density, and as a result, local compression brings the local density closer to the critical density (at least until the local F exceeds Fc), thus increasing the compressibility of that local volume and assisting additional compression. In contrast, the bulk density of state 1 (Fr ) 1.5) is above the critical density, and hence local compression moves the local density away from the critical point, thus decreasing the local compressibility and hindering further compression. That low bulk density states will exhibit greater compression than those of higher density has already been observed experimentally by Carlier and Randolph.60 These authors found that the maximum solvent density enhancement around a solute occurs for supercritical states having bulk densities of about 1/2Fc to 2/3Fc, rather than for those with F ) Fc where the bulk compressibility is maximized. This behavior, which has also been observed in simulation studies,42,65 has now been given a fundamental basis,59 see below. Ion-Dipole Complex; CC Vs MD. In Figure 4 we compare radial density distributions from the CC model and MD simulations for the ion-dipole complex at SCW states 2 and 3. The radial distributions for the ion-dipole complex are centered on the more negatively charged chlorine atom. As a result, for small r the spherical volume element 4πr2 dr contains part of the solute cavity, which is inaccessible to the solvent, and the radially averaged value F(r) is lower than the local density outside of the solute cavity. This excluded cavity effect is the cause of the nonuniformity of the IC radial density distribution and also of the slight structure present in both the CC and IC distributions. Comparing the density distribution from the CC and MD calculations for state 3, we find the relationship to be similar
J. Phys. Chem. B, Vol. 101, No. 6, 1997 1067
Figure 4. Radial density distribution around the negatively charged Cl in the ion-dipole complex: (a) state 3, (b) state 2.
to that found for the chloride ion for this same SCW state. Comparing with the IC model results, we find that the average local density still exceeds that of the bulk at 12 Å. The density distributions for state 2 illustrate important properties of this supercritical state, so they are discussed in this context below. Ion-Dipole Complex; State Dependence. Comparing the CC density distributions for state 2 and state 3, which must both decay to approximately the same bulk density, it is evident that both the local density enhancement and the range of this enhancement are much greater for state 2. This relationship is expected because the temperature of the state 2 is closer to the critical value than that of state 3, and consequently the reduced compressibility of this state is larger. Inspection of the CC and MD density distributions for state 2 suggests that the density enhancement predicted by the MD simulation decays to the bulk value more slowly than does that predicted by the CC model, although noise in the simulation results precludes a definite conclusion. Such a deviation is indeed expected for states near to the critical point, as “critical phenomena" which are not incorporated into the CC model become important. This point is most readily understood by considering the partial molar volume of solvation, which is a measure of the solute-induced density enhancement. As shown by Chialvo and Cummings,59 the partial molar volume is comprised of two components. The first component depends upon the direct solvent-solute correlation function and is thus understood to be a direct result of the attractive solvent-solute potential interactions. In electrostatic systems, this phenomenon is just electrostriction, and both the CC model and the MD simulations account for it. The second component, which is proportional
1068 J. Phys. Chem. B, Vol. 101, No. 6, 1997
Luo and Tucker TABLE 2: Free Energies of Solvation (in kcal/mol) IC rA (Å) state 1 8.00 1.44 0.00 state 2 8.00 1.44 0.00 state 3 8.00 1.44 0.00
Figure 5. Compression-induced solvation free energy.
to the compressibility, represents the phenomenon whereby the local density enhancement generated by the first, direct component itself induces density enhancements over a longer range as a result of long-range solvent-solvent density correlations which diverge as the critical point is approached. This second phenomenon is not included in the CC model, while it is included in the MD simulations, at least to the extent that the simulation box length is large enough to encompass these longrange correlations. Thus, we attribute the apparently extended range of the MD density distribution relative to the CC distribution to the solvent-solvent correlation effects which are not taken into account in the CC model. Along these same lines, Wood and co-workers found, for SCW at Tr ) 1.03 and Fr ) 1.0, that the CC model grossly underestimates the partial molar volume of NaCl, indicating that the solvent-solvent correlation effects neglected by the CC model dominate the activation volume at this near critical state.45 B. Reaction-Path-Dependent Solvation Energetics. Compression-Induced SolVation Energy. In order to examine the effect of electrostriction on the solvation energies for the Cl+ CH3Cl reaction, we introduce the compression-induced solvation free energy,
∆Fc ) ∆Fel - ∆FIC el
(4)
which is just the difference between the solvation free energy calculated with and without solvent compression. This quantity is plotted in Figure 5 for various points along the Cl- + CH3Cl reaction path for each of the SCW states. For states 2 and 3, the compression-induced stabilization of the reaction complex varies in magnitude along the reaction path from about 3 kcal/ mol at the transition state to over 10 kcal/mol at 8 Å separation. The effect is much smaller for state 1, in agreement with the the lesser degree of compression observed for this state in Figure 2. That the degree of electrostriction and the resultant compression-induced stabilization are greater for the charge-localized reactant complex than for the more delocalized transition state is consistent with our previous findings for the anisole hydrolysis reaction.37 Note that the most compressible State 2 also exhibits the greatest reaction-path-dependent change in the compressioninduced stabilization. CC Vs MD SolVation Energies. Table 2 gives the absolute solvation energies for three points on the reaction path, the transition state, ion-dipole complex, and reactant, calculated by both the CC and IC models,66 for all three SCW states. We see that the decrease in compression-induced stabilization as the reaction proceeds from the reactant to the transition state
CC
∆Fel
∆Fr
∆Fel
∆Fr
MD ∆Ar
-57.2 -50.3 -40.1
0.0 6.9 17.1
-59.1 -50.8 -40.4
0.0 8.3 18.7
0.0 7.9 19.1
-39.9 -34.7 -27.4
0.0 5.2 12.5
-50.2 -42.1 -30.6
0.0 8.1 19.6
0.0 7.3 19.4
-35.4 -30.9 -24.8
0.0 4.5 10.6
-43.9 -37.8 -27.8
0.0 6.1 16.1
0.0 6.0 15.9
compounds an already substantial decrease in the incompressible solvation energy along this path. Also shown in Table 2 are the computed CC model and IC model relative solvation free energies, ∆Fr(rA) ) ∆Fel(rA) - ∆Fel(rA ) 8 Å), along with the MD values42 for the same quantities (which are labeled ∆Ar in Table 2 because these calculations were performed at constant N, V, T ; see the Appendix). The relative solvation energies, ∆Fr(rA), are shown as a function of the reaction coordinate in Figure 6. At each SCW state point, the reaction path dependence of the free energies predicted by the CC model is in excellent agreement with those from simulation. In fact, the largest difference is only about 1.3 kcal/mol, with the average deviation from the MD results being about 0.5 kcal/mol, well within the error bars of the CC method. The CC model provides a major improvement over the IC model, which consistently underestimates the reaction path dependence of the solvation free energies. For example, the IC model yields relative free energy increases which are 2.0, 6.9, and 5.3 kcal/mol too low at the transition state for SCW states 1, 2, and 3, respectively, while the corresponding CC model predictions are off by only 0.4, 0.2, and 0.2 kcal/mol. Our results demonstrate that, as was suggested by Bennett et al.,41 the poor predictions of traditional IC models for this reaction in SCW are indeed a result of electrostriction, which is neglected in those models but is treated quite accurately by our CC solvation model. The astute reader may be wondering whether the CC model should predict such accurate free energies for SCW state 2 if the neglected solvent-solvent density correlations are indeed important in determining the long-range behavior of the density distribution for this state (see section IV.A). However, an insensitivity of solvation free energies to the long-range solvent-solvent density correlations is consistent with preliminary CC model studies which show that the range over which the density enhancement significantly affects the solvation free energies can be substantially shorter than the range of the density enhancements themselves.67 It follows that the neglected solvent-solvent density correlations, which serve primarily to extend the range of the density enhancement, should have little effect on the solvation free energies. Additional confirmation that such insensitivity is expected has been provided by Chialvo and Cummings,59 who showed that solvation free energies may be expressed solely in terms of the short-range, direct, component of the solute partial molar volume, which is itself expressed in terms of the nondivergent direct correlation functions. It follows that the long-range divergent component of the solute partial molar volume, which reflects the diverging solventsolvent density correlations, does not significantly affect the solvation free energies. These results suggest that the CC model may provide reasonable estimates of solvation energies even for supercritical states very near the critical point, where the
Cl- + CH3Cl Reaction in Supercritical Water
a
b
c
Figure 6. Reaction path-dependent differential solvation free energies for (a) state 1, (b) state 2, and (c) state 3: CC (circle); IC (square); MD (line).
diverging solvent correlation length would make molecular simulations impossible. Lastly, it is important to keep in mind that the excellent agreement found between the CC model and MD simulation is for relatiVe solvation free energies along a reaction path, a situation in which molecular interaction effects neglected in the continuum model tend to cancel. We have not examined the accuracy of the absolute solvation energies predicted by the CC model, but we expect they will demonstrate similar accuracies (and cavity size dependences) as do incompressible continuum models under ambient conditions.
J. Phys. Chem. B, Vol. 101, No. 6, 1997 1069 V. Conclusions Continuum solvation models offer the advantage of rapid computation, making them ideal for studying large or complex systems, as well as for problems requiring a multitude of calculations. An important problem of this type is the study of solute reaction energetics in supercritical fluids, for which solvation energies are required at many points along the solute reaction path for a number of thermodynamic states, including states near the critical point for which computer simulation is infeasible. Yet, conventional continuum solvation models often perform poorly when applied to supercritical fluids. Here, we have demonstrated the usefulness of the compressible continuum solvation model by showing that, compared to molecular dynamics simulation, this model provides an accurate picture of the Cl- + CH3Cl exchange reaction in SCW. In particular, this model was found to provide reasonable estimates of the solvent-solute clustering in these systems. The radial density distributions predicted by the continuum model at three SCW state points of varying bulk compressibility were found to agree rather well with the underlying, nonstructural behavior of those predicted by MD simulations. At the state point closest to the critical point, it was found that the density distribution predicted by the MD simulations appeared to be of longer range than that predicted by the compressible continuum model. This result suggests that at this state point the long-range solvent-solvent density correlationsswhich diverge as the critical point is approachedshave become important, as these effects are included in the MD simulations, but not in the continuum model. The reaction path dependence of the solvation free energies for this reaction are very accurately predicted by the compressible continuum model, with the deviation from the MD predictions being less than 1.5 kcal/mol at all of the reaction path points considered for all three SCW states. This agreement is in direct contrast with the incompressible continuum model predictions, which err by as much as 7 kcal/mol relative to the MD results. Thus, we have demonstrated that the observed failure of these standard models in the supercritical regime is indeed a result of neglected solvent-compression effects and that the compressible continuum model effectively accounts for these effects. Note that the CC model is unique in its ability to isolate the solvent effect of compression on solvation energetics. The fact that this continuum model accurately predicts the relative solvation energies at state point 2, despite the apparent importance of long-range solvent-solvent density correlations at this state, is consistent with earlier suggestions that these longrange, divergent correlations do not significantly affect solvation energies. These results suggest that the compressible continuum model may be used to evaluate solvation energies at supercritical states for which the diverging correlation length prohibits computer simulation. Acknowledgment. We thank Keith Johnston, Peter Rossky, and co-workers for providing us with their simulation data for the chloride plus methyl chloride system in supercritical water. This material is based on work supported by the National Science Foundation under Grant CHE-9307679. S.C.T. gratefully acknowledges an NSF Young Investigator Award and a Camille Dreyfus Teacher-Scholar Award. Appendix. Is the Electrostatic Work a Helmholtz or a Gibbs Free Energy? The relationship between the electrostatic work done on a volume of fluid, F′el, and the Helmholtz, A, and Gibbs, G, free energies has been given by Frank as44,68
1070 J. Phys. Chem. B, Vol. 101, No. 6, 1997
Luo and Tucker
dA ) S dT - P dV + E d(VP) + µ dN
(5)
dG ) S dT + V dP + E d(VP) + µ dN
(6)
and
where the variables S, T, V, µ, and N have their usual meanings. The electrostatic free energy appears in these expressions as
dF′el ) E d(VP)
(7)
where E is the electric field, V is the volume of the fluid in the presence of the field, and P is the polarization of the fluid in the field. For a compressible fluid, the polarization may be changed by the field in two ways. First, there may be an increase in the orientational order of the dipoles, an effect treated herein by the traditional linear response approximation, P(E) ) χE where the susceptibility χ ) 0(s - 1) and s is the solvent’s static dielectric constant (which also incorporates the effect of any electronic polarization). Second, field-induced compression of the fluid will yield an increase in the dipole density. This nonlinear effect on the polarization is treated by introducing a field-dependent solvent dielectric constant, i.e., P(E) ) 0[s(E) - 1]E.38,44 From eqs 5 and 6, it is clear that if the charging is conducted under conditions of constant N, V, and T (as were the Cl- + CH3Cl MD simulations), then dF′el ) dA, whereas if the conditions are constant N, P, and T, then dF′el ) dG. Therefore, whether the free energy provided by the CC model is the Helmholtz or the Gibbs free energy is determined by the conditions under which dF′el ) E d(VP) is computed by this method. To begin, consider the case of small condenser plates immersed in a large reservoir of fluid. As a charge separation is built up on the condenser plates, the volume of fluid in the resulting field between the plates remains fixed. However, fluid may flow freely into this volume from the reservoir which is assumed to be sufficiently large that it acts as a particle bath of constant chemical potential µ. Hence, the effect of the charging process on the constant, specific volume in the field V occurs under conditions of constant µ and T.44 In the solvation problem, we define a volume around the solute cavity over which the solute-induced field is extant. For the continuum calculation, this volume is taken to be comprised of a sum of smaller, specific volume elements, V, each of which is treated analogously to the volume of fluid between the condenser plates; only now the relevant E field is the local value of the soluteinduced field.37 Clearly, the effect of charging the solute on each specific volume element of the fluid, V, can be treated as occurring under conditions of constant µ and T if and only if the solute is effectively at infinite dilution. In order to clearly understand the effect of electrostriction on the thermodynamic properties of the fluid, it is useful to look briefly at the effect of introducing an E field on a specific volume of fluid while requiring the number of molecules in that volume, n, to remain fixed. Note that this means that µ may vary but that the fluid density in the field, F, must remain constant upon charging. As given by Frank,44 the electrostriction pressure required to compensate the effect of the E field on the specific volume element is then
( ) ∂P ∂E2
[( )
∂s 1 ) - 0 F 2 ∂F n,V,T
E,T
]
- (s - 1)
(8)
where the first term on the right-hand side is considered to be a “cohesive” pressure and the second term a “tractive” pressure. Frank showed that the cohesive term is related to the field-
induced alteration of the chemical potential of the fluid, µ, in the specific volume V. That is, in order for there to be no tendency for fluid to leave the specific volume upon charging, the pressure on the fluid outside the specific volume (i.e., outside of the field), Π, would have to decrease according to44
( )
( )
∂s 1 ∂Π ) - 0F 2 2 ∂F ∂E
E,T
(9)
to keep the chemical potential external to the field equal to that of the specific volume in the field. The right-hand side of eq 9 is just the cohesive component of eq 8, demonstrating that the cohesive component indeed represents the effect of the field dependence of the solvent chemical potential. These equations also show that the mechanical pressure difference between the fluid inside and outside of the field, i.e., ∂(P - Π)/∂E2, is just the tractive component of eq 8. Now, reconsider the case of constant chemical potential and temperature. In this case, there will be no change in Π, the pressure on the fluid external to the field, because the fluid density in the field will increase so as to keep the chemical potential constant. When the specific volume is also held constant, this change of density transpires as a result of fluid flow from the reservoir into the specific volume. Frank found this change in density to be given by38,44
[ ( )]
( ) ∂F ∂E2
∂s 1 ) FκT 0F 2 ∂F µ,T
E,T
(10)
where κT is the isothermal compressibility, from which it can be seen that this density change exactly compensates for the cohesive component of the electrostriction pressure. As a result, this term no longer appears in the mechanical pressure required to compensate the E field effect on the specific volume element,44 and thus the mechanical pressure is just
( ) ∂P ∂E2
µ,T
1 ) 0(s - 1) 2
(11)
which is again the tractive term (of eq 8). In the CC solvation model, each specific volume element, V, experiences the solute-induced field under conditions of constant µ and T; these conditions are required by the use of eq 10 in this model. To determine the type of free energy to which the electrostatic work E d(VP) will correspond, we must determine what external conditions these internal constraints correspond to. First, T must be fixed, and second, we restrict our attention to cases in which the total number of solvent molecules in the entire solution, N, is fixed. Now, what about the pressure external to the solution? The local pressure on each specific volume element is dependent upon the local field acting on that element as given by eq 11, and thus, there will be a pressure gradient compensating the electric field. However, beyond the range of the field, this tractive electrostriction pressure will go to zero (eq 11), and the pressure on the fluid outside of the fieldswhich is just Πswill be unchanged by the electrostriction. Hence, upon charging, the pressure external to the solution will be fixed, and the electrostatic free energy as calculated herein may be equated with the Gibbs free energy. What about the external, or total, volume, V? The condition of constant µ was established by assuming infinite dilution, a circumstance under which the amount of solvent remaining outside of any solute-induced field is sufficiently large that it may act as a particle resevoir. This implies that the number of solvent molecules outside of the solute field is effectively constant, and since the pressure outside of the field is also
Cl- + CH3Cl Reaction in Supercritical Water constant, it follows that the volume of solvent outside the field is effectively constant. Since the sum of the specific volume elements in the field is also constant, the total volume must be effectively constant. Therefore, the electrostatic work computed from this model may be equated with the Helmholtz free energy for infinitely dilute systems. Note, however, that to model a constant volume system of insufficient size or dilution would require the use of a variable-µ analog of eq 10 for evaluating F(E). References and Notes (1) Savage, P. E.; Gopalan, S.; Mizan, T. I.; Martino, C. J.; Brock, E. E. AIChE J. 1995, 41, 1723. (2) Subramaniam, B.; McHugh, M. A. Ind. Eng. Chem. Res. 1986, 25, 1. (3) Ikushima, Y.; Saito, N.; Yokoyama, T.; Hatakeda, K.; Ito, S.; Arai, M.; Blanch, H. W. Chem. Lett. 1993, 109. (4) Johnston, K. P.; Haynes, C. AIChE J. 1987, 33, 2017. (5) Peck, D. G.; Mehta, A. J.; Johnston, K. P. J. Phys. Chem. 1989, 93, 4297. (6) Hrnjez, B. J.; Mehta, A. J.; Fox, M. A.; Johnston, K. P. J. Am. Chem. Soc. 1989, 111, 2662. (7) Antal Jr., M. J.; Brittain, A.; DeAleida, C.; Ramayya, S.; Roy, J. C. In Supercritical Fluids; Squires, T. G., Paulaitis, M. E., Eds.; ACS Symposium Series 329; American Chemical Society: Washington, DC, 1987. (8) Narayan, R.; Antal, Jr., M. J. J. Am. Chem. Soc. 1990, 112, 1927. (9) Paulaitis, M. E.; Alexander, G. C. Pure Appl. Chem. 1987, 59, 61. (10) Wu, B. C.; Klein, M. T.; Sandler, S. I. Ind. Eng. Chem. Res. 1991, 30, 822. (11) Klein, M. T.; Mentha, Y. G.; Torry, L. A. Ind. Eng. Chem. Res. 1992, 31, 182. (12) Tucker, S. C.; Gibbons, E. M. In Structure and ReactiVity in Aqueous Solution; Truhlar, D. G., Cramer, C. J., Eds.; ACS Symposium Series; American Chemical Society: Washington, DC, 1994. (13) Haar, L.; Gallagher, J. S.; Kell, G. S. NBS/NRC Steam Tables, Hemisphere: Washington, DC 1984. (14) Moore, J. W.; Pearson, R. G. Kinetics and Mechanism; Wiley: New York, 1981. (15) Chandler, D. J. Chem. Phys. 1978, 68, 2959. (16) Tucker, S. C. In New Trends in Kramers Reaction Rate Theory; Understanding Chemical ReactiVity, Talkner, P.; Ha¨nggi, P., Eds.; Kluwer Academic: Dordrecht, 1995. (17) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University: Oxford, 1989. (18) Chandrasekhar, J.; Smith, S. F.; Jorgensen, W. L. J. Am. Chem. Soc. 1985, 107, 154. (19) Jorgensen, W. L. Acc. Chem. Res. 1989, 22, 184. (20) Kollman, P. A.; Jr., K. M. M. Acc. Chem. Res. 1990, 23, 246. (21) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids: Academic Press: London, 1986. (22) Mouritsen, O. G. Computer Studies of Phase Transitions and Critical Phenomena; Springer: Berlin, 1984. (23) Martinez, H. L.; Ravi, R.; Tucker, S. C. J. Chem. Phys. 1996, 104, 1067. (24) Davis, M. E.; McCammon, J. A. Chem. ReV. 1990, 90, 509. (25) Cramer, C. J.; Truhlar, D. G. In ReViews in Computational Chemistry; Boyd, D. B., Lipkowitz, K. B., Eds.; VCH: New York, 1995. (26) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 1978. (27) Tomasi, J.; Persico, M. Chem. ReV. 1994, 94, 2027. (28) Gilson, M.; Sharp, K. A.; Honig, B. J. Comput. Chem. 1988, 9, 327. (29) Rick, S. W.; Berne, B. J. J. Am. Chem. Soc. 1994, 116, 3949.
J. Phys. Chem. B, Vol. 101, No. 6, 1997 1071 (30) Hyun, J.-K.; Babu, C. S.; Ichiye, T. J. Phys. Chem. 1995, 99, 5187. (31) Cossi, M.; Mennucci, B.; Tomasi, J. Chem. Phys. Lett. 1994, 228, 165. (32) Ehrensen, S. J. Phys. Chem. 1987, 91, 1868. (33) Abraham, M. H.; Liszi, J.; Meszaros, L. J. Chem. Phys. 1979, 70, 2491. (34) Beveridge, D. L.; Schnuelle, G. W. J. Phys. Chem. 1975, 79, 2562. (35) Johnston, K. P.; Bennett, G. E.; Balbuena, P. B.; Rossky, P. J. J. Am. Chem. Soc. 1996, 118, 6746. (36) Luo, H.; Tucker, S. C. J. Am. Chem. Soc. 1995, 45, 11359. (37) Luo, H.; Tucker, S. C. J. Phys. Chem. 1996, 100, 11165. (38) Wood, R. H.; Quint, J. R.; Grolier, J.-P. E. J. Phys. Chem. 1981, 85, 3944. (39) Klein, M. T.; Torry, L. A.; Wu, B. C.; Townsend, S. H.; Paspek, S. C. J. Supercrit. Fluids 1990, 3, 222. (40) A fourth pressure was not considered because the change in the rate with pressure was too small to yield a conclusive comparison with theory. (41) Bennett, G. E.; Rossky, P. J.; Johnston, K. J. Phys. Chem. 1995, 99, 16136. (42) Flanagin, L. W.; Balbuena, P. B.; Johnston, K. P.; Rossky, P. J. J. Phys. Chem. 1995, 99, 5196. (43) Balbuena, P. B.; Johnston, K. P.; Rossky, P. J. J. Phys. Chem. 1995, 99, 5196; J. Am. Chem. Soc. 1994, 116, 2689. (44) Frank, H. S. J. Chem. Phys. 1955, 23, 2023. (45) Wood, R. H.; Carter, R. W.; Quint, J. R.; Majer, V.; Thompson, P. T.; Boccio, J. O. J. Chem. Thermodyn. 1994, 26, 225. (46) Klapper, I.; Hagstrom, R.; Fine, R.; Sharp, K.; Honig, B. Proteins 1986, 1, 47. (47) Wood, R. H.; Quint, J. R. J. Phys. Chem. 1989, 93, 936. (48) Hill, P. G. J. Chem. Phys. Ref. Data 1990, 19, 1233. (49) Archer, D. G.; Wang, P. J. Phys. Chem. Ref. Data 1990, 19, 371. (50) Chandrasekhar, J.; Smith, S. F.; Jorgensen, W. L. J. Am. Chem. Soc. 1984, 106, 3049. (51) Huston, S. E.; Rossky, P. J.; Zichi, D. A. J. Am. Chem. Soc. 1989, 111, 5680. (52) Jorgensen, W. L.; Madura, J. D.; Swensen, C. J. J. Am. Chem. Soc. 1984, 106, 6639. (53) Jean-Charles, A.; Nicholls, A.; Sharp, K.; Honig, B.; Tempczyk, A.; Hendrickson, T. F.; Still, W. C. J. Am. Chem. Soc. 1991, 113, 1454. (54) The reaction path dependence of the R′ set was taken to be the same the R set. (55) Berendsen, H. J. C.; Potsma, J. P. M.; von Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, C., Ed.; Reidel: Dordrecht, 1981. (56) de Pablo, J. J.; Prausnitz, J. M.; Strauch, H. J.; Cummings, P. T. J. Chem. Phys. 1990, 93, 7355. (57) Kaski, K.; Binder, K.; Gunton, J. D. Phys. ReV. B 1984, 29, 3996. (58) Debenedetti, P. G. Chem. Eng. Sci. 1987, 42, 2203. (59) Chialvo, A. A.; Cummings, P. T. AIChE J. 1994, 40, 1558. (60) Carlier, C.; Randolph, T. W. AIChE J. 1993, 39, 876. (61) Cummings, P. T.; Cochran, H. D.; Simonson, J. M.; Mesmer, R. E.; Karaborni, S. J. Chem. Phys. 1991, 94, 5606. (62) Sun, Y. P.; Fox, M. A.; Johnston, K. P. J. Am. Chem. Soc. 1992, 114, 1187. (63) Petsche, I. B.; Debenedetti, P. G. J. Chem. Phys. 1989, 91, 7075. (64) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (65) Knutson, B. L.; Tomasko, D. L.; Eckert, C. A. Debenedetti, P. G.; Chialvo, A. A. In ACS Symposium Series 488; Bright, F. V., McNally, M. E., Eds.; American Chemical Society: Washington, DC, 1992. (66) The IC results are as recalculated by us. Deviations from those presented in ref 41 falls within the numerical error of our finite difference calculation. (67) Luo, H.; Tucker, S. C. Work in progress. (68) We give all equations in SI units; 0 is the permittivity of free space.
