J. Phys. Chem. B 2001, 105, 6611-6619
6611
Stability of Ionic and Radical Molecular Dissociation Pathways for Reaction in Supercritical Water† Robin E. Westacott,‡,§ Keith P. Johnston,*,‡ and Peter J. Rossky*,§ Department of Chemical Engineering and Institute for Theoretical Chemistry, Department of Chemistry and Biochemistry, The UniVersity of Texas at Austin, Austin, Texas 78712 ReceiVed: January 3, 2001; In Final Form: March 6, 2001
Molecular dynamics simulations have been used to examine the effect of supercritical water solvent density on the competition between reaction pathways for the dissociation step of a model SN1 reaction. The effects are investigated using an empirical valence bond theory that explicitly includes the effects of solvation, particularly those on the diabatic ionic state. At low supercritical water densities, the solvent stabilization is insufficient to give rise to a local minimum on the free energy surface corresponding to a contact ion pair intermediate, although the free energy surface is completely ionic in character to solvent densities less than 0.05 g cm-3. The nature of the surface is also changed by solvent density; the change from a mostly covalent (80%) molecule to completely ionic dissociation products is decreasingly rapid as supercritical water density is decreased. Radial density functions reflect how solvation changes along the reaction coordinate and how local density enhancement provides the solvation required to stabilize the ionic products. These calculations indicate that the diabatic ionic state is lowest in free energy until extremely low supercritical water solvent density (∼0.03 g cm-3), considerably lower than would be expected if local density enhancement were ignored, as in a simple Born model calculation. The free energy difference between the two pure states at the dissociation plateau indicates that covalent products may be expected to reach approximately 22% of the total at the lowest density (0.0435 g cm-3) considered here.
1. Introduction Recently, motivated by industrial interest, supercritical fluids have received a great deal of attention as solvents in a variety of processes ranging from oxidation of organic wastes to treatment of radioactive materials. The base for current understanding of supercritical water was set over the last three decades and includes experimental determination of calorimetric1 and electrochemical2,3 properties. More recently, spectroscopic methods have been used to investigate structural properties of supercritical water.4,5 Up-to-date reviews of the understandings of supercritical fluids can be found in a recent issue of Chemical ReViews.6 The nature of supercritical water is such that, above the critical temperature, the density can be varied at constant temperature over an enormous range from ambient-like to vapor density using only small changes in the pressure. This means that the solvating power afforded by the solvent can be changed from highly polar to that of an organic solvent to that of a gas. As the solvent density decreases, so do the bulk dielectric constant and, therefore, the stabilization afforded to ions and dipolar species. At the critical point for water (Tc ) 647.15 K, Pc ) 22.1 MPa, Fc ) 0.322 g cm-3), the dielectric constant is reduced to ) 5.3 (cf. dimethyl ether ) 5.02, water ) 80, both at 298 K and atmospheric pressure). For the development of industrial processes involving supercritical fluids, it is essential to understand the effects that solvent †
Part of the special issue “Bruce Berne Festschrift”. * Corresponding authors. ‡ Department of Chemical Engineering. § Institute for Theoretical Chemistry, Department of Chemistry and Biochemistry.
density has on the progress of chemical reactions. The experimental study of solvation of ions and molecules in supercritical fluids is well documented, and studies of chemical reactions in these solvents have recently been performed.7-10 In many cases, the Born model11 has been used to determine the free energy of solvation in supercritical water. But it has become clear that since there are potentially large fluctuations in the local solvent density, a simple continuum description of the supercritical fluid is not sufficient.12 A substantial body of work has been published on the density inhomogeneities in supercritical fluids and the nature of density enhancements (known as electrostriction) around solute species.13,14 Several workers have described more complex inhomogeneous continuum models that take into account the local solvent density around solutes in supercritical fluids.15 However, the most complete treatment of the solvent involves explicit inclusion of molecular solvent effects. As the density of supercritical water is decreased, it behaves less like a liquid and increasingly takes on the properties of a hightemperature gas. It is thus important to understand the effects that changes in the properties of the supercritical water have on the mechanism of chemical reactions performed in this fluid, especially those involving charge transfer, which requires solvent stabilization of ionic species. The use of supercritical water as a reaction medium has received a great deal of attention over the past few years;16 for a further example, see ref 17. Hydrolysis reactions, in which the water solvent takes part in the reaction, have received particular attention, since they can be used for the treatment of organic wastes. Notable studies involve the hydrolysis of methylene chloride10 and ethers8,9 and the oxidation of methane18 and methanol19 in supercritical water. Several workers have indicated that reactions in supercritical water solvent may
10.1021/jp010005y CCC: $20.00 © 2001 American Chemical Society Published on Web 04/11/2001
6612 J. Phys. Chem. B, Vol. 105, No. 28, 2001 progress through competing reaction paths. Antal et al.7 showed that reactions involving simple organic molecules in supercritical water progress through different mechanisms, which correlate with the ionic product, Kw, of the supercritical solvent. The ionic product is temperature- and density-dependent and can be as low as 10-24 at low density. They indicated that when Kw > 10-14 reactions proceed through heterolytic mechanisms and when Kw < 10-14, the reaction proceeds through homolytic mechanisms. Townsend et al.8 have shown that ether hydrolysis progresses in competition with pyrrolysis in supercritical water. They investigated a series of benzyl and phenyl ethers and concluded that the selectivity for the hydrolysis pathway increases with supercritical water solvent density. Their results indicate a density at which there is a crossover in preferred reaction pathways. This crossover supercritical water solvent density is dictated by the properties of the reactant/solute species and occurs in the range Fw ) 0.09-0.20 g cm-3. In another study by the same group, Huppert et al.9 studied the hydrolysis of o-methoxyphenol in supercritical water. They showed that the hydrolysis rate constant increased with increasing solvent polarity and that the propensity for products of hydrolysis reactions increased with solvent density. It is therefore clear that changes in supercritical water solvent density can have a profound effect on the competition between reaction pathways. These effects should be correlated with the properties of the solvent, but it is likely that the nature of the solute plays a substantial role, especially where solvent-solute interactions are significant.20 There is considerable interest in reactions in supercritical water that follow ionic heterolytic reaction pathways rather than homolytic pathways, as the latter tend to be less selective. In the work described in this paper, we will use computer simulation of a simple model reaction having competing pathways to investigate supercritical water solvent effects. A preliminary report of this work has been published elsewhere.21 The use of computer simulation methods to investigate chemical reactions is widespread, and there are already several studies of chemical reactions in supercritical fluids. Balbuena et al.22 have used molecular dynamics simulations to investigate the SN2 reaction of Cl- and CH3Cl in supercritical water. They found that even at very low densities, the barrier to reaction was intermediate between the gas phase and ambient water limits, although closer to the ambient water case. It was also observed that the hydration numbers around the reactants were similar to that seen for ambient water, but there was less hydrogen bonding largely due to the increase in temperature. They also showed that quite extreme conditions (i.e., greater than 760 K and density lower than 0.08 g cm-3) are required to remove half of the water molecules from the first solvation shell of the chloride ion.12 Rate constants for the SN2 reaction 9-12 orders of magnitude greater than for ambient water were also estimated, largely due to temperature effects. Marrone et al.23 used ab initio methods with Kirkwood theory to describe solvation effects on the hydrolysis of methylene chloride in supercritical water. These calculations showed that the rate of hydrolysis decreases by about 3 orders of magnitude when the solvent is heated from subcritical to supercritical temperatures. They also showed that the separation distance between the reacting species becomes smaller as the dielectric constant, and therefore density, of the solvent decreases. Hydrolysis reactions are an important class of laboratory reaction that are, at the same time, ubiquitous in such industrial processes as the production of alcohols and ethers and the destruction of waste organic materials. Here we will investigate
Westacott et al. the prototypical hydrolysis of tert-butyl chloride. For steric reasons, hydrolysis in this system follows the SN1 mechanism, and it is one of the main model systems used in the development of methods directed toward understanding aspects of physical organic chemistry and solvation. It is generally rationalized by an ion pair mechanism in which the reactant, RX, ionizes to form a stable intermediate contact pair R+X-, followed by a solvent separated pair R+//X-, and then isolated ions R+ + X-.24 Any of these products may then react with a nucleophile to produce substitution products. The ionization is conventionally the rate-determining step for these types of reactions and can be used to demonstrate the solvating power of a solvent. In particular, the ion products are stabilized by polar solvents. In the gas phase, the lack of solvation leads to dissociation with the formation of radicals. Thus, the degree to which the solvent can stabilize the ions dictates the progress of the reaction. The strong solvent dependence of this reaction is illustrated by the fact that the reaction rate decreases by a factor of 90 if the solvent is changed from water to a 50:50 v/v mixture of ethanol and water and by a factor of 1.2 × 106 if the solvent is changed to pure ethanol.25 The rate constant increases by a factor of 6.6 × 1017 from the gas phase to ambient water.25 For supercritical water, for which the density and, therefore, dielectric constant vary greatly with only small changes in pressure, the environment of reacting species also changes dramatically. It is worth noting here that the complementary elimination reaction (E1) also requires stabilization of ionic intermediates to proceed but requires a basic environment for removal of a proton from the cation that leads to the formation of an alkene. The aqueous tert-butyl chloride system under ambient conditions has been the subject of several detailed theoretical studies, including those by Jorgensen et al.,26 Ford and Wang,27 Hynes et al.,28,29,30 Keirstead et al.,31 Hartsough and Merz,32 and Takahashi et al.33 The potential of mean force (PMF) for the association of the tert-butyl cation and chloride anion was determined by Jorgensen et al.26 using free energy perturbation and Monte Carlo simulations. They concluded that the solvent separated ion pair (t-Bu+//Cl-) is more stable than the contact ion pair (t-Bu+Cl-) by about 2 kcal mol-1. They also found that for t-Bu+Cl-, a water molecule approached very close to the back of the tert-butyl cation, implying that the solvent did take the role of nucleophile in these reactions. Keirstead et al.31 used an empirical valence bond treatment that included solvent effects to study the formation of the contact ion pair in liquid water. The PMF was characterized by a free energy barrier to ion formation of about 20 kcal mol-1, which is in good agreement with that determined by Abrahams by experiment and classical thermodynamics.34 Further, several workers have used semiempirical or ab initio quantum chemical methods to investigate this system. Ford and Wang27 have used the AM1 and MNDO semiempirical methods with reaction field to calculate solvation energies for molecules, cations, anions, and reacting systems, including t-BuCl. The PMF generated using the AM1 method had a reaction barrier of 20.5 kcal mol-1 at 3.0 Å and a barrier to recombination of t-Bu+Cl- of 0.9 kcal mol-1. Hartsough and Merz32 used the quantum mechanics/ molecular mechanics (QM/MM) method using the PM3 Hamiltonian for the tert-butyl chloride solute, which was surrounded by 483 water molecules. They obtained a deep minimum at the equilibrium bond length and a barrier of 17.8 kcal mol-1 for dissociation to t-Bu+Cl-. The contact ion pair was seen to be stabilized in a well of depth 8.7 kcal mol-1, with respect to separated ions, and a shallow minimum corresponding to the solvent separated ion pair was observed at 6.4 Å. Takahashi et
Ionic and Radical Molecular Dissociation Pathways al.33 used the MP2 level of theory and solvent reaction field to calculate the PMF for tert-butyl chloride dissociation. They obtained an activation energy of ∼15 kcal mol-1 and a very shallow minimum for solvent separated ion pair but no minimum for the contact ion pair. In the present work, we have chosen to use a valence-bondbased model, with an explicit solvent analogous to that of ref 31. Such models are particularly well suited to the theoretical study of charge-transfer reactions. Polanyi and co-workers35,36 used valence bond theory to describe charge transfer as arising from a mixture of two diabatic electronic states. This picture was developed in more detail by Pross and Shaik37 and introduced into simulation in the form used here by Warshel et al.38 This approach used the valence bond theory to describe the electronic ground state of a molecule as a linear combination of diabatic states. For the tert-butyl chloride system, the two most important states are a covalent state with a bonding interaction and no charge separation and an ionic state with full charge separation as tert-butyl cation and chloride anion. The weight of each form, cc and ci, respectively, in the most stable state (adiabatic ground state) is determined by solving the secular equation for the model. For tert-butyl chloride in the gas phase, the dipolar molecule (ground state) is predominantly covalent but has substantial ionic character that gives rise to a dipole moment of 2.13 D (implying ci2 ) 0.20). At larger separations, the ground state of the system is completely covalent (cc2 ) 1.0), as radicals form when the carbon-chlorine bond breaks. When the tert-butyl chloride molecule is dissolved in a polar solvent, the solvent stabilizes the ionic state, and the corresponding diabatic energy surface is lowered in energy. With sufficient solvation, the ionic state lies below the covalent surface at large separations and potentially quite small separations, as well. This leads to the formation of ions when the molecule dissociates. In supercritical water, where density and corresponding dielectric constant can be changed dramatically and continuously above the critical temperature, the solvent stabilization afforded to ions is variable. As seen in the experimental work described above, this leads to the possibility of competition between reaction mechanisms, or, put another way, changes in the reaction free energy surface that make alternate reaction pathways accessible. For the simple case of tert-butyl chloride, the competing reaction mechanisms are ionic and radical dissociation. Here, we will employ the empirical valence bond (EVB) approach due to Warshel38 that includes solvent effects, previously used by Hynes and co-workers28,29,30,31 to investigate tertbutyl halides. Molecular dynamics simulations with free energy perturbation39 have been performed to obtain the adiabatic reaction free energy surface, or potential of mean force (PMF), along the reaction coordinate as a function of supercritical water density and for ambient water. This is the condensed phase equivalent of the gas-phase free energy surface. For supercritical water, we expect variation between the limiting gas phase and ambient water cases, described above. The remainder of this paper is organized as follows. Section 2 contains a complete description of the models, potentials, and methodology used in this work. It describes the valence bond model used to describe the tert-butyl chloride system and the potentials used to describe the two pure diabatic states. The molecular dynamics simulations are also described in detail. Section 3 contains the results of this study, which are presented in terms of thermodynamics and structural features and a discussion of their implications. Finally, conclusions are drawn, and indications of future directions are given in Section 4.
J. Phys. Chem. B, Vol. 105, No. 28, 2001 6613 2. Methodology Hamiltonian. The form of the adiabatic Hamiltonian described by Keirstead et al.31 has been used to describe the tertbutyl chloride in water system. This arises through a two state EVB treatment of the tert-butyl chloride molecule, including solvation effects. The ionic and covalent states are considered diabatic and electronically coupled via the coupling element Vel. Solving the secular determinant with neglect of off-diagonal overlap gives rise to an expression for the energies of the ground and first excited states
1 1 2 E0,1 ) (Hcov + Hion) - (∆H2 + 4Vel )1/2 2 2
(1)
and expressions for the weights of each pure state
ccov ) 2Vel2{4Vel + [∆H + (∆H2 + 4Vel )1/2]2}-1/2 (2) 2
ccov ) (1 - ccov2)1/2
2
(3)
where Hcov and Hion are the Hamiltonians for the pure diabatic covalent and ionic states, respectively, ∆H ) Hcov - Hion, and Vel is the electronic coupling between the two pure states (discussed later). The classical adiabatic Hamiltonian is defined from the ground-state energy E0 in eq 1 using the negative root and is explicitly given by
1 1 Had ) Kw + Ks + Vw + (Vis + Vcs ) + (Viinter + Vcinter) 2 2 1 {[(Vis - Vcs ) + (Viinter - Vcinter)]2 + 4Vel2}1/2 (4) 2 where Ks and Kw are the kinetic energies of the solute and the water, respectively, Vw is the potential energy of the water, Vis and Vcs are the potential energies of the solute in the ionic and covalent states, respectively, and Viinter and Vcinter are the watersolute potential energies for the ionic and covalent states, respectively. Potentials and Parameters. The potentials used in this study have been substantially modified from those used by Keirstead et al.31 in order to answer the questions of interest in this work. In particular, since a wide range of solvent environments is of interest, it is important that the model used here reproduces the gas-phase dipole moment of the tert-butyl chloride system over the full range of separations. In the following subsections, potentials for the pure ionic state, pure covalent state, and water are presented, and the differences from those used by Keirstead et al.31 are discussed. Further, the electronic coupling between the two pure states has been assigned an exponential form in line with expectations and the work of Hynes et al.28-30 Solute Potential. The solute intramolecular interactions have been considered using a five-site model for tert-butyl chloride.26,31 Here, the tert-butyl group is represented by three united atom methyl groups trigonally arranged around the central carbon and on a plane perpendicular to the carbon-chlorine bond. The tert-butyl group was considered to have this fixed geometry throughout the simulations, and the parameters described below provide the equilibrium molecule having this planar tert-butyl arrangement with the same properties as the real tetrahedral molecule. The potential models used to describe the pure ionic and covalent states have been modified from those used by Keirstead et al.31 in order to mimic the gas-phase dipole moment of the tert-butyl chloride molecule within the EVB model. The pure ionic state is described by a combination of a Buckingham (exp-6) potential and a Coulombic term for
6614 J. Phys. Chem. B, Vol. 105, No. 28, 2001
Westacott et al.
TABLE 1: Parameters for Ionic Solute Intramolecular Interactions atom
q (e-)
σ (Å)
(kcal mol-1)
C Me Cl
0.4 0.2 -1.0
2.250 3.000 4.417
0.050 0.100 0.118
TABLE 2: Parameters for Solute-Solvent Interactions
interactions of the central carbon and methyl groups of the tertbutyl ion with the chloride ion
Vi(r) ) Voffset + ACCle-RCClrCCl + 4CCl
( ) ( )
3AMeCle-RMeClrMeCl + 12MeCl
qCqCl σCCl + + rCCl rCCl qMeqCl σMeCl 6 +3 (5) rMeCl rMeCl 6
where Me refers to the methyl groups in the tert-butyl species, C to the central carbon in the tert-butyl species, and Cl to the chlorine. Voffset is the difference in potential energy between the ionic and covalent states at r ) ∞ and is equal to the sum of the first ionization energy of tert-butyl and the electron affinity of chlorine.31 For this work, we have used the same value as that of Keirstead et al.,31 where Voffset ) 71.19 kcal mol-1. Here, the tert-butyl group is taken to be planar for all C-Cl separations, and therefore
rMeCl ) (rCCl2 + b2)1/2
(6)
where b is the carbon to methyl group bond length. The attractive r-6 term in the Buckingham potential and the Coulombic term are identical to that used by Keirstead et al.,31 which was taken from the work of Jorgensen et al.,26 with the same parameters being used (Table 1). The exponential terms were used to provide a more justifiable representation of the repulsive interactions. This was found to impart more ionic character into the tert-butyl chloride molecule in the gas phase by reducing the gradient of the repulsive wall near the equilibrium bond length. The parameters for the repulsive part of this potential are ACCl ) AMeCl ) 24 100 kcal mol-1 and RCCl-1 ) RMeCl-1 ) 0.531 Å-1. For the pure covalent state, a Morse potential was used
Vc(r) ) V0(e-2(r-r0)/a - 2e-(r-r0)/a)
(7)
where the separation is based on the distance between the central carbon of the tert-butyl group and the chlorine atom. The values for r0 (1.83 Å) and V0 (60.72 kcal mol-1) were adjusted from those values used by Keirstead et al.31 In effect, reduction of V0 from the value of Keirstead et al. acknowledges that the carbon to chlorine bond in the mixed state will have some ionic character, and this serves to increase the gas-phase dipole moment of the molecule. Correspondingly, r0 was increased from the equilibrium bond length of the molecule so that the bond length for the adiabatic ground-state remains at the experimental value. The value of the a parameter, which describes the curvature of the exponential function, was calculated from (2V0/4π2µc2ω2)1/2, where ω ) 560 cm-1 is the experimental bond vibration frequency. Here, we have assumed that since the tert-butyl chloride molecule is ∼80% covalent, the curvature of the ground state at the global minimum is equivalent to the curvature of the pure covalent state. The electronic coupling between the two pure diabatic states has been investigated by Hynes et al. and argued to be exponential-like.28,29,30 This is consistent with electronic overlap
atom
q (e-)
σ (Å)
(kcal mol-1)
Cl Clt-Bu t-Bu+ O H
0 -1 0 +1 -0.8476 +0.4238
3.35 4.417 5.3304 5.3304 3.1655 -
0.3448 0.118 0.1695 0.1695 0.6501696 -
effects. In this work, an exponential functional form has been used to describe the change in electronic coupling with separation
Vel(r) ) Ae-βr
(8)
The form of the exponential was chosen to reproduce two predetermined points: the inferred value of Vel (15 kcal mol-1)31 at the crossing point of the diabatic pure states31 and an acceptable value of the gas-phase dipole moment while leaving the position and depth of the adiabatic minimum at the experimental values. We have used A ) 4765 kcal mol-1 and β ) 1.767 Å-1. The calculated dipole moment for this model is 1.7 D (cf. 2.13 D experimental). SolVent Potential. Water was described using the SPC/E potential,40 which has previously been shown to give a good representation of bulk water under supercritical41 conditions. The critical properties of SPC/E water have been determined by Guissani and Guillot42 and were found to be 640 K, 0.29 g cm-3, and 16 MPa. This model describes the interactions between water molecules using one Lennard-Jones site positioned on the oxygen atom and a point charge located on each of the three atoms. The parameters for SPC/E water are presented in Table 2. SolVent-Solute Potential. The solvent interacts with a twosite solute species, where the tert-butyl group is treated as a united atom, following Keirstead et al.31 The interactions between the water molecules and the tert-butyl chloride molecule were described using electrostatic and Lennard-Jones potentials as 3Nw 2
Vinter )
∑ ∑ i)1 χ)1
{ [( ) ( ) ] } σiχ
4iχ
riχ
12
-
σiχ riχ
6
+
qiqiχ riχ
(9)
The Lennard-Jones parameters were derived as the geometric average of the individual components. The coupled parameters are presented in Table 2. Simulations. The reaction coordinate was defined as the distance, rA, between the central carbon atom of the tert-butyl group and the chlorine atom. Molecular dynamics (MD) simulations with statistical mechanical perturbation theory39 were used to calculate the change in Helmholtz free energy in the solution at each step δrA along the reaction coordinate. The change in Helmholtz free energy was obtained from the ensemble average at fixed rA, when the system was perturbed from rA to rA + δrA according to the exact expression
δA(rA) ) -kBT ln〈exp[-(Ug+s(rA + δrA) Ug+s(rA))/kBT]〉rA (10) where Ug+s is the sum of the water-solute and solute-solute potential energies and 〈.....〉 indicates a thermal average taken for the solution at fixed rA when the system is perturbed from rA to rA + δrA. MD was performed at sequential intervals of 0.25 Å along the reaction coordinate. At each distance, rA, 20 ps of MD was
Ionic and Radical Molecular Dissociation Pathways
J. Phys. Chem. B, Vol. 105, No. 28, 2001 6615
performed after equilibration, in which the average required in eq 10 was evaluated for δrA ) (0.125 Å. In all cases, the time step was 2 fs. Cubic periodic boundary conditions were applied, and the Ewald sum was used to calculate the long-range forces. A total of 500 water molecules was used, together with one solute molecule. The size of the simulation box was chosen to be appropriate for the required water density. Simulations were performed at the 673 K and at solvent densities of 0.29, 0.087 and 0.0435 g cm-3. The calculations were also performed for ambient water (298 K and 1 g cm-3) for comparison. The PMF was calculated as r
W(r) )
∑ δA(rA) r)∞
(11)
where W(r) is the change in Helmholtz energy of the system at rA, which results from the sum of all the contributions δA(rA) when the tert-butyl and chloride moieties are brought together in small steps δrA from infinite separation to the distance rA. The EVB method captures the effect of the solvent on the relative energies of the pure diabatic states. It is expected that the solvent will have a relatively negligible effect on the nonpolar pure covalent state, which will therefore be similar to the gas phase, but that the effects on the ionic state will be dramatic. Thus, we set the relative free energies of the two states
{
0 ; neutral radicals W(|rA| ) ∞) ) V offset + ∆Asolv ; ions
(12)
where ∆Asolv is the Helmholtz solvation free energy of the isolated ions, and
Voffset ) Eion + Eaff
(13)
where Eion and Eaff are the ionization energy of tert-butyl and the electron affinity of chlorine, respectively. The values of ∆Asolv for the chloride anion have been taken from the literature,12 except for the lowest supercritical water density considered here. The latter and all values for the tert-butyl cation have been calculated using the Born model for solvation of ions following the method described in ref 12. 3. Results and Discussion In this work, the reaction free energy surface along the reaction coordinate for the dissociation of tert-butyl chloride has been calculated for the gas phase, in ambient water and three supercritical water densities. The effect of solvation is expected to have a dramatic effect on the character of these surfaces. A key property of a solvent that characterizes its solvating power is the dielectric constant. For SPC/E water, the dependence of the dielectric constant on temperature and density can be calculated from43
(F) ) 1 +
β(F) C(F) + 2 T T
B(F) ) b1F + b2F2 + b3F3 C(F) ) c1F + c2F2 + c3F3
(14)
where b1 ) 3690.27, b2 ) 5346.95, b3 ) -1269.30, c1 ) 2.39747E6, c2 ) 4.02446E6, and c3 ) -923962. For the temperature (673 K) and solvent densities used in this work, the dielectric constants are 5.44, 2.06, and 1.5, respectively. The ambient value is 80.
Figure 1. Gas-phase reaction energy surface for the tert-butyl chloride system. The line styles represent the adiabatic ground state (solid line), the pure covalent state (dashed line), the pure ionic state (dotted line), and electronic coupling (dot-dashed line).
The PMF as a function of reaction coordinate for a range of solvent densities has previously been reported by us.21 In the gas phase, the adiabatic ground-state energy is very similar to that published by Keirstead et al.31 and has a single minimum corresponding to the equilibrium bond distance and a depth corresponding to the bond dissociation energy (see Figure 1). As expected, the main differences from this earlier work occur for the pure diabatic states, due to the different potential forms used. At small separation distances between the tert-butyl and chlorine, the covalent state is lowest in energy. At intermediate distances, the two surfaces cross twice, and between these points, the ionic state is lowest in energy. At larger separations, the covalent state again has lowest energy, indicating that radicals form when the carbon-chlorine bond breaks. Both crossing points of the diabatic surfaces are avoided by the adiabatic surface. Since the electronic coupling decreases exponentially to zero with increasing separation in the present model, the covalent state is completely dominant after the second crossing point of the diabatic curves, indicated by the fact that the adiabatic curve and the covalent curve exactly coincide at large separations. At small separations, the adiabatic PMF is a combination of both pure diabatic states (∼20% ionic). At intermediate separations, the degree of ionic character increases until it is almost 100%, but at larger separations, it decreases rapidly to zero (Figure 2). This suggests that as the carbonchlorine bond is stretched, the electronegativity of chlorine pulls the electron pair in the bond onto itself, but as the bond starts to break, one of the electrons is transferred back to the tertbutyl group since there is insufficient stabilization (i.e. no polar solvent) of ions. The PMF in ambient water is shown in Figure 3. Compared to that for the gas phase (Figure 1), the pure covalent curve is essentially identical on this scale. The main differences are due to the effect of the solvent on the pure ionic curve; solvation leads to a stabilization by ∼80 kcal mol-1 at the minimum of the pure ionic state. Consequently, the adiabatic ground state is stabilized by ∼14 kcal mol-1 at its minimum compared to the gas phase. The adiabatic PMF has a minimum at r )1.8 Å and W(r) ) -98 kcal mol-1. The height of the barrier at 2.35 Å is 23 kcal mol-1. The activation energy is in good agreement with the experimental values obtained by Winstein and Fainberg,44 who obtained the value of 19.5 kcal mol-1 and the values obtained through simulation methods, ∼20,31 ,32 and 18.1 kcal mol-1.27 The position of the transition state is in excellent
6616 J. Phys. Chem. B, Vol. 105, No. 28, 2001
Figure 2. Fraction of Ionic character for the tert-butyl chloride system in the gas phase.
Figure 3. Reaction free energy surface for tert-butyl chloride in ambient water. The key is the same as that for Figure 1, with electronic coupling omitted.
agreement with previous simulation results, 2.28 Å31 and 2.3 Å,32 although much smaller than the values obtained by Okuno (3.075 Å)45 and Takahashi (3.0 Å).33 There is a shallow, broad minimum corresponding to the contact ion pair at ∼3.1 Å. The position of this minimum is in excellent agreement with that found by Keirstead et al.,31 Hartsough and Merz,32 Jorgensen et al.,26 and Ford and Wang.27 The contact ion pair is stabilized by 7 kcal mol-1 relative to the transition state and therefore lies 16 kcal mol-1 above the undissociated molecule. The depth (7 kcal mol-1) is in reasonable agreement with these earlier simulation results, and it is in excellent agreement with the experimental estimate of Abraham.34 The shape of the PMF and the energy at large separations are dramatically different from the gas phase. The steep repulsive wall of the ionic surface at small separation is affected very little by solvation, but the surface at large separations is substantially lowered in energy, leading to the broad, shallow minimum. Compared to the gas phase, the second crossing point of the two diabatic curves is removed in ambient water, and the ionic surface is the lower in energy at all separations beyond the first crossing (see Figure 3). The main differences between the gas phase and solution phase occur, as expected, due to changes in the pure ionic curve. To see this effect in the water systems, we calculated the set of solvated ionic curves. To do this, we use the ionic potential
Westacott et al.
Figure 4. Effect of supercritical water density on the pure ionic state of tert-butyl chloride (673 K) at 0.29 g cm-3 (dotted line), 0.087 g cm-3 (dashed line), and 0.0435 g cm-3 (dot-dashed line). The ionic state in ambient water (solid line) and the covalent curve (solid line with filled symbols) are shown for comparison. The results are approximate for r e 2.8 Å.
and the adiabatic ground-state simulations in eq 10. This assumes that the adiabatic ground state solvent distribution is a good representation for the pure ionic state. This is not the case when the ground state is predominantly covalent, i.e., in the gas phase or for the solvated tert-butyl chloride molecule (r e 2 Å), but for separations at which the degree of ionic character is higher (r > 2.8 Å), this is satisfactory. The result is shown in Figure 4. We show the full curve, noting the inaccuracy for r e 2.8 Å. For r e 2.5 Å, the repulsive solute potential dominates. Hence, the region affected by the expedient used here to approximate the purely ionic curve influences the result significantly only over a very small region of separation. As solvent density, and therefore dielectric constant, increases, the stabilization afforded to the ionic state by the solvent increases. This can be seen clearly in the figure, where the ionic curves have been set according to eq 12, and the covalent curve is shown for reference. It is clear that increasing solvent density has a dramatic effect on the ionic curve. Even at the lowest solvent density investigated here (0.0435 g cm-3), the effect of the solvent moves the ionic curve to lower energy, and only one crossing point occurs with the covalent curve (cf. two crossings in the gas phase). This indicates that ionic dissociation products are more stable (cf. covalent for the gas phase). Further increase in solvent density moves the ionic curve to lower energy and the minimum to larger r. The single crossing point with the covalent curve is moved to lower energy and smaller separation distances. It is this crossing point at high density that gives rise to the transition barrier at r ) 2.35 Å and the contact ion pair minimum for ambient water (Figure 3). Also, as solvent density increases, the minimum for the ionic curve becomes broader, although this broadening has very little effect on the adiabatic ground state. For the supercritical water cases (Figure 5), a distinct change in shape of the PMF is observed as the solvent density decreases. In particular, the energy of the global minimum for the solvated dipolar t-BuCl becomes more positive, by ∼34 kcal mol-1, over the range of systems studied. The energy of the transition state for the highest supercritical density (0.29 g cm-3) is ∼16 kcal mol-1 higher than that for the ambient case, and consequently, the activation energy increases by 22 kcal mol-1, and the separation distance is increased slightly to about 2.5 Å. A contact ion pair minimum which is significantly shallower than that
Ionic and Radical Molecular Dissociation Pathways
J. Phys. Chem. B, Vol. 105, No. 28, 2001 6617
Figure 5. Adiabatic reaction free energy surfaces for tert-butyl chloride in supercritical water (673 K). The key is the same as that for Figure 4. Asymptotic limits (W(∞)) shown adjacent to each curve on right of figure.
Figure 6. Solute dipole moment along the adiabatic reaction free energy surfaces of tert-butyl chloride in supercritical water (673 K) and ambient water. The key is the same as that for Figure 4.
for the ambient case is observed at ∼2.8 Å. The two lower density cases follow a similar trend, having increasingly shallower wells for the tert-butyl chloride molecule. Of note, at the lowest two densities, there is no local minimum on the surface corresponding to the contact ion pair, although there is stabilization due to the solvent relative to the gas phase. Hence, for these low-density cases, there is no transition state at ∼2.5 Å, but a distinct inflection occurs on these surfaces at a point that corresponds to the transition state for the high density cases. Thus the conventional stepwise mechanism noted in the Introduction does not apply at low density. It is interesting to note that the PMFs for the supercritical water systems are similar to those obtained when ab initio or semiempirical methods with continuum solvent are used to calculate the PMF for ambient water.33,45 Thus, local solvation around the transition state even in ambient water is not adequately represented by the continuum models used with these methods or insufficient water molecules are included when a cluster model is used. Extrapolation of the results presented in Figures 4 and 5 indicates that a change from ionic to covalent character of the adiabatic ground state at large separations would occur at supercritical water densities of ∼0.03 g cm-3. This compares with 0.08 g cm-3 calculated using the simple Born model. The effect of decreasing solvent density on solute polarity can be observed directly from the probability amplitudes for both diabatic states, yielding a solute dipole moment. These were calculated using eqs 2 and 3 and averaged over every time step of the simulation. In Figure 6, the molecular dipole moment is plotted along the reaction coordinate for the four solvent systems. It is clear that the solute changes from being mostly covalent at small separations to completely ionic at large separations. In ambient water, for separations close to the transition state, the change in ionic character is evidently extremely rapid. Between separations of 2 and 2.5 Å, the degree of ionic character changes from 0.2 to 0.95, and at the transition state separation (2.35 Å), it has a value of 0.74. Compared to the ambient water system, the transition from mostly covalent character to completely ionic character is increasingly delayed in distance for the supercritical systems. The rate of change also decreases with solvent density. At the transition state (or inflection on the free energy surface for the two lowest density cases), the degree of ionic character is much less than that for the ambient water case. For the 0.29 g cm-3 case, the degree of ionic character is very close to 0.65 at the transition state.
At the point of inflection on the two lower surfaces, the degree of ionic character is 0.6 (0.087 g cm-3) and 0.54 (0.0435 g cm-3). Although the adiabatic ground state is completely ionic in nature at large separations for all solvent densities studied here, the energy gap between the pure diabatic states at low solvent density becomes very small at the dissociation plateau. In this region, the reaction is likely to proceed with some probability on both the ground and excited state surfaces. It is therefore interesting to consider at what density covalent products become a significant fraction of the total. Through consideration of the free energy difference between the two pure states at the asymptotic limit, it is possible to estimate the relative populations of the two states and therefore the proportion of products. Using this approach, the percentage of products that are covalent are insignificant in ambient water and at 0.29 g cm-3, whereas at lower supercritical water solvent densities, we find the fraction is 3 × 10-4% (0.087 g cm-3) and 22.4% (0.0435 g cm-3) at 673 K. At higher temperatures, it would be expected that covalent products would become increasingly more significant at higher bulk solvent densities. To illustrate local density augmentation here, we have calculated the average local solvent density around the tertbutyl and chlorine species from the molecular dynamics trajectories at important points along the reaction coordinate. These points correspond to the equilibrium molecule, the transition state, the contact ion pair and separated ions for the four solvent densities. The local solvent density around the chlorine moiety is shown in parts a-c of Figure 7; in part d, the tert-butyl cation is also included. These figures show the local density of solvent and at large separations should approach the bulk density. Figure 7a clearly shows that there is no significant density augmentation above the bulk value around the tert-butyl chloride molecule at any of the four densities. This is a reflection of the relatively low dipole moment of the solute (∼2 D). However, as the molecule is dissociated, solvent molecules are attracted to the solute particles. At the transition state (Figure 7b), there is no increase in solvent density for the ambient case, but there is considerable enhancement in the supercritical cases, by a factor of 2.5 at 0.29 g cm-3, a factor of 7 at 0.087 g cm-3, and more than 1 order of magnitude at 0.0435 g cm-3. This shows that solvent density augmentation is responsible for stabilizing the transition state in supercritical water, and the degree of enhancement increases dramatically
6618 J. Phys. Chem. B, Vol. 105, No. 28, 2001
Westacott et al.
Figure 7. Radial solvent densities around the chlorine species of tert-butyl chloride in supercritical water at important points on the reaction free energy surface: (a) upper left: the equilibrium molecule, (b) upper right: the transition state (or point of inflection for the two lowest density cases), (c) lower left: the contact ion pair (or equivalent separation distance for the two lowest density cases), and (d) lower right: the separated ions (tert-butyl cation is indicated by 9). The key is the same as that for Figure 4.
with decreasing solvent density. For the contact ion pair (or equivalent separations at the two lowest densities), a similar but magnified effect is observed as seen in Figure 7c. For the case of the chloride ion of the contact pair in ambient water, the plot looks very similar to that observed for the chloride ion by Yamaguchi and Soper46 using neutron diffraction, with significant density enhancement in three distinct solvation shells at 3.4, 5, and 7.8 Å, with densities of 2.7, 1.3, and 1.1 g cm-3, respectively. For the supercritical solvents, only one solvation shell was observed, with a local density of ∼1.0 g cm-3 for both the 0.29 and 0.087 g cm-3 cases, and 0.7 g cm-3 for the 0.0435 g cm-3 bulk densitysan enhancement of more than 15 times. For the isolated ions, the local densities around both the chloride and tert-butyl ions are shown in Figure 7d. Starting with the chloride ion, at 0.29 g cm-3, the first solvation sphere, at 3.4 Å, has a local density of 0.92 g cm-3. A small second peak, at 6 Å, occurs with a local density of ∼0.4 g cm-3. For the two lower bulk density cases, a single distinct solvation sphere is observed, with local densities of 0.82 and 0.6 g cm-3, respectively. For the tert-butyl cation, there is only one peak for each bulk density, and the local density enhancement is approximately 50% of that for the chloride ion. This is a reflection of the much larger size of the tert-butyl cation relative to the chloride ion.
Conclusions In this paper, we have considered the effect of supercritical water solvent density on the dissociation step of a model SN1 reaction. Using an empirical valence bond model, we have investigated the competition between ionic and covalent reaction pathways for tert-butyl chloride. Molecular dynamics simulations with free energy perturbation have been used to construct the reaction free energy surface for the dissociation of tert-butyl chloride in three supercritical water densities and in ambient water for comparison. The results of the simulations indicate that local density enhancement provides the solvation energy required to stabilize the ionic state at very low supercritical water solvent densities. At very low solvent densities, the local minimum on the free energy surface corresponding to the contact ion pair intermediate is lost. These results can be explained by changes in the solvation of the pure ionic state, which is significantly lowered in energy as solvent density, and therefore dielectric constant, is increased. These results also indicate that the crossover to the covalent reaction pathway would occur at the remarkably low density of about 0.03 g cm-3, significantly lower than that indicated by simple continuum models of the solvent. However, at low solvent density, the free energy difference between the pure ionic and covalent states is small, and the reaction is likely to proceed on both the ground state and first excited state surfaces. We estimate that radical products
Ionic and Radical Molecular Dissociation Pathways are likely to be formed in significant fraction at supercritical water solvent densities below ∼0.08 g cm-3, with an estimated 22% at 0.0435 g cm-3 and equal amounts at 0.03 g cm-3. Acknowledgment. The authors would like to acknowledge the financial support of the Robert A. Welch Foundation and the U.S. Department of Energy. References and Notes (1) Wood, R. H.; Smith-Magowan, D. In Thermodynamics of Aqueous Systems with Industrial Applications; Newman, S. A., Ed.; American Chemical Society: Washington, DC, 1980; Vol. 133, pp 569-581. (2) Mesmer, R. E.; Sweeton, F. H.; Hitch, B. F.; Baes, C. F. In High Temperature High Pressure Electrochemistry in Aqueous Solutions; Jones, D. d. G., Staehle, R. W., Eds.; National Association of Corrosion Engineers: Houston, TX, 1976; pp 365-374. (3) Marshall, W. L.; Frantz, J. D. In Hydrothermal Experimental Techniques; Ulmer, G. C., Barnes, H. L., Eds.; Wiley: New York, 1987; Chapter 11. (4) Flarsheim, W. M.; Bard, A. J.; Johnston, K. P. J. Phys. Chem. 1989, 93, 4234. (5) Johnston, K. P.; Balbuena, P. B.; Xiang, T.; Rossky, P. J. Simulation and Spectroscopy of Solvation in Water from Ambient to Supercritical Conditions. In InnoVations in Supercritical Fluids; Hutchenson, K. W., Foster, N. R., Eds.; American Chemical Society: Washington, DC, 1995; Vol. 608, pp 77-92. (6) Kajimoto, O. Chem. ReV. 1999, 99, 355 and other reviews in the same issue. (7) Antal, M. J.; Brittain, A.; DeAlmeida, C.; Ramayya, S.; Roy, J. C. In Supercritical Fluids. Chemical Engineering Principles and Applications; ACS Symposium Series 329; American Chemical Society, Washington D.C., 1987; Chapter 7, p 77. (8) Townsend, S. H.; Abraham, M. A.; Huppert, G. L.; Klein, M. T.; Paspek, S. C. Ind. Eng. Chem. Res. 1988, 27, 143. (9) Huppert, G. L.; Wu, B. C.; Townsend, S. C.; Klein, M. T.; Paspek, S. C. Ind. Eng. Chem. Res. 1989, 28, 161. (10) Salvatierra, D.; Taylor, J. D.; Marrone, P. A.; Tester J. W. Ind. Eng. Chem. Res. 1999, 38, 4169. (11) Born, M. Z. Phys. 1920, 1, 45. (12) Johnston, K. P.; Bennett, G. E.; Balbuena, P. B.; Rossky, P. J. J. Am. Chem. Soc. 1996, 118, 6746. (13) Tucker, S. C. Chem. ReV. 1999, 99, 391. (14) Song, W.; Biswas, R.; Maroncelli, M. J. Phys. Chem. A. 2000, 104, 6924. (15) Balbuena, P. B.; Johnston, K. P.; Rossky, P. J.; Hyun, J. K. J. Phys. Chem. B 1998, 102, 3806. (16) Rossky, P. J.; Johnston, K. P.,Chemistry in Supercritical Water: Insights from Theory and Simulation. In Steam, Water and Hydrothermal
J. Phys. Chem. B, Vol. 105, No. 28, 2001 6619 Systems: Physics and Chemistry Meeting the Needs of Industry; Tremaine, P. R., Hill, P. G., Irish, D. E., Balakrishnam, P. V., Eds.; NRC Research Press: Ottawa, Canada, 2000; pp 391-401. (17) Chandler, K.; Liotta, C. L.; Eckert, C. L.; Schiraldi, D. AIChE J. 1998, 44, 2080. (18) Webley, P. A.; Tester, J. W. Energy Fuels 1991, 5, 411. (19) Brock, E. E.; Oshima, Y.; Savage, P. E.; Barker, J. R. J. Phys. Chem. 1996, 100, 15834-15842. (20) Luo, H.; Tucker, S. C. J. Phys. Chem. 1996, 100, 11165. (21) Westacott, R. E.; Johnston, K. P.; Rossky, P. J. J. Am. Chem. Soc. 2001, 123, 1006. (22) Balbuena, P. B.; Johnston, K. P.; Rossky, P. J. J. Phys. Chem. 1995, 89, 1554. (23) Marrone, P. A.; Arias, T. A.; Peters, W. A.; Tester, J. W. J. Phys. Chem. A 1998, 102, 7016. (24) March, J. AdVanced Organic Chemistry: Reactions, Mechanism and Structure, 3rd ed.; Wiley: New York, 1985. (25) Catala´n, J.; Diaz, C.; Garcia-Blanco, F. J. Org. Chem. 1999, 64, 6512. (26) Jorgensen, W. L.; Buckner, J. K.; Huston, S. E.; Rossky, P. J. J. Am. Chem. Soc. 1987, 109, 1891. (27) Ford, G. P.; Wang, B. J. Am. Chem. Soc. 1992, 114, 10563. (28) Kim, H. J.; Hynes, J. T. J. Am. Chem. Soc. 1992, 114, 10508. (29) Kim, H. J.; Hynes, J. T. J. Am. Chem. Soc. 1992, 114, 10528. (30) Mathis, J. R.; Kim, H. J.; Hynes, J. T. J. Am. Chem. Soc. 1993, 115, 8248. (31) Keirstead, W. P.; Wilson, K. R.; Hynes, J. T. J. Chem. Phys. 1991, 95, 5256. (32) Hartsough, D. S.; Merz, Jr., K. M. J. Phys. Chem. 1995, 99, 384. (33) Takahashi, O.; Sawahala, H.; Ogawa, Y.; Kikuchi, O. J. Mol. Struct. (THEOCHEM) 1997, 393, 141. (34) Abraham, M. J. Chem. Soc., Perkin Trans. 2 1973, 2, 1893. (35) Ogg, R. A.; Polanyi, M. Trans. Faraday Soc. 1935, 31, 604. (36) Baughan, E. C.; Evans, M. G.; Polanyi, M. Trans. Faraday Soc. 1941, 37, 377. (37) Pross, A.; Shaik, S. S. Acc. Chem. Res. 1983, 16, 363. (38) Warshel, A.; Weiss, R. M. J. Am. Chem. Soc. 1980, 102, 6218. (39) Zwanzig, R. W. J. Chem. Phys. 1954, 22, 1420. (40) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (41) Balbuena, P. B.; Johnston, K. P.; Rossky, P. J. J. Phys. Chem. 1996, 100, 2706. (42) Guissani, Y.; Guillot, B. J. Chem. Phys. 1993, 98, 8221. (43) Neumann, M. Computer Simulation of Water and the Dielectric Equation of State. In Physical Chemistry of Aqueous Systems; White, H. J., Sengers, J. V., Neumann, D. B., Bellows, J. C., Eds.; Begell House, Inc.: New York, 1995. (44) Winstein, S.; Fainburg, A. H. J. Am. Chem. Soc. 1957, 79, 5937. (45) Okuna, Y. P. J. Phys. Chem. A 1999, 103, 190. (46) Yamaguchi, T.; Soper, A. K. J. Chem. Phys. 1999, 110, 3529.
