5196
J. Phys. Chem. 1995, 99, 5196-5205
Temperature and Density Effects on an s N 2 Reaction in Supercritical Water Lewis W. Flanagin,?Perla B. Balbuena,? Keith P. Johnston,*g? and Peter J. Rossky*3* The University of Texas at Austin, Austin, Texas 78712 Received: December 6, 1994; In Final Form: January 23, 1995@
Molecular dynamics computer simulation is used to relate the thermodynamic properties along the reaction coordinate to microscopic solvation for the sN2 reaction of the chloride ion with methyl chloride as a function of temperature and density. Extreme conditions (e.g. reduced densities of 0.05 and 0.3 for reduced temperatures of 1.0 and 1.3, respectively) are found to be necessary to remove half of the water molecules in the first solvation shell about the chloride ion. As the temperature is increased and density decreased, the number of Cl--water hydrogen bonds decays faster than the coordination number. By analogy to adsorption phenomena, augmentation in the local solvent density relative to the bulk (clustering) is interpreted in three regions corresponding to gas, near-critical, and liquid-like densities: The lifetime of a water molecule in the first coordination sphere is found to be about 4 times shorter than under ambient conditions; the lifetime of a hydrogen bond between C1- and water decreases by a comparable factor of about 6. The values of AA,AE, and -TAS associated with conversion from the reactant state to the transition state are explained in terms of the variations in the average coordination numbers and hydrogen bonding. The combined effects of changes in temperature and solvation lead to an increase in the rate constant by 9-12 orders of magnitude under supercritical conditions compared to ambient conditions.
1. Introduction
Supercritical water (SCW) is of interest for hydrothermal breeding of crystals, as a benign solvent for reaction and separation processes, and in hydrothermal oxidation of organic molecules, also called supercritical water oxidation, for waste destruction.i.2 SCW exhibits a unique blend of solvent properties, which can be tuned over a wide range with temperature and density. In essence, it can be used to bridge the fields of aqueous and nonaqueous chemistry over a continuum. Even when the dielectric constant is low, as for nonaqueous solvents, water does not surrender its ability to solvate ions, e.g., when the temperature is increased from 25 “C to the critical temperature, 374 “C, and the bulk density is decreased by a factor of 6.3-4 Despite this large decrease in bulk density (er = O S ) , which lowers the dielectric constant E to 2.65,5 the coordination number about a chloride ion changes little, according to molecular dynamics s i m ~ l a t i o n . ~This ~ ~ .pronounced ~ increase in the local versus bulk density of water in a region near the solute is sometimes termed solvent-solute clustering or microheterogeneous solvation; for ionic solutes, the term electrostriction is also used. This electrostriction can have a large effect on the equilibrium and rate constants for chemical reactions, as shown in theoretical studies of the sN2 substitution reaction of chloride and methyl chloride3s8and of the ion pairing of Na+ and Cl-.437 To date, electrostriction has been prevalent in chloride SCW solutions for the densities studied by simulation. Our objective is to raise temperature and lower density to study the degree of loss of this electrostriction and to examine how this influences the rate of a chemical reaction. Experimentation has elucidated little about the structure of SCW solutions, as spectroscopic studies are rare in this challenging environment. Ion pairing has been characterized in nitrate solutions with Raman spectroscopy9and in 2-naphtholate solution with UV-vis spectroscopy.1° Furthermore, UV-vis
* Author to whom correspondence should be addressed. Department of Chemical Engineering. Department of Chemistry and Biochemistry. @Abstractpublished in Advance ACS Abstracts, March 15, 1995.
0022-3654/95/2099-5 196$09.00/0
solvatochromic measurements indicate that the local densities of water molecules about benzophenone and acetone are augmented above the bulk densities.” In the present study, as has been done previously, l 2 solvent-solute “clustering” is operationally. defined by the existence of a distinct difference between the local and bulk densities. On the basis of spectroscopic studies in both C0213 and water,” solvent-solute clustering has been characterized in three density regions: gasphase solute-solvent clustering, clustering in the near-critical region, and “liquid-like” solvation. The nature of these three regions will be examined in detail in the Results and Discussion section, and an analogy to adsorption on surfaces will be presented. Molecular dynamics computer simulation is an especially important tool for studying SCW solutions, given the small amount of experimental data currently available. In supercritical fluids, computer simulation has evolved rapidly from simple Lennard-Jones mixtures to polar mixtures, to pure water, to ionic aqueous solutions, and recently to ion pairine7 and chemical reactions in water including ions and dipole^.^^^ The chemical reaction studied in this report is the sN2 reaction of methyl chloride with a chloride ion: C1- CH3Cl- ClCH3 C1-. The model used here has been described earlier.14si5It incorporates an asymmetric reaction coordinate defined as rA - rccr - reel-. the difference between the distances from the methyl carbon to each of the chlorine atoms (denoted C1- and Cl’). When the neutral methyl chloride molecule and the negatively charged chloride ion are separated by large distances, the system’s net negative charge of -le is localized on the chloride ion. As the reaction proceeds and the gap between the two solutes closes, the chloride ion transfers charge to the methyl chloride. The specific geometry at rA = 1.441 8, corresponding to a minimum in the gas-phase internal energy is defined as the ion-dipole complex (IDC). In the transition state complex (TSC), defined as the system geometry at rA = 0, both the system’s geometry and charge distribution are completely symmetrical with the methyl group (represented in the model by a single &el5) bearing a partial positive charge
+
0 1995 American Chemical Society
+
sN2 Reaction in Supercritical Water (0.492e) and each chloride bearing a partial negative charge of -0.746e.14.15 The charge dispersal in passing from the reactant state to the transition state accounts for large changes in solvation along the reaction coordinate. Accordingly, the chloride nucleophile becomes less solvated, and the leaving chloride more solvated, as the reaction proceeds to the transition ~ t a t e . ~ . ~ ~ Previously, we have simulated the Sp~2reaction of the chloride ion and methyl chloride to relate the free energy surface to solvation at the molecular level for a reduced temperature T, = 1.0 and reduced density e, = 1.5.338 The free energy barrier at the transition state was nearly as high as in ambient water, corresponding to the small reduction in the coordination number about C1-, which stabilizes the reactants. As density decreases, electrostriction leads to maintenance of a nearly constant firstshell coordination number, even for a reduced density as low as To break up this clustering or microheterogeneous structure, more extreme conditions are required. In light of these observations, the objectives of this study are (1) to characterize when and how water loses its ability to solvate species of varying polarity at high temperatures and low densities, (2) to provide a clear description of solvent-solute clustering and associated lifetimes by comparing new simulation results with previous spectroscopic studies and concepts from adsorption on surfaces,16 and (3) to understand the relationship between solvation at the molecular level and the free energy surface (including the energy and entropy) for this sN2 reaction at low densities. We will emphasize the combination of both the typically aqueous and nonaqueous characteristics of SCW. The system under study is particularly rich, as there are three species of varying polarity including a dipole and an ion, and the solvent forms hydrogen bonds with both the solutes and itself. Thus, we can examine how hydrogen bonding and clustering vary as a function of the charge distribution on the solute species, as these charges change from the reactant state to the IDC to the TSC for a wide range in temperature and density. In the following sections, we describe the methods used in this study and then address the three primary objectives in the order described above.
2. Computational Method Previous reports3.*describe in detail the method employed to compute the free energy of the reaction. Ab-initio calculat i o n ~ ' of ~ ~the ' ~gas-phase reaction provided parameters for the pairwise additive Coulombic and Lennard-Jones intermolecular potential functions acting between the various sites used for both solute-solvent and solvent-solvent interactions. The SPC17 model for water was used; the critical values of temperature and density calculated with this model are 587 K and 0.27 g/cm3, respectively.18 By definition, spans from infinity (separated reactants or products) to 0 (transition state complex). Via a functional form fitted to ab-initio results (see Table 3 in ref 14 and Table 2 in ref 15), all parameters vary smoothly along the reaction coordinate. The incremental change in free energy at each step along the reaction coordinate is calculated by applying statistical mechanical perturbation theory: l9
Here 6A is the change in free energy when the system is drA, U,, is the water-solute perturbed from rA to rA configurational energy, and the brackets indicate an average taken for the solution at fixed r A . The corresponding potential of mean force (PMF) is given by
+
J. Phys. Chem., Vol. 99, No. 14, 1995 5197
where AA(r.4) is the integral change of free energy in the solution at rA, which results from the sum of d l the contributions 6A(rA) when the solutes are brought together in small steps drA from infinite separation to the distance rA. The second term, hUgaS(rA), is the corresponding change in the gas phase. In the calculations reported, we have assumed as a reference point that w(lrAl = 8 A) = 0. This choice is justifiable because the value at ambient conditions is less than 1 kcaUm01.~~ For each point along the reaction coordinate, molecular dynamics simulations were performed for about 20 ps to equilibrate the system. Then the system was perturbed for 1525 ps by drA = f0.125 A to obtain the average in eq 1. The time step was in all cases 1 fs. A total of 484 water molecules was used to solvate this reaction. The size of the simulation box was L x L x 1.5L. For each thermodynamic state condition, the bulk density was obtained by changing the volume of the box, keeping constant the number of water molecules. Thus, L was the appropriate value in angstroms to obtain the desired bulk density in regions far from the solute, for example L = 28.84 A for T, = 1 and er = 1.5.8 In all cases, it was found that the solvent reaches bulk density within 14 A from the solute. As in our earlier ~ o r k , periodic ~,~ boundary conditions were applied and the forces were calculated using the minimum image convention. For each value of rA, the solutes were constrained to be fixed and collinear with the longest axis of the simulation box. The energetic contribution L?&(IA)to the free energy was calculated from the simulations by a method based on finite difference derivatives with respect to temperature20 for all values of ?-A. The entropy change, -Ths(rA), is given by AA - AE.
3. Results and Discussion 3.1. Local Densities. In this section we examine the behavior of the local density, e(r), as a function of r, where r is the radial distance from the center of a solute atom. Here local density is defined as the equilibrium density at a point r, equal to the product of the pair correlation function g ( r ) and the bulk density, e b . Thus, it is the ensemble average of the real number of molecules, dN(r), present in a spherical shell between the radii r - 6 r and r dr divided by the volume of the spherical shell, 6V:
+
The density function eclo(r),which gives the local density of the oxygens in the water molecules about the chloride nucleophile, may be used to substantiate the postulate that solvation decreases with charge dispersal. From simulations performed for 11 test cases of the IDC at four reduced densities (0.05,0.1,0.3,0.5) and at three reduced temperatures (1.0, 1.1, 1.3), two representative conditions were chosen for further study to include the reactants and the TSC: case a with T, = 1.0 (587 K) and e, = 1.5 (0.405 g/cm3) and case b with T, = 1.3 (763.1 K) and e, = 0.5 (0.135 g/cm3). Local densities for reactants, IDC, and TSC at the selected supercritical conditions are compared (Figure 1). The solvent attains bulk density by 14 A, which is beyond the cutoff for these graphs. For both cases and all geometries, the first solvation shell extend! from 3 to 5 A, with a peak occumng in the range 3.25-3.75 A. The height of the peak for this shell represents a maximum in local density, em=.Thus, in passing from reactants to the transition state, emax decreases from 1.80 to 0.66 g/cm3for the case shown
5198 J. Phys. Chem., Vol. 99, No. 14, 1995
Flanagin et al.
1.6
2.0
(a> 0 0
p, = 1.5 (0.405 g / C C )
1.5
0 0
+Reactants
\
L
++TSC
3 1.0 v
I
1.2
+pr = 0.5 (0.135 g/cc)
‘ cn 1
CJ)
2 0.8 L
-ou 0.6
V
Q
0.5
-8- pr =
0.3 (0.081 g/cc)
+pr=
0.1 (0.027 g/cc)
+pr = 0.05 (0.014
v
-0
i
T, = 1.0
1.4
T, = 1.0
g/cc)
0.4 0.2
T, = 1.3
1.2
1
p, = 0.5 (0.1 35 g k c )
0 0
1.0
0
2
\
cn
. 0.8
2 0.6 L
n
5 0.6
v
8 0.4
0
2
0.8
0.4
Q
0.2
0.2
0.0 0
0 2
6
4
8
10
rr A
0
2
4
6
8 1 0 1 2
rr A
Figure 1. Local density gclo(r) for the reactants, IDC, and TSC at (a) T, = 1.0, er = 1.5 and (b) T, = 1.3, e, = 0.5. Note that (a) and (b) have different vertical scales.
Figure 2. Effect of bulk solvent density on the local density eclo(r) for the IDC at: (a) T, = 1.0 and (b) T, = 1.3. Note that (a) and (b) have different vertical scales.
in Figure la, for which @b = 0.405 g/cm3, and from 1.26 to 0.39 g/cm3 for the case in Figure lb, for which @b = 0.135 g/cm3. Furthermore, in case b, even though the temperature is more than 2.5 times greater than in ambient water and the density is 0.135 times smaller, the chloride ion is still highly solvated in the reactant and IDC states, with local densities comparable to the bulk density of ambient water. The reactant state’s local density displays a clearly observable second peak at about 6 8, for both conditions. The corresponding transition state’s local density also peaks near 6 8, for case a (Figure la), whereas for case b (Figure lb), it levels to a plateau following the first shell. A decrease in @clo(r)signifies a loss of electrostriction; thus, the data confirm that solvation diminishes extensively in transforming reactants to the TSC. The ion-dipole complex displays behavior between the two extremes yet similar to that exhibited by the separated reactants, given the larger charge still present on the chloride nucleophile (-0.973e).14,15 A comparison of the plots of local density at the selected conditions reveals that higher temperatures and lower densities within this range inhibit solvation and that the more charged species are less affected than the less charged ones. This is evident in the ratio @‘“x(case b)/emax(casea); one finds for reactants, IDC,and TSC, respectively, 0.70, 0.63, and 0.59.
These observations are in accord with behavior noted in our previous study.3 Additionally, the first solvation shell is broader at the higher temperature and lower density. The kinetic energy imparted by higher temperatures causes the molecules to distribute further into otherwise less favorable radial distances, thereby broadening the solvation shell. The effects of temperature and bulk density on the local density @clo(r)are isolated individually for the IDC in Figure 2. The local density as a function of radial distance is plotted at the model critical temperature and at a reduced temperature of 1.3 for reduced densities in the broad range 0.05-0.5 (0.0135-0.135 g/cm3). Balbuena et d 3included a similar plot at the critical temperature but for reduced densities greater than 0.5. It was observed there that at these higher densities the local density curves are essentially superimposable for the first solvation shell. Although a detailed discussion of the behavior at the critical temperature is included in our previous work,3 it is important to note that substantial local densities persist at this temperature even at much lower bulk densities. However, as can be seen from the data taken at T, = 1.3, the ion solvation does begin to break down as the bulk density is decreased below half of the critical density. It is remarkable that such extreme conditions are required to break down the solvent structure near the solute. Although the precise value of this turning point is
S N Reaction ~ in Supercritical Water
J. Phys. Chem., Vol. 99, No. 14, 1995 5199
system dependent due to the nature of solute-solvent interactions, a significant reduction in solvation between reduced densities of 0.3 and 0.5 has been observed spectroscopically for a variety of supercritical fluids, such as water,ll CF3H,I3 and C O Z . ~ ~ Thus far, we have examined the effects of temperature and density on the solvation about the chloride ion, a highly charged species. We now consider species with less charge, in particular the leaving group Cl'. Previo~sly,~ it was shown that at, e, = 1.5, gcro in the IDC exhibits little structure at T, = 1.0, and the structure disappears at T, = 1.1. The desolvation due to decreases in density is more pronounced for the less charged C1' than for the more ionic C1-. The pair distribution functions obtained for the new set of temperatures and densities (not shown) exhibit the same trends as described above. Desolvation due to temperature increases and density decreases is substantially greater for the less polar species, due to the weaker solute-solvent interactions. 3.2. Coordination Number and Number of ChlorideWater Hydrogen Bonds. Having explored the effects of temperature and bulk density on the local density of solvent around the chloride ion, we now examine solvation in terms of two important characteristics: the coordination number, nc1, and the number of hydrogen bonds, nHB. The coordination number, ncl, is defined by the expression (4) where R " is the distance corresponding to the first minimum after the fist peak in the chloride-water oxygen pair correlation function gclo(r), i.e. enclosing the first solvation shell. The number of ion-water hydrogen bonds is represented by the second quantity, nHB, which here is defined analogously to na, but with gclo replaced by gClH; earlier work3 shows that this geometric definition corresponds closely to an energetic one for this system. For our present calculations, R " was 4.75 8, for nc1 and 3.00 8, for ~ H B . In Figure 3, ncl and ~ H calculated B for the IDC are plotted versus reduced density for three supercritical temperatures. The ncl curves approach the same asymptotic limit corresponding to a full monolayer, at a value corresponding to the value observed in ambient water. At constant temperature, as density decreases in the "near-critical" region (0.5 < e, < 2), ncl and ~ H vary B relatively linle for T, i 1.1. (This region corresponds to the plateau of the Langmuir isotherm that we will discuss below.) In the low-density, "gaslike" region (e, O S ) , both nc1 and ~ H fall B sharply with decreasing density at constant temperature. Thus, for this chloride-water system, the strength of solute-solvent interactions is able to maintain solvent structure in the fist coordination shell down to half of the critical density; beyond this point, highly ordered structure gradually collapses. Temperature effects on ncl and nm in the vicinity of the critical temperature, i.e. from T, = 1.0 to 1.1 (587 to 645.7 K), are small for Qr < 0.5. However, in raising the reduced temperature from 1.0 to 1.3 (587 to 763.1 K for the model), the solvent structure is reduced by almost a factor of 2. At the highest reduced density considered (er = Z), the number of chloride-water hydrogen bonds drops by almost an entire unit (from 4.5 to 3.6) in passing from T, = 1.0 to T, = 1.3. It is interesting to note that at 485 "C and 0.1 1 g/cm3 our simulations give ncl= 4.6,compared to a value of 4.5 obtained using a lattice model that incorporates hydrogen bonding.22 An alternative method for determining ~ H B ,based on an energetic rather than geometric criterion, comes through integrating the fist peak of the water-solute pair energy distribution
8 7
6 0
C
5 4
3 2 I
1
4
3.5 m
I C
3 2.5
I / 2t/ A
0
I
I
I
0.5
1
1.5
2
Pr Figure 3. Effects of solvent density and temperature for the IDC on (a) coordination number ncl, and (b) the number of hydrogen bonds nm.Note that (a) and (b) have different vertical scales.
to the first minimum (see Figure 4). The number of water molecules interacting with the chloride nucleophile in the IDC is given as a function of the energy of pair interactions. The main peak centered near 0 kcdmol is characteristic of the large number of water molecules that are located far away from the solute. As the density decreases at constant temperature, or as the temperature increases at constant density, the distribution of molecules in this central peak range narrows, indicating the decrease in the number of molecules that feel the presence of the solute. The peaks with a large negative energy (about -12 kcal/mol for the IDC) represent the hydrogen bond between C1- and H20. As seen from the degree of deviation among the curves shown in Figure 4a, the amount of hydrogen bonding as determined by the energetic criterion remains nearly constant at the lower temperature (T, = 1.0) over most of the range of densities (0.3 i e,. I2.0) but begins to decline at the lowest density (e, = 0.05). The energy distribution at T, = 1.1 (not shown) is quantitatively and qualitatively very similar to that at the critical temperature. However, at the higher temperature (T, = 1.3), the peak corresponding to the hydrogen bonding has diminished from its value at the lower temperature and continues to decline as density is lowered (Figure 4b). This behavior graphically represents the transformation of hydrogen bonds to less specific interactions that occurs as temperature is raised significantly beyond the critical temperature. Thus, whether using a geometric or energetic basis for measuring the
Flanagin et al.
5200 J. Phys. Chem., Vol. 99, No. 14, 1995
30
2
2
-a -u
flJ 1.5
m
Q
9
20
11
15
b
Q
E
25
1
2
Y-
Y
0
10
0
0
C
0.5
5
10
-3 -
1.5
3 0 Q
2
5
E o 1 Y-
O
0
0
0.5
-5 0 -15
0 -10
-5
0
5
10
Energy, kcal/mol Figure 4. Effects of solvent density on the distribution of pair interaction energies for water molecules with the nucleophile for the IDC at (a) T, = 1.0 and (b) T, = 1.3.
degree of hydrogen bonding, one observes a large loss of hydrogen bonding at T, = 1.3 (763.1 K). In addition to calculating the coordination number for the first shell, we found it instructive to consider effects beyond the first coordination shell for qualitative comparison with spectroscopic measurements. In particular, solvatochromic shifts for UV-visible probes” are influenced by dipole-dipole, dipole-induced dipole, and hydrogen-bonding interactions, which one expects are governed primarily by the first two solvent shells. As noted earlier, these probes indicate three regimes of solvation. The excess number, N x , of solvent molecules through the second coordination shell is given by
where ncl(R2) is the two-shell coordination number. Here, RZ is the distance (7.05 A) corresponding to the minimum after the second peak in gclo(r),and Rcl- is the radius of the chloride ion for the IDC; we use the value RCI-= 2.34 A. This excess number of solvent molecules, which is similar to the integral of the function @local - 4t,12.21-23,24 is more pronounced in supercritical fluids than in liquid solvents and is often used to characterize solvent-solute clustering.
0.5
1.5
1
2
Pr Figure 5. Effects of solvent density and temperature on the total excess number N, of solvent molecules in the first and second coordination shells (eq 5 ) for the IDC.N, is given in (b); ncl through the second shell is shown in (a) for reference.
The variation of ncl(R2) and Nx with solvent density, depicted in Figure 5 , may be compared qualitatively to spectroscopic measurements of local density augmentation. Figure 5a shows that ncl(R2) initially rises quickly with density until it reaches a middensity region (0.3 5 er 5 O S ), where at TI = 1.0 it begins to level off; for T, 2 1.1 it continues to rise, but more gradually at higher density. As seen from Figure 5b, Nx passes through three characteristic regions defined earlier: the gas-phase region, the near-critical middensity region, and the high density region with “liquid-like” solvation. In the low-density, gasphase region, nC1 and N x increase rapidly. In the middensity region (0.3 5 el I OS), ncl begins to level off, and Nx reaches a maximum. In the high-density, liquid-like region, the slope of ncl increases, especially at T, = 1. This increase would be more apparent at even higher densities.13 In this region, N x decreases steadily, finally attaining a value of 0 at about Qr = 2 for the temperatures under consideration. The maximum value of N x decreases with temperature from 14.3 at T, = 1.0 to 10.9 at T, = 1.1, and then to 7.7 at T, = 1.3. Thus, for this system, the excess number of solvent molecules is greatest for lower, near-critical temperatures in the region between gas-phase and liquid-phase behavior. To understand these results for local solvation or solventsolute clustering, we use an analogy to adsorption phenomena.16,25When a fluid is physically adsorbed on an open surface or in a confined medium (pore) by electrostatic and/or Lennard1~13325
S N Reaction ~ in Supercritical Water Jones forces, the following notions are helpful: “average densities”, defined as the mean number of molecules adsorbed per unit volume and also referred to as the total adsorption per unit volume,z6and “excess densities”, defined as the difference between the average density (of molecules adsorbed) and the bulk density (of molecules not adsorbed). When plotted as a function of bulk density, the excess densities always pass through a single maximum and then decrease as bulk density is raised. It is found that the maximum is reached when the first layer is completed. For a given system, the position of the maximum depends on the size of the confined region; that is, lower bulk densities are needed to saturate a smaller region. In our study, the solute is analogous to the pore or open surface. The “confined region” must be assigned a definition, e.g. the first or first and second coordination shell(s). The adsorption isotherm (average density in the pore versus bulk density) corresponds in our case exactly to the behavior of the actual number of molecules in a sphere of a given radius surrounding the solute versus the bulk density. Since the average density or ncl is an absolute property, no maximum is found, and a plateau is reached when the first layer, or the first coordination shell, is completed. The shape of the isotherm from zero bulk density to the density for which the first layer is nearly filled depends on the relative strength of the solutesolvent interactions. This dependence gives rise to different types of isotherms.27 The Langmuir type, which corresponds to attractive surface-fluid interactions, is appropriate for describing the first coordination shell in our case of the strongly attractive ion-water system. The property N, corresponds to the excess density and is found to go through a maximum, as expected. The position of the maximum is not at the critical density, but depends upon the strength of solute-solvent interaction~.’~-~~,~~ With adsorption phenomena as a guideline, it is interesting to see how two different approaches to the same problem, spectroscopy and computer simulation, complement each other. In spectroscopic experiments such as solvatochromic measurements and electron paramagnetic resonance spectro~copy,~~ the range of the cybotactic region is not known a priori. The results of computer simulation allow one to determine on a more solid basis the specifics of structural organization in the fluid. By comparing the solvatochromic measurements with the structural detail available from computer simulation, one can infer the extent of the active solvation region. All of the trends in the simulation results agree qualitatively with spectroscopic studies of microscopic solvation for solvents such as C0213 and water.” A layered isotherm is obtained for sufficiently low reduced temperatures, as in adsorption phenomena. The three regions of solvent-solute clustering appear in each. For nc1 including two (or even more) shells (Figure 5a), there is a sharp increase in the gas-phase region, a plateau in the midregion, and a further increase in the high-density region as the outer levels become more dense. At higher temperatures (Tr = 1.1 and T, = 1.3), the curves exhibit less sharp transitions from low to high densities in both simulation and spectroscopy. As in adsorption, increasing the temperature causes the relatively sharp layering transitions to become smoothed out. Each of these regions has been observed in spectroscopic studies.13,21-25 Thus, the simulation results suggest that the UV-vis and fluorescence spectroscopic probes are influenced primarily by the first-shell solvation shell and then to a lesser extent by the outer shells. Recently,28as an alternative to the above analysis, a new type of excess property was defined as the difference between the local density of solvent about solute minus that of solvent about
J. Phys. Chem., Vol. 99, No. 14, 1995 5201
solvent. It was claimedz8that this new excess property is more physically meaningful than eq 5 because the new excess property satisfies three necessary conditions that Nx (as defined by eq 5 ) purportedly does not, namely, (1) it goes to 0 when the solute becomes the solvent, (2) it goes to 0 for an ideal gas, and (3) it is more related to solute spectroscopic measurements of local density augmentation. All three points of this claim are invalid. Firstly, if the solute becomes the solvent, there is no reason that the N x of eq 5 should be zero. If solvent-solvent interactions are attractive, then there will be positive values of Nx especially in regions where the fluid is compressible because of the first two peaks in g(r). Secondly, for an ideal gas, the molecules are far apart, and there are no molecules in the first few coordination shells, such that Nx is 0, despite the above claim to the contrary. In fact, N , must be 0 when @b is zero since @b appears in each term (see eqs 4 and 5 ) . Thirdly, the spectroscopic studies do not obviously measure either excess property N,, defined by eq 5 , or the proposed alternative. The solvent causes a change in the spectra according to the number and type of solute-solvent interactions, but the density augmentation of the solvent around another solvent molecule cannot be computed from solute spectra. Moreover, since the solvent is spectroscopically transparent (UV-vis, fluorescence, and EPR), the local density of solvent about solvent is not measurable. Hence, the newly proposed excess property28cannot be determined from experiment. 3.3. Lifetimes of Water Molecules in the First Shell of the Chloride Ion. Until this point, we have discussed only static (time-averaged) values for the coordination number and the number of chloride-water hydrogen bonds. However, other important issues are the lability of the molecules in the first coordination shell and the lifetimes of the H-bonds represented in Figure 3. The characterization of the dynamic behavior of molecules located in a given region can be focused on two alternative aspects: the identity of individual molecules present at any time or the fluctuations in the mean values, regardless of the identity of molecules. Here we will discuss only the first aspect; both will be studied in more detail in a future publication. To study the time evolution of a property of a set of tagged molecules present at a given initial time to, samples were taken every At picoseconds. Here, the interval At was varied between 0.2 and 2 ps. Physically, we follow the time evolution of a set of molecules which are identified as present in a specified region at an initial time to. Then, these tagged molecules are considered again at time t later. If at this t = nAt (n = 1, 2, ...) they are again in the region, they have “survived.” The following time autocorrelation function is defined to provide the survival probability: N
N
n
j = I m=O
+
j= I
where Sj (to m a t ) = 1 if a tagged water molecule j is present in the defined region at time to mAr, and 0 otherwise. The sums are taken over all N water molecules, the indicated average is with respect to initial times to, and the denominator normalizes c(t) so that it is initially unity. As defined, c(t) has no contributions from a molecule for any time after that on which it has been found absent. For the coordination number of the chloride ion, we define the region of interest as a sphere of radius R from the center of the solute; the distance R was taken as 4.45 %, for the reactants and 4.75 %, for the IDC. For the H-bonds, a value of R = 2.95 8, was used in all cases. The results for c(t) were averaged
+
5202 J. Phys. Chem., Vol. 99, No. 14, 1995
Flanagin et al.
I\\
0.8
-+-
0.8
I T, = 1 .O,p, = 0.5
-+e- T, = 1.3, p, = 0.5
n c, v
0.6
+T,
0.6
= 1.3, p, = 0.3
n
w
v
u
0
0.4
0.4 I
4\ \
0
0.5
I
0.2 0
0
8
4
12
16
time, ps
Figure 6. Survival probability of molecules in the first coordination shell (2 ps sampling interval) at T, = 1.O and Qr = 1.5 for the reactants, IDC, and TSC. I
I
I
1
1
1.5
2
0.8
n
0.6
c, 0
U
0.4 0.2
0 0
0.5
time, ps Figure 7. Survival probability of Cl--water H-bonds (0.2 ps sampling interval) at T,= 1.0 and er = 1.5 for the reactants, IDC,and TSC.
over 100-200 initial times to chosen randomly from the simulation data. Figure 6 shows the time decay of the correlation c(t) of molecules in the first coordination shell at er = 1.5 and Tr = 1 for different points along the reaction coordinate: reactants, IDC, and TSC. The sampling time interval At was 2 ps. The reason for this choice was to include contributions from molecules that temporarily leave the shell without returning to the bulk. The absolute quantitative results are somewhat sensitive to the value of At (differing at most by only about a factor of 2), but the relative values are not. For the choice of At = 2 ps, fitting the decay to an exponential of the form c(t) = exp(-r/t), we obtained lifetimes (designated t)of 3.2 ps (reactants), 3.2 ps (IDC), and 1.8 ps (TSC). The longer lifetimes for the reactants and IDC are expected due to the stronger solute-solvent interactions. As shown in Figure 7, a similar trend is observed for the lifetimes of the set of H-bonds between water molecules and the chloride for the same density and temperature. The absolute results for lifetimes in this case were much more sensitive to the choice of A t , but the relative values are meaningful in any
1
1.5
2
2.5
time, ps Figure 8. Effect of solvent density and temperature on survival
probability for the hydrogen bonds for the IDC. case. In this report, we present results corresponding to At = 0.2 ps, which is chosen so as to be sufficiently long to avoid apparent hydrogen bond breaking caused solely by solvent librational motion. The sensitivity of the results to the value of At and its relation to experimental measurements of lifetimes will be addressed in a future publication. The curves from Figure 7 were also fitted by an exponential form. The lifetimes obtained were 0.62 ps for the reactants, 0.61 ps for the IDC, and 0.27 ps for the TSC. Again, the stronger hydrogen bonds for the reactant and IDC states lead to longer lifetimes. The effects of temperature and density on the lifetime of the H-bonds in SCW are observed in Figure 8 for the IDC. Lower densities and higher temperatures decrease the stability of the H-bonds. In ambient water, we obtained an average lifetime of 4 ps (for the same At = 0.2 ps), Le., about 6 times longer than at the highest supercritical density we have studied (er= 1S, not shown). The same analysis for the coordination number (not shown) reveals that under supercritical conditions the lifetime of the water molecules in the first coordination shell remains the same (approximately 3 ps) for a wide range of densities and temperatures (erbetween 0.3 and 1S ; Tr between 1.0 and 1.3), while for ambient water with the same model and choice of At the lifetime is on the order of 10-14 ps. The higher kinetic energy at higher temperatures relative to the interaction strength shortens the lifetimes. 3.4. Free Energy Calculations. Figures 9-12 relate the calculated changes in free energy in solution of the system and the concomitant changes in internal energy and entropy. The uncertainty in the calculated changes in energy progressively increases from reactants to the TSC because the simulation time required for the solvent to rearrange around a charged species increases as the charge decreases. We note here that our experience is that systems dominated by electrostatic forces have larger free energies of solvation but can be simulated with the same precision in less time than those where the main interactions are of the van der Waals type. The reason might be that the strong Coulombic interactions rearrange solvent molecules faster than in the van der Waals case, where packing and translational diffusion are the predominant mechanisms for the solvent rearrangement^.^^ Thus, we find that, for these reasons, there is more uncertainty in the exact value for the energy barrier as one approaches r A = 0 than for the more charged reactants to the IDC. The changes in internal energy are calculated by
S N Reaction ~ in Supercritical Water
J. Phys. Chem., Vol. 99, No. 14, 1995 5203 I
10
-E"m
5
3 a
15
0 Y
T, = 1.O,p, = 0.5
10
0 1
2
3
4
5
6
7
8
0
Figure 9. Effects of solvent density and temperature on changes in free energy of solvation as a function of the reaction coordinate r A .
D
' I
1
'
q
1 '
'
I
"
' ' 1
' 1 1 ,
'
1
I
=1.0,pr=1.5
- t T
= 1.0, p ~ 0 . 5
+T
= 1.3, p = 0.5
6
2
E b
2
0
5
&
o
a
m"
e 1
2
3
4
5
6
7
4
5
6
7
8
o -2 -4
1-
-6 0
3
1
4
I
a
phase
1
Figure 11. Effects of solvent density and temperature on changes in solvent internal energy in solution as a function of the reaction coordinate rA.
''1
+T
*gas
Y
z e
I
I
5
0
LL-
I
I
+
15
0
20
I
I
20
\
25 %
I
8
0
1
2
3
+Tr
= 1 .O,p, = 0.5
+TI
= 1.3, p, = 0.5
4
5
6
7
8
r ,A
A Figure 10. Effects of solvent density and temperature on potential of mean force (total free energy change) as a function of the reaction coordinate rA.
(formally) differentiating the free energies; therefore, the statistical error in internal energy is even greater than that of the free energy.20 Overall, then, we find that these plots are more useful for discussing general trends in behavior rather than exact values of these quantities, particularly AE and -TAX The change in free energy of the solvent system along the reaction coordinate is depicted in Figure 9. Due to the reduction in cnlvatinn nf Pl- alnnu the r e n r t i n n rnnrdinate a c rharue ir
distributed more equally between the solutes, the free energy change is, as expected, positive, indicating a barrier for the transformation of reactants to TSC. As solvent molecules are added, the reactant state is stabilized more than the transition state because of stronger interactions with the solvent. The barrier at the transition state is insensitive with respect to density in the region er = 0.5-1.5 (0.135-0.405 g/cm3) but is affected by the temperature. As the reduced temperature is increased from 1.0 to 1.3, the barrier decreases from 19 to about 16 kcaY mol. The reduction in the barrier is a consequence of the
Figure 12. Effects of solvent density and temperature on changes in entropy in solution as a function of the reaction coordinate r A .
reduction in the relative strength of the solvent's interaction with solute and solvent. The effect of temperature and density on the formation of the IDC is seen better by examining the potential of mean force (PMF) curve (Figure 10). The PMF previously has been calculated for the gas phase, in ambient water, and in SCW.8J4,15 A deep well in the PMF coincides with the formation of the IDC and has been evident in SCW and the gas phase but not in ambient water. The PMF curve depicted in Figure 10 shows that the minimum in the well, corresponding to the IDC, is virtually identical for the two conditions at the critical temperature but drops by about 2 or 3 kcdmol for the higher temperature. The PMF may be used also to obtain values for the equilibrium constant and hence the rate constant via transition state theory. The equilibrium constant, K,is dominated by AW = - W",the PMF activation barrier for the reaction, so that one has the following proportionality:
wsc
5204 J. Phys. Chem., Vol. 99, No. 14, 1995
Flanagin et al.
SCW: T = 763 K, p = 0.135 g k c
(7)
15.9 kcal/mol 17.5 kcal/mol -T ASscw= - 1.6 kcal/mol AA,,,
which can be used to estimate differences in rates. Because of this exponential relationship, small changes in the PMF tremendously affect the extent of reaction. For e, = 0.5, the ratio K / K A w , where the subscript AW denotes ambient conditions, is, based on AWIRT alone, about 4 x lo9 for TI = 1.0 and about 2 x lo1*for TI = 1.3. The corresponding values for AW differ by about 2 kcaymol: 26.07 and 24.13 kcaymol, respectively (much less than the difference in AA at r A = 0 since both the TSC and IDC are affected). Thus, the equilibrium constant at T, = 1.3 (763.1 K) is about 600 times greater than that at T, = 1.0 (587 K). In transition state theory, since k is proportional to K , the corresponding ratios in the rate constant MAWare approximately lo9 and 10l2, respectively. As in the case of the free energy, the internal energy monotonically increases as the reaction proceeds to the transition state (Figure 11). The energy barrier accompanying the transformation of reactants to products is about 18 kcdmol for T, = 1.3 and about 23-24 kcallmol for the two cases at TI = 1.0. Again the temperature has a significant impact on the internal energy of the system. Higher temperatures cause the solvent to lose solvating power; thus, less energy is required for the loss of solvation in going from reactants to the TSC. The effect of bulk density in the region near the IDC is that the difference in heat of solvation, EIDC- Ereactants, is greater at the lower density because the IDC, in comparison to the reactants, is less solvated at the lower density. The entropic behavior (Figure 12) at 8, = 0.5 differs significantly from that at e, = 1.5 over much of the range in TA, although the values are similar at the TSC. We must allow the possibility that this is a result of statistical error, although we have no evidence for that at the moment. For e, = 0.5, temperature has little effect on -TAS for r A down to 1 A; however, a significant difference is evident for smaller values of r A including the TSC. At the TSC, -TAS is about -3 to -4 kcallmol for both conditions at the critical temperature (T, = 1.0) but rises to -1.6 kcallmol for the higher temperature (I", = 1.3). Because the magnitude of -TAS is smaller at the higher temperature, increasing the temperature causes the water molecules near the reactants to become more disordered, approaching the same level of disorder as in the TSC. As seen from the radial distribution functions (Figure l), the difference in structure between the reactants and TSC decreases with higher temperatures, thereby signaling a loss in electrostriction about the reactants. This loss of electrostriction, in turn, causes the entropy difference to drop as well. In a separate we have performed free energy calculations for this system using an electrostatic model that treats the solvent as a c ~ n t i n u u m . ~These ~ - ~ ~ results will be reported elsewhere.30 Qualitatively, the preliminary results obtained from this model are very encouraging,considering that no adjustable parameters are used. These results indicate that, quantitatively, the model misses the effects of local solvation; at lower densities, the model significantly underestimates AA. One expects that the use of bulk properties in the continuum model will fail to account for local solvation, which raises AA, and the available results imply that, in fact, the error in the continuum model becomes larger as the excess solvation grows, when the reduced density is lowered from 1.5 to 0.5 (see Figure 5b). The results of the present study are summarized succinctly in Figure 13. The values of AA, AE, and -TAS are given in the horizontal direction at the bottom for ambient water (AW) and at the top for SCW at 763 K and 0.135 g/cm3 (TI= 1.3, e1
=
AEscw =
[Reactants]
[Reactants]
*
-
[TSCI,,,
I
,, AA,,
= 22.7
AE,,
=
-T AS,,
0'30 [TSCIAW
kcal/mol
22.0 kcal/mol = 0.7 kcal/mol
AW: T = 298 K, p = 1.0 glcc Figure 13. Themdynamic cycle for the energetic changes. Vertical paths: AW to SCW. Horizontal paths: reactants to TSC.
= 0.5, for the model solvent). The ratios of first-shell coordination numbers and hydrogen bonds at SCW compared to AW are given in the vertical direction for both the reactants and the TSC. The diagram illustrates that the coordination number about the chloride nucleophile decreases at the supercritical conditions relative to ambient conditions, but only by about 20-25%. Moreover, the amount of hydrogen bonding is significantly lower at the high-temperature conditions, especially for the TSC. The additional thermal energy gained and density lost in changing from ambient to supercritical conditions break almost a third of the H-bonds in the reactant state and two-thirds of the H-bonds for the TSC. This result is consistent with the smaller charges and hence weaker bonds in the TSC. It is instructive to compare Figure 13 with an analogous plot (Figure 12 of ref 3) for TI = 1.0 and eI = 1.5. In the earlier study,3 the coordination numbers RCI changed little from AW to SCW, and likewise AE for the reaction (horizontal paths) did not vary. At the higher temperature and lower density in the new study, the modest loss in nc1 and substantial loss in hydrogen bonding produce the 20% drop in AE at the TSC. The change in -TAS under these conditions is negligible, unlike the earlier study at TI = 1.0, where the change in AA from AW to SCW was entropy driven. As expected, the decrease in AA from AW to SCW is considerably larger (6.8 versus 3.4 kcall mol) at T, = 1.3 and PI = 0.5 than in the earlier study at a lower temperature and a higher density. 4. Conclusions
Previous MD simulation studies demonstrated the great affinity of water for the C1- ion, in that the coordination number changes little from AW to SCW where T, = 1 and e, = 1.5. The new results suggest that rather extreme conditions are required for desolvation. At T, = 1, ncl and the free energy barrier for the sN2 reaction change little as the density is lowered from 0.4 to 0.14 g/cm3. To remove half of the water in the first coordination shell, the bulk density must be reduced to 0.08 g/cm3 at 760 K. Water continues to solvate ions in the first shell even when the bulk dielectric constant is drastically reduced. For a density of 0.135 g/cm3. an increase in T, from 1.O to 1.3 reduces the solvation enough to lower the net barrier by 2 Wmol. The rate constant for the reaction in the range T,
S N Reaction ~ in Supercritical Water = 1.0-1.3 is, however, 9-12 orders of magnitude greater than in ambient water. At T, = 1.3, -TAS changes little from AW to SCW, so that changes in AA mimic those in AE. In contrast, changes in -TAS are larger than those in AE from AW to SCW at Tr = 1.0and er = 1.5. A detailed analysis of the number of water molecules surrounding the chloride ion in the first and second shells provides clear evidence of the persistence of a well-defined solvent structure (clustering) for SCW densities down to half of the critical density, closely paralleling the inferences based on spectroscopic studies. More modest conditions lead to a decay in the number of water-chloride ion H-bonds. The three regions of clustering at gas, near-critical, and liquid-like densities are also apparent in both spectroscopy and simulation. A clear understanding of the phenomena emerges by an analogy to the behavior of total and excess densities in the context of adsorption isotherms. The stability of water-ion clusters is of primary interest in determining their influence in dynamic processes. Here we show that, over a range of densities and temperatures (0.3 < e, < 1.5; 1.0 < T, < 1.3),solute-solvent clusters in SCW are short-lived compared to their ambient counterparts, with an average lifetime of 3 ps versus an average ambient lifetime of 12 ps for the present model. Water molecules forming H-bonds with the chloride ion in SCW yield lifetimes reduced by a comparable amount, about 6 times shorter than those in ambient water.
Acknowledgments are made to the U.S. Army Research Office for University Research Initiative Grant No. DAAL 0392-6-0174and DAAH 04-93-6-0363and to the Separations Research Program at the University of Texas. Partial support of this work by a grant from the R. A. Welch Foundation (F0761)to P.J.R. is also acknowledged. We thank the University of Texas System Center for High Performance Computing and Cray Research Inc. for computational support. Partial support of this work was also provided by the Army Research Office contract number DAAL03-89-C-0038 with the University of Minnesota Army High Performance Computing Research Center (AHPCRC) and the DoD Shared Resource Center at the AHPCRC. This material is based upon work supported under a National Science Foundation Graduate Research Fellowship awarded to L.W.F. References and Notes (1) Shaw, R. W.; Brill, T. B.; Clifford, A. A.; Eckert, C. A,; Franck, E. U. Chem. Eng. News 1991, 69, 26.
J. Phys. Chem., Vol. 99, No. 14, 1995 5205 (2) Tester, J. W.; Holgate, H. R.; Armellini, F. J.; Webley, P. A,; Killilea, W. R.; Hong, G. T.; Bamer, H. E. In ACS Symposium Series; Tedder, D. W., Pohland, F. G., Eds.; American Chemical Society: Washington, DC, 1993; Vol. 518, p 35. (3) Balbuena, P. B.; Johnston, K. P.; Rossky, P. I. J. Phys. Chem. 1995, 99, 1554. (4) Cui, S. T.; Harris, J. G. Chem. Eng. Sci. 1994, 49, 2749. (5) Archer, D. G.; Wang, P. J. Phys. Chem. Ref:Datu 1990, 19, 371. (6) Cummings, P. T.; Cochran, H. D.; Simonson, J. M.; Mesmer, R. E.; Karabomi, S . J. Chem. Phys. 1991, 94, 5606. (7) Gao, J. J. Phys. Chem. 1994, 98, 6049. (8) Balbuena, P. B.; Johnston, K. P.; Rossky, P. J. J. Am. Chem. SOC. 1994, 116, 2689. (9) Spohn, P. D.; Brill, T. B. J. Phys. Chem. 1989, 93, 6224. (10) Xiang, T.; Johnston, K. P. J. Phys. Chem. 1994, 98, 7915. (11) Bennett, G. E.; Johnston, K. P. J. Phys. Chem. 1994, 98, 441. (12) Kim, S.; Johnston, K. P. Ind. Eng. Chem. Res. 1987, 26, 1206. (13) Sun, Y.-P.; Fox, M. A.; Johnston, K. P. J. Am. Chem. SOC. 1992, 114, 1187. (14) Chandrasekhax, J.; Smith, S. F.; Jorgensen, W. L. J. Am. Chem. SOC. 1985, 107, 154. (15) Huston, S. E.; Rossky, P. J.; Zichi, D. A. J. Am. Chem. SOC. 1989, 111, 5680. (16) Kajimoto, 0.;Futakami, M.; Kobayashi, T.; Yamasaki, K. J. Phys. Chem. 1988, 92, 1347. (17) dePablo, J. J.; Prausnitz, J. M.; Strauch, H. J.; Cummings, P. T. J. Chem. Phys. 1990, 93, 7355. (18) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Intermolecular Forces; Reidel: Dordrecht, 1981. (19) Zwanzig, R. W. J. Chem. Phys. 1954, 22, 1420. (20) Fleischman, S. H.; Brooks, C. L., IIJ J. Chem. Phys. 1987,87,3029. (21) Knutson, B. L.; Tomasko, D. L.; Eckert, C. A.; Debenedetti, P. G.; Chialvo, A. A. In ACS Symposium Series; Bright, F. V., McNally, M. E. P., Eds.; American Chemical Society: Washington, DC, 1992; Vol. 488, p 60. (22) Gupta, R. B.; Johnston, K. P. lnd. Eng. Chem. Res. 1994,33,2819. (23) O’Brien, J. A.; Randolph, T. W.; Carlier, C.; Ganapathy, S . AIChE J. 1993, 39, 1061. (24) Tom, J. W.; Debenedetti, P. G . Ind. Eng. Chem. Res. 1993, 32, 2118. (25) Carlier, C.; Randolph, T. W. AIChE J. 1993, 39, 876. (26) Tan, Z.; Gubbins, K. E. J. Phys. Chem. 1990, 94, 6061. (27) Balbuena, P. B.; Gubbins, K. E. Langmuir 1993, 9, 1801. (28) Chialvo, A. A.; Cummings, P. T. AIChE J. 1994, 40,1558. (29) Kollman, P. Chem. Rev. 1993, 93, 2395. (30) Bennett, G. E.; Johnston, K. P.; Rossky, P. J. In preparation. (31) Sharp, K. Columbia University, 1988. (32) Gilson, M. K.; Honig, B. Proteins 1988, 4, 7. (33) Nicholls, A.; Sharp, K. A.; Honig, B. Proteins 1991, 11, 281 JP943232R
